Minimally Modified Gravity fitting Planck data better than $\Lambda$CDM

We study the phenomenology of a class of minimally modified gravity theories called $f(\mathcal{H})$ theories, in which the usual general relativistic Hamiltonian constraint is replaced by a free function of it. After reviewing the construction of the theory and a consistent matter coupling, we analyze the dynamics of cosmology at the levels of both background and perturbations, and present a concrete example of the theory with a $3$-parameter family of the function $f$. Finally, we compare this example model to Planck data as well as some later-time probes, showing that such a realization of $f(\mathcal{H})$ theories fits the data significantly better than the standard $\Lambda$CDM model, in particular by modifying gravity at intermediate redshifts, $z\simeq743$.


I. INTRODUCTION
The ΛCDM model for cosmology depends on 6 parameters [1,2]. Various features of our universe are successfully described by the ΛCDM model. Setting aside the cosmological constant problem and other theoretical issues, it provides us the simplest description of the universe with the cosmic microwave background radiation, the large-scale structure and the accelerated expansion, the last of which we probe for example by supernovae observations. Up to now, this model is still considered the best fitting model to the cosmological data sets.
Despite the success of the ΛCDM model to describe our universe, it faces tensions in explaining a few parameters from various cosmological data sets. This situation motivates us to consider models beyond ΛCDM. The most famous and significant tension lies between the estimation of today's Hubble expansion parameter H 0 from Planck data [3], and its measurement from observations of the local universe, such as those from the Hubble Space Telescope (HST) [4], H0LiCOW [5], Megamaser Cosmology Project (MCP) [6], Carnegie-Chicago Hubble Program (CCHP) Collaboration [7], etc. The significance of this first tension goes up to 5.3σ [5]. This discrepancy can be understood as a tension between late time data, which do not assume any prior cosmological model, and early universe data, which assume a prior model, i.e. the ΛCDM. It comes to light as the early universe and the late time observations (HST) have achieved greater precision in their measurements. At present, there are several approaches to resolve this particular tension, including modified gravity (see [8] for example), or adding new types of matter content (as early dark energy [9] for example).
The number two tension lies between the estimation of the growth of structure S 8 from Planck data and that from the redshift space distortion data. The disagreement has a statistical significance of up to 3.2σ [10]. There are several investigations to address this issue from the perspective of modified gravity theories [11,12].
Two other results could fit the bill as tensions within the ΛCDM model. The Planck 2018 result has reported a preference for a higher lensing amplitude [13,14]. This anomaly has been argued to be related to a new possible tension in the ΛCDM model of cosmology, known as Ω K tension, which is investigated in [15][16][17]. In fact the Planck collaboration also has reported that the Planck data alone prefers a closed universe [18]. In [17], it has been found that the curvature parameter of the universe is in tension with the Planck+late time data and BAO data sets. In their investigation they have found better fit to closed universe than ΛCDM.
The cosmological tensions arising within ΛCDM can be interpreted as the indication that ΛCDM may be only a first approximation to a more accurate theory of cosmology. In this context it is worth investigating the cosmology of modified theories of gravity to address these various tensions. As mentioned, there have been several attempts which rely on modified theories of gravity. In this work, we focus on a certain class of modified theories of gravity in which there are only two local gravitational degrees of freedom propagating [19][20][21], dubbed minimally modified gravity (MMG) theories. These theories break four dimensional diff-invariance keeping the three dimensional diff-invariance so that the standard framework of cosmological studies can be applied. Even though breaking four dimensional diff-invariance leads, in general, to new degrees of freedom in addition to the usual two gravitational modes, MMG theories have Hamiltonians linear in the lapse function so that they come with a Hamiltonian constraint that will reduce the dimension of the physical phase space and leave only two tensorial modes in the theory.
Recently there has been a Hamiltonian construction for a class of MMG theories, dubbed f (H) theories [22], in which the standard Hamiltonian constraint H of general relativity (GR) is modified to a free function f (H). This leads to changing the kinetic structure of the theory from GR in a background-dependent way and provides an opportunity to address the cosmological tensions.
In this article we study the cosmology in f (H) theories, considering a model that we have named "kink" model as it has a kink in the first derivative of the free function f (H). We introduce a consistent coupling to matter which relies on a gauge fixing in the Hamiltonian as it was done in [23,24], and we study linear perturbations so that we can exhibit the no ghost conditions. We then implement the scalar linear perturbation equations in the Boltzmann code called CLASS [25], with covariantly corrected baryon equations of motion [26]. Subsequently a Monte Carlo sampling, using MontePython [27,28], is done against various cosmological data sets. We consider data of Planck 2018 with planck_highl_TTTEEE, planck_lowl_EE, planck_lowl_TT, planck_lensing polarization [18,[29][30][31][32], of HST observations [4] consisting in the single data point of the Hubble constant H 0 = 74.03 +1. 42 −1.42 , of the baryon acoustic oscillation (BAO) from 6dF Galaxy Survey [33] and the Sloan Digital Sky Survey [34,35], and of the joint light curves (JLA) comprised of 740 type Ia supernovae [36]. We refer to all these data sets as Planck2018+HST+BAO+JLA. The chains are then analyzed using the well-suited GetDist package [37]. We find that there is a remarkable improvement with respect to ΛCDM in the likelihood-parameter χ 2 with a difference of ∆χ 2 = 16.6 for the chosen data sets. Although there is only a minimal improvement for the H 0 tension, the other tensions that have appeared within ΛCDM have a chance to be addressed by the f (H) theory.
This paper is organized as follows. In section II, we review f (H) theories from their inception, and discuss the particular matter coupling that has to be adopted to avoid the propagation of unwanted degrees of freedom. In section III, we use the Lagrangian introduced in section II to derive the dynamics of the cosmological background, as well as its perturbations. In section IV, we propose a concrete model where the function f depends on three free parameters and has a kink in its first derivative. As we will show subsequently, such a model has a very good fit to cosmological data. Then, in section V we show that the kink model consistently gives a better fitness parameter than ΛCDM, considering both early and late time cosmological data sets together. This comparison is the main result of this paper. Finally, we conclude in section VI with a summary of results and a discussion.

II. MINIMALLY MODIFIED GRAVITY AND f (H) THEORIES
Minimally modified gravity (MMG) theories are modifications of four-dimensional GR with two local gravitational degrees of freedom. A systematic construction of gravitational theories with only (up to) two degrees of freedom has been initiated in [20,38] (see also [24] for a different perspective). The idea consists in renouncing the invariance under four dimensional diffeomorphisms but keeping the three dimensional (spatial) diff-invariance. As generically Lorentz-breaking gravity theories have more than two degrees of freedom, one has to find the conditions for the theory to possess enough constraints that would kill the extra degrees of freedom, which would leave us with (at most) two gravitational modes only. This section is devoted to review these conditions following [22].
In the first part, we will quickly recall the basis of the construction of MMG theories as it was done in [22] where the Hamiltonian point of view was adopted. Then, we will show how f (H) theories naturally emerge as simple but interesting examples of MMG theories. Finally, we will explain how to couple matter to these theories without introducing new degrees of freedom in the gravitational sector, following the construction developed in [20,23].

A. Modifying the phase space
The construction of the f (H) class of MMG theories presented in [22] relies on the Hamiltonian formulation of GR. The idea consists in modifying the phase space of GR, and not directly the Lagrangian, in such a way that the modified theory remains invariant under spatial diffeomorphisms only but still propagates two tensorial degrees of freedom.
Hence, we start with the Arnowitt-Deser-Misner (ADM) parametrization of the metric, in terms of the lapse function N , the shift vector N i and the induced spatial metric γ ij . Hereafter, we will denote by D i the covariant derivative compatible with γ ij , we will lower and raise spatial indices by the metric γ ij and its inverse γ ij , and we will use the notation γ for the determinant of the spatial metric. The phase space is parametrized by the usual ten pairs of conjugate variables, where π ij , π i and π N are momenta, and δ(x − y) is the 3-dimensional δ-distribution.
In GR, N and N i are Lagrange multipliers, thus their conjugate momenta π N and π i vanish, which produces a set of four primary constraints, We are using the notation ≈ for the weak equality in the phase space, i.e. the equality up to constraints. Concerning the momenta conjugate to γ ij , they are easily related to the extrinsic curvature tensor by, with K ≡ K i i being the trace of the extrinsic curvature. As a consequence, one can immediately compute the Hamiltonian of GR which takes the very well-known form, where the Hamiltonian constraint H 0 and the vectorial (momentum) constraints H i are given by, with R being the spatial curvature of the metric γ ij . The conservation under time evolution of the primary constraints (3) leads to the secondary constraints H 0 ≈ 0 and H i ≈ 0, which form together with (3) a set of first class constraints. The Hamiltonian analysis closes here: as we started with 10 pairs of variables (2) and we found 8 first class constraints, we end up with the expected 2 tensorial degrees of freedom of GR.
In [22], one proposed a deformation of the phase space of GR requiring that the modified theory remains invariant under spatial diffeomorphisms only yet still propagates two tensorial degrees of freedom. In this approach, one starts with the same (non-physical) phase space parametrized by the usual ten pairs of conjugate variables (2) and one looks for a total Hamiltonian of the form, where S is a three-dimensional scalar constructed from the variables (γ ij , π ij , R ij , N ) and their spatial derivatives. It was implicitly assumed that neither N nor N i are dynamical variables. Requiring that the theory propagates only two degrees of freedom (or less) leads necessarily to the condition that S is an affine function of N , i.e.
which is, of course, compatible with the fact that N is not dynamical. Hence, H 0def can be viewed as a deformation of the usual Hamiltonian constraint of general relativity whereas V is a new term. However, the functions V and H 0def are not arbitrary and must satisfy extra conditions to ensure that no degrees of freedom other than the two tensor modes propagate [22]. Even though these conditions have not been solved in full generality and rigor, it was shown that any theory whose Hamiltonian is given by (7) with while V is totally free, defines a MMG theory. In other words, (9) is a sufficient condition for the theory defined by the Hamiltonian (7) to propagate (at most) two degrees of freedom. Finding all the functions H 0def (γ ij , π ij , R ij , N, D i ) which satisfy (9) seems a priori to be a highly complicated problem. The reason is that, for the theory to propagate gravitational waves, H 0def must depend on both π ij (which contains time derivative of γ ij ) and three-dimensional curvature terms (which contain gradients of γ ij ), and the Poisson bracket between two such terms has, in general, a very complex expression. Fortunately, this problem admits a solution that we know very well, that is the Hamiltonian constraint H 0 of GR (6) which satisfies, where H 0 [N ] and H i [N i ] are smeared constraints defined, for any function N and any vector field N i , by the integrals, As it is very well-known, the property (10) is intimately linked to the invariance under diffeomorphisms of GR. An immediate consequence is that any deformed Hamiltonian constraint of the form where f is an arbitrary function, is also a solution of (9), because it satisfies, for any N 1 and N 2 . The MMG theories defined by (12) have been dubbed f (H) theories with reference to f (R) theories. However, contrary to f (R) theories, f (H) theories do not propagate a scalar mode in addition to the tensor modes. More precisely, f (H) theories were defined with the additional condition that V = 0 in the Hamiltonian (8), which we will also assume hereafter.
By a Legendre transformation, one can easily compute the corresponding action. Indeed, the equation of motion for γ ij ,γ enables us to relate the momenta π ij to the extrinsic curvature K ij as follows, from which we can implicitly obtain π ij in terms of K ij because, in general, this equation is non-linear in π ij . Nonetheless, one can compute the action which, after a simple calculation, is given by where C is formally obtained by solving the equation Hence, we obtained the gravitational part of the action of f (H) theories. Now we are going to explain how to couple consistently this action to matter.
C. Coupling to matter: problem and solution As was shown first in [23,24], one cannot couple matter minimally at this stage of the construction (except in the trivial case f (H 0 ) = H 0 ) precisely due to the deformation of the Hamiltonian constraint. As long as there is no matter, the deformed Hamiltonian constraint f (H 0 ) satisfies the deformed diffeomorphisms algebra (13) which ensures that f (H 0 ) remains first class. However, if one directly adds the matter part to the gravitational constraints as it is done in GR, one should consider the gravity+matter constraints and one shows that the new deformed Hamiltonian constraintH 0def is, in general, no longer first class because In short, {H 0def [N 1 ],H 0def [N 2 ]} is not weakly vanishing and then the Hamiltonian constraint is downgraded to a second class constraint, which leads to the existence of an extra degree of freedom in the theory. A solution to this problem was suggested in [23]: it consists in adding a gauge fixing term at the level of the vacuum Hamiltonian, that will split the Hamiltonian constraint into a pair of second class constraints, which then allows to couple matter minimally. We will follow this procedure, and give an explicit example for a gauge-fixing in the next subsection. 1

D. Consistent gauge-fixing in the Hamiltonian
In this subsection, we specify a gauge fixing procedure that can be used to couple matter fields consistently. By adding the gauge fixing term thanks to the Lagrange multiplierλ i , to the Hamiltonian of the f (H) theory, we obtain a new form of the gravitational Hamiltonian given by, where we used a simple integration by part in the gauge fixing term. Now, the Hamiltonian equation of motion for γ ij leads to which can be solved with respect to π ij as Here H 0 is viewed as a function of the velocities K ij (and no more of the momenta π ij ), and is obtained from the following implicit equation, that generalizes (17) in the presence of the gauge fixing term. Performing a Legendre transformation, we immediately obtain the Lagrangian density where H 0 is still given by (24). In order to implement the relation (24) directly in the Lagrangian, we consider, instead of (25), the equivalent Lagrangian densitỹ where the Lagrange multiplier λ 1 ensures that the new field C satisfies C = H 0 . At this stage, we can safely introduce matter fields minimally coupled to the metric (1) and then we will use this Lagrangian density in the remainder of the paper.

III. COSMOLOGY: BACKGROUND AND LINEAR PERTURBATIONS
This section aims at studying cosmological properties of f (H) theories whose dynamics is governed by the Lagrangian density (26) supplemented with a matter Lagrangian density L mat of a perfect fluid which is given by the Sorkin-Schutz Lagrangian density [40] (see also [26,[41][42][43] for example). In the first subsection, we introduce some notations, compute the equations of motion for a cosmological (homogeneous and isotropic) background and give preliminary results concerning linear perturbations about this background. This enables us to extract in particular linear stability conditions for the tensor modes and the scalar mode. In the last two subsections, we explain how to implement the equations for the background and the linear perturbations in the Boltzmann code solver in order to solve them numerically and to compare to Planck data and other experimental probes, which we are going to do in the subsequent sections.

A. Cosmology of f (H) theories: notations and preliminary results
As we have just announced in introducing this section, we consider the total Lagrangian density where M P is the Planck mass,L grav the f (H) Lagrangian density (26) and the matter is supposed to be a perfect fluid whose dynamics is governed by the Schutz-Sorkin Lagrangian density L mat , where ϕ is a scalar field, ρ is the energy density, which depends on the number density n defined from the vector J µ and the metric g µν according to (28). A detailed study of this Lagrangian density can be found in [26].
To study both the dynamics of the cosmological background and of the linear perturbations, we consider a threedimensional spatial metric given by a perturbed flat Friedmann-Lemaître-Robertson-Walker (FLRW) geometry, i.e. in the form where a(t) is the scale factor, ζ and E are scalar perturbations, β i is a vector perturbation and h ij a tensor perturbation. Notice that we have already decomposed the perturbations following the usual scalar-vector-tensor decomposition which implies that β i is divergenceless whereas h ij is transverse and traceless on the FLRW background, i.e., The remaining components of the metric, the lapse function and shift vector, are parametrized by where α and χ are scalar perturbations whereas χ i is a divergenceless vectorial perturbation, i.e. δ ij ∂ i χ j = 0. To finish with the variables entering in the gravitational Lagrangian density, one should also perturb the auxiliary fields λ 1 ,λ i and C according to, where, again, δλ 1 , δλ 2 , δC are scalar perturbations, and δλ 2i is a divergenceless vectorial perturbation, i.e. δ ij ∂ i δλ 2j = 0. Note that, as usual, the background quantities are time-dependent only, and the values ofλ i and N i vanish due to their vectorial nature and the background symmetry.
Concerning the variables entering specifically in the matter Lagrangian density (28), they are parametrized as follows where δJ, δj, and v are scalar perturbations. Now, we have introduced all the necessary variables to study the dynamics of the background and of the linear perturbations. For that purpose, we expand the total Lagrangian density (27) up to the second order in perturbations. The Lagrangian density linear in the perturbations provides the following equations of motion for the background, where f ,C and f ,CC are respectively the first and second derivatives of f with respect to C, H =ȧ/(aN ) is the Hubble constant and the so-called total particle number N tot is an integration constant obtained from integrating the conservation constraint ∇ µ J µ = 0 on a flat FLRW geometry. The last equation of (34) and the last two equations of (35) are the modified Friedmann equations and the usual conservation equation for the minimally coupled matter.
It is easy to verify that, in the limit f ,C → 1, we recover the GR equations of motion coupled to a perfect fluid.
Taking the time derivative of the first equation of (35), and using the modified Friedmann equations together with the conservation equation, we find the useful relation, Therefore, the equations of motion enable us to express the four quantities f,Ḣ, C andĊ in terms of f ,C , f ,CC and the other standard cosmological quantities a, H, ρ, . . . only. Hence, we can eliminate the variables f,Ḣ, C andĊ by a direct substitution, which we are going to do from now on. Expanding the Lagrangian density at second order in perturbations and Fourier transforming as usual with respect to the spatial comoving coordinates, we can then find, for high k ≡ | k| (with k the comoving momentum), the no-ghost conditions, together with the speeds of propagation. Tensor perturbations give so that letting f ,C > 0 at all times is sufficient to avoid the presence of ghost-like tensor modes. Furthermore, the famous constraint on the speed of propagation of gravitational waves today [44] (at redshift z = 0) imposes that c T ≈ 1 to the accuracy of order O(10 −15 ), which leads to the condition that (f ,C ) 2 ≈ 1 to the accuracy of order O(10 −15 ) at z = 0 as well. The analysis of scalar perturbations gives, Hence, we find that as long as f ,CC = 0, the speed of propagation changes for any matter component, dust included. Finally, and as expected, there are no propagating vector perturbations in the gravity sector.

B. Implementation of background equations of motion
Now, we are going to explain in detail the way we have implemented the background equations of motion in the code. First, we extend the theory to the case where several fluids are coupled to gravity (in order to model radiation, dust and eventually dark energy components) by adding to the total Lagrangian density (27) different matter Lagrangian densities. In that case, the two modified Friedmann equation and the conservation equation given in (34,35) become, where the subscript i parametrizes the different fluids, while the remaining equations of motion are unchanged, On using the equation for C(t) in (41), the first modified Friedmann equation in (40) can be rewritten as where we used the notations, , Hence, the modification of GR shows up, in the first Friedmann equation, as a fluid of effective density ρ f . Following a similar procedure we can also find the expression for the effective pressure P f or equivalently p f ≡ P f /(3M 2 P ). Indeed, using the definitions of f and i given in (43), the second modified Friedmann equation in (40) can be reformulated as follows, where we used the notations, Hence, p f is the (normalized) pressure of the effective fluid which describes the modifications of GR in a cosmological background. Finally, we also define the effective equation of state w f ≡ p f / f . Now, we are ready to show how to implement these equations of motion in the code. We start by choosing the lapse function so that it coincides with the scale factor, i.e. N = a, and then the time t becomes the conformal time τ . Hence, from now on, an overdot represents the derivative with respect to τ , i.e., for any function F . This choice is motivated by the fact that most of the Boltzmann solvers are written in terms of the conformal time. It is also convenient, from this section on, to introduce the derivative with respect to N e ≡ − ln(a/a 0 ), where a 0 is the scale factor today (at z = 0) for any function F . The τ and N e derivatives are related bẏ Then, we differentiate the first Friedmann equation in (40) with respect to τ and, after using the conservation equations (35) for each fluid, we obtain the dynamical equation, which is equivalent to As a consequence, the dark energy component does not contribute to the derivative of C and, if we denote by ρ m the dark matter density and by ρ r the radiation density, we have where H 0 and Ω i0 are the values of H and Ω i today (at z = 0). The second differential equation has been introduced so that we obtain an autonomous system of Ordinary Differential Equations (ODEs). We need to provide the initial conditions, namely C(N e = 0) = C 0 and a(0) = 1 to fully integrate the system 2 . Therefore, in order to have C at any value of the redshift z (z > 0) we just need to integrate the system of ODEs up to N e = ln(1 + z). Finally, on solving the system of ODEs (50) for a given function f (C) we are also able to find at any redshift the value of every relevant background quantity.
Let us now discuss the initial conditions. From (42), we know that today Ω f 0 = 1 − i Ω i0 . On combining this relation with the first Friedmann equation in (40), we find that Furthermore, since we impose from observations that [f ,C (z = 0)] 2 ≈ 1 to the accuracy of order O(10 −15 ) (see discussion below Eq. (37)), we find that C 0 ≈ −6H 2 0 , and we can safely 3 set this value as the initial condition for C.

C. Implementation of perturbation equations
In this subsection, we are going to compute the equations for the linear scalar perturbations and reduce them to a minimal but complete system where we have eliminated the redundant variables. We start with (7 + 3 s) scalar perturbations (where s is the number of matter components): for the metric, we introduced ζ, E, α, χ in (29)-(31), for the matter components we introduced δJ i , δj i and v i in (33), and finally we introduced δλ 1 , δλ 2 , δC in (32) for the remaining variables. As usual, we need to expand the total Lagrangian density for f (H) theories coupled to matter fields up to second order and we perform a Fourier transformation with respect to the spatial comoving coordinates to obtain the equations of the perturbations. In the sequel, we will denote by E Q = 0 the equation of motion of the perturbation variable Q obtained by deriving the quadratic Lagrangian density with respect to Q. For instance, we find and similar equations can be computed for the other variables. Of course, we make use of the equations of motion for the background to simplify several expressions, in particular, the background lapse function still coincides with the scale factor, namely N (τ ) = a(τ ), where τ is the conformal time. The first step consists in eliminating the redundant variables. We can see that we can integrate out the auxiliary variables δJ i , δj i in the matter Lagrangian density, only leaving the fields v i and δ i for each matter field components (see [26] for details). In the gravitational sector, we make use of the spatial gauge freedom (the invariance under three-dimensional diffeomorphisms) to fix E = 0 as it is often done. Furthermore, as we are going to see later on, it is convenient to use the same gauge invariant combinations of fields as the ones used in ΛCDM to describe perturbations in the Newtonian gauge. Hence, we make the following field redefinitions, Again, these combinations, in the language of GR, correspond to the gauge-invariant combinations which reduce to the Newtonian fields (when we adopt the Newtonian gauge).
2 As we need to find the background evolution at very high redshifts, we need a stable integrator. We found it in a Runge-Kutta Cash-Karp (4, 5) method. Instead for the Boltzmann code involving also the equations of motion for the perturbations we have used the default integrator of CLASS. 3 One could solve the initial condition equation C 0 = −6H 2 0 /[f ,C (z = 0)] 2 for C 0 iteratively. At lowest order, one finds C 0 = −6H 2 0 , and at next order C 0 = −6H 2 0 /[f ,C (C = −6H 2 0 )] 2 . However, for the models we have considered, the second solution is numerically indistinguishable from the value C 0 = −6H 2 0 . Now, we have all the ingredients to compute and simplify the equations for the perturbations. We start with the matter equations which are given by, where n i i,nn / i,n =ṗ i /˙ i and σ i is the shear perturbation for each fluid. We see that the form of the equations of motion in the matter sector is identical to the ones in GR in the Newtonian gauge. This is a consequence of the fact that the auxiliary field C in the f (H) theory does not have any direct coupling with matter fields, and matter is minimally coupled to the gravitational sector.
Then, we study the equations E δλ1 , E δλ2 and E δC . After a direct calculation, we see that we can integrate out the variables δλ 1 , δC and δλ 2 according to, Furthermore, on considering the equation of motion E χ = 0, we can also find χ as follows, Finally, it remains to study the equations which determine the dynamics for the perturbations in the gravitational sector, namely for φ and ψ. We start by substituting the Newtonian gauge fields (53) directly in the total Lagrangian density L so that we obtain a new equivalent Lagrangian density L for these new variables. Interestingly, we find that the equation for χ, that we denoteĒ χ = 0, is a new one and can be written as where we have introduced the notations, In the limit where f ,C → 1 and f ,CC → 0, we recover the expected results of GR and, (59) represents one of the two fundamental equations which involve the metric perturbation variables, that Boltzmann solvers need to solve besides the equations of motion for the matter variables. At this stage, we need another equation for ψ. We find it by the following procedure. First of all we consider E ψ = 0. This equation can be written as Then, we computeĖ ψ and we obtain an equation that involves the termsφ,δ i andθ i . Therefore from linear combinations of this equation with E 1 and the matter equations of motion, as well asĒ ψ , we can remove each of these derivatives. Finally we obtain the remaining expected equation, which, in the GR limit, reduces to the correct shear perturbation equation.
We have now all the dynamical equations of motion necessary to implement the dynamics of the whole system: the matter equations (54) and (55) together with the "gravitational" equations (59) and (62) couple the perturbation variables θ i , δ i , φ and ψ. Notice that this theory does not add any new propagating degree of freedom and this is the reason why we do not need to add any new equation to the Boltzmann code. The main difference between this theory and GR lies, at the level of perturbations, on the two modified equations of motion for the metric perturbation variables.

IV. A CONCRETE EXAMPLE: THE "KINK" MODEL
In the following, we will consider one explicit example of f (H) theories by choosing a particular class of functions for f . Before presenting the model in details, we give our motivations and also some of the most important results we have obtained with it.
This model has been built so that it is consistent with the constraint c T = 1 at low redshifts, which means f ,C (z = 0) ≈ 1, and it possibly addresses some so-far unsolved cosmological puzzles. We have found that the model we are to introduce is capable of fitting the chosen data sets better than ΛCDM. Even if this model does not improve much the fit to the H 0 measurement, it achieves a better fit especially for Planck data. Furthermore, although in the context of spatially curved (Ω K0 = 0) ΛCDM, one is able to fit Planck data alone better, as soon as we introduce late time data, BAO data in particular, Ω K0 = 0 is strongly constrained. Some authors [17] are led to state that ΛCDM is ruled out in the light of such behavior. However, the model we are going to introduce, does not feel such a strong constraint from late time data. Before and after the insertion of late time data, Planck data are considerably fit better by the model we are about to discuss. Hence, this study opens a window on using modified gravity models to address not only late time cosmology but even cosmology at intermediate/high redshifts. Before defining our model, let us remind that ΛCDM corresponds to the case where f (C) is an affine function of C, i.e. of the form f (C) = C − C 0 + f 0 with f 0 ≡ f (C 0 ). The model we are now considering consists in choosing f (C) such that its derivative is of the form, where a 1 , a 2 and a 3 are free positive real parameters. The fact that a 1 > 0 ensures f ,C > 0 which is the condition to avoid the presence of ghost-like tensor modes (37). The hyperbolic tangent function gives f ,C a kink shape and, for this reason, we dub it the "kink"-model. A direct integration with respect to the C variable leads to where f 0 is the integration constant, C 0 is the value of C today, and we have made use of the softplus function 4 . Some of the free parameters are constrained by the initial conditions (51). Indeed, as C 0 ≈ −6H 2 0 and f ,C ≈ 1 at z = 0, then the parameters a 2 and a 3 must satisfy the condition (a 2 − 6)/a 3 1 which implies that a 2 a 3 . Furthermore, on using the Friedmann equation, we already know that f 0 = −6H 2 0 (1 − Ω f 0 ). Now, let us study some properties of the model at very early times when C < 0 and |C| a 2 H 2 0 . In that case, the function f (C) simplifies and becomes, As a consequence, at early times, the theory is equivalent to a GR theory with f (C) an affine function of C, i.e. of the form f (C) = αC + β, where however the constants α and β are different with respect to ΛCDM, i.e. α = 1, in general. Hence, on using the expression of C given in (35) and the definitions (43), we find at early times that which means that f has non-trivial contributions at early times. In particular, this implies that Ω r does not go to unity at early times, in general. Then one may wonder if this is enough to rule out at once the model. To see this is not the case, let us study the behavior of the effective equation-of-state for the f (C) component at early times, for instance during radiation domination when Ω r Ω m . A direct calculation shows that where we have used Eq. (66). Therefore this model gives an effective contribution to radiation at early times.
Therefore when w f ≈ 1/3, then f will contribute to the effective radiation energy density eff r = r + f so that, on the background, the effective Ω eff r will have an extra component due to Ω f and lim z→∞ Ω eff r = 1, as we expect. We want to make clear that this model is not a model of "dark radiation" for at least two reasons: 1) it gives an effective radiation component only during radiation domination, whereas at late times, ρ f , on the background, behaves instead as an effective cosmological constant; 2) no extra (massive or massless) degree of freedom is introduced in the theory at any time, even during radiation domination.

V. RESULTS
In this section, we present the results obtained from running a Monte Carlo sampling of the parameter space for the kink model. Sampling the parameter space helps understanding the features of the model, as otherwise we would find it hard to predict which part of the parameter space is more appealing with respect to cosmology. The fitness parameter χ 2 is an important quantity, which tells how much the data prefer the model under investigation in the parameter estimation. Whenever a model has a better, i.e. lower, χ 2 than the standard cosmological model ΛCDM, that model deserves attention in the context of cosmological tensions. As we already wrote in the previous section, we will see that the kink model (64) consistently gives a better fitness parameter than ΛCDM. For this reason, we report on the kink model.
We present here our study of the behavior of both ΛCDM and the kink model for several data sets constraining the evolution of the background and perturbations both at high and low redshifts. We used a Monte-Carlo sampling in order to find the bestfit parameter points for both models which minimize the χ 2 . We gave very broad priors on the free parameters of the kink model. In particular we set −1 < a 1 < 50 (with a flat prior), 0.5 < log 10 a 2 < 18 (with log-flat prior), and −9 < log 10 β < −0.75 (with log-flat prior) where a 3 = β a 2 . We set such large priors as the model was to be applied to a large redshift range, without a priori knowing whether the procedure would find any good fit at all in the allocated sampling time. However, in the limit a 1 → 0 the kink model reduces to ΛCDM. We knew therefore that there should have been at least some non-zero good-fit interval for the parameters. Despite those very broad priors, the procedure successfully converged to a good fit.

A. Comparison of χ 2
We found that the kink model performs better than ΛCDM: the fitness parameters χ 2 (total and respective to each experiment) are compared to those obtained in the ΛCDM case in Table I. The difference in χ 2 between the two models is ∆χ 2 = 16.6. Even though we have three more parameters, we find that the kink model is preferable and fits the data much better than ΛCDM. To be more precise, according to the Akaike Information Criterion (AIC), for each new parameter introduced by a model, effectively the χ 2 should be roughly speaking increased by a value ≈ 2 [45,46], in order to allow a more accurate comparison. Even doing so, as we have 3 new parameters, we still find a significative difference in χ 2 given by ∆χ 2 AIC ≈ 16.6 − 2 × 3 = 10.6, which shows that the kink model has a much larger probability to be a better model to fit the data, compared to ΛCDM. This result would become even stronger in total χ 2 = 3473.27 in total χ 2 = 3456.67 Table I. Comparison of the total χ 2 for the bestfit for both ΛCDM and the kink model for all the considered data sets.
if we would be fixing the parameter log 10 β to some value within the 2σ constraint. In fact, as we will show later on, we have a large degeneracy on the parameter β, which can therefore be fixed a priori. Of course, as an alternative, it is still possible that some yet unknown systematics are tilting the balance in disfavor of ΛCDM, especially for early time data. Before studying the bestfit of the model, we want to make a few statements about this result. First of all we can see that the kink model, performs better than ΛCDM on the HST data point, thus alleviating the tension in measurements of the value of today's Hubble factor to some degree. The model performs in a way very similar to ΛCDM for the other late time data, consisting mainly of BAO and Supernova Type Ia data (JLA). The fact that the model is similar to ΛCDM on these last experiments is afterall maybe not a surprise because the kink model was constructed as to reduce to ΛCDM at late times (with possibly different values for the background parameters, like H 0 for instance).
However, there is an important difference at early time when f ,C = 1. In fact, we find that the main contribution to a lower value of the χ 2 takes place in the redshift interval, which goes from today up to values of redshift sensitive to Planck results, and the kink model performs consistently better than ΛCDM for any of the Planck experiments: ∆χ 2 TTTEEE = 12.53, ∆χ 2 lowl EE = 1.01, and ∆χ 2 lowl TT = 1.55. The reason for such an improvement is due to both a modification of the background (an effective radiation component, which at late times changes into a cosmological constant contribution, and, as we will see, to a fast change of the value of H at intermediate redshift z 743), and to a different dynamics for the perturbations (whenever f ,C = 1). With these considerations, it is now time to focus our attention to the constraints and the bestfit of the model in order to understand its characteristics.

B. Two dimensional likelihoods and bestfit
On analyzing the chains obtained after Monte-Carlo sampling we show the two-dimensional marginalized likelihoods for the cosmological variables of interest in Fig. 1. The results for the one-dimensional 2σ constraints are summarized in Table II First of all, it is interesting to see the bound on a 1 , which evidently states that a 1 is different from 0 at 2σ, i.e. the data suggests non-negligible deviation from ΛCDM. The bound on a 2 is also interesting as we have a lower and an upper limit. In particular, the value of a 2 selects the energy at which f ,CC changes significantly. We will comment later on at which redshift this change occurs. Finally let us briefly discuss about the parameter β, which has been defined above as a 3 = β a 2 . The value of β does not affect the value of the χ 2 when it is sufficiently small. In fact, in order to study better the dependence of the χ 2 on the parameter β, we have performed the experiment of changing the value of β for the best fit and confirmed that its value does not affect the minimum of χ 2 , as shown in Fig. 2 which plots the behavior of the likelihood when β varies, as long as log 10 β −4.1, for higher precision. This means that we have a degeneracy for small values of β. In fact, the allowed lower bound was hitting the lower prior limit (which was set to be 10 −9 ). Much smaller values would in fact require a sufficiently high precision which is not allowed in our computations. The meaning of a small β is anyhow simple, the redshift range for which f ,CC changes significantly cannot be too large, i.e. the transition between the two different f ,C cannot be adiabatic. Later on, when we study in detail the minimum for the whole data sets, we will also try to understand the reason why large values of β seem to be excluded.
As already stated above for these all-redshift-range data sets, we can see that the bestfit kink model alleviates the H 0 tension slightly, but, mostly, it gives a better fit on Planck data compared with ΛCDM. In fact, a value for ∆χ 2 = 16.6 is large enough (even on considering three more parameters for the kink model) as to exclude ΛCDM at 2σ (because a 1 is different from 0 at 2σ, and the a 1 → 0 limit is a smooth one, as we have no other propagating mode than ΛCDM). This result is encouraging for the kink model and at least tells that flat-ΛCDM does not necessarily give the better fit to Planck data, or we may conclude that, if flat-ΛCDM is the only allowed model a priori, Planck data indicate some unknown systematics. A simple solution, in this sense following Occam's razor, is that models different from ΛCDM have an arena to test gravity already in the Planck data sets and at high redshifts.
In summary, we have found a model which fits Planck data better than flat ΛCDM. The fact that there are models which fit Planck data better than flat ΛCDM is not new (see e.g. [16,17,47], the last one describing a model in the context of Horndeski theories). For example, it is now well known that Planck data prefer Ω K = 0 [16,18]. However, a non-flat model is then in strong tension with BAO and inflationary paradigm at the same time [17]. On the contrary, we have shown here that BAO data do not constrain much this model (or they constrain it as much as flat-ΛCDM) and, as a consequence, do not spoil the fact that the kink model still performs better than flat ΛCDM.

C. Background evolution
In the following we study in detail the behavior of the minimum of the χ 2 starting with the background evolution. We have fixed the various background parameters to their bestfit values and integrated the background equations of motion. The results of the evolution of the three Ω i functions are given in Fig. 3. 10 100 1000 10000 100000 We can see that some non-trivial dynamics happens at intermediate redshifts (z 743). In order to understand better what happens we further study the dynamics of the background for several variables. One variable of interest is ε H ≡ −Ḣ/(N H 2 ) (which can be calculated via Eq. (44)), as it shows whether the universe accelerates (when ε H < 1), decelerates (approaching 2 during radiation domination), and if some non-trivial singularity is present (besides the big-bang). In Fig. 4, we can see that around z 743, there is a fast but smooth transition at which ε H grows but remains finite. Therefore no singularity is present. Furthermore, we can see that this behavior takes place in an interval ∆z < 1. Finally this variable shows that the universe starts being radiation dominated, then matter dominated, and finally the universe starts accelerating.  On the left panel we can see that today the universe accelerates (since εH < 1), whereas at early times the universe is radiation dominated (as εH → 2). On the right panel, we zoom around the redshift z 743, and we show that the transition is fast but smooth, i.e. no singularity is present (and this also implies that numerics are stable).
Another variable of interest is represented by w f = p f /ρ f , i.e. the effective equation of state for the f -component. Its behavior is shown in Fig. 5, where we can see that the f -component, as already mentioned before, behaves as radiation at early times, but as a cosmological constant at late times.   Figure 6. Evolution of the squared speeds of propagation for the scalar-matter modes and tensor modes. On the left panel we can see that the squared speeds of propagation for both dust (c 2 m ) and radiation (c 2 r ) change, and their values become negative during the transition. For a short time gradient instabilities take place, but as we will see later on, it is not so catastrophic (because the transition is short enough) as to make the matter perturbations totally unstable. This implies that β (the parameter which determines how fast the transition is) must be small enough so that the instability-era is negligible. On the right panel, we show instead that at early times the speed of tensor modes was not unity, but today it is indistinguishable from unity as to pass the multi-messenger constraints.  Figure 7. Dynamics of the perturbation variable δγ, the photon energy density contrast. It is clear from both panels (the right one describing its evolution during the transition-era), that such a variable still evolves in a way similar to GR after this transition ends). For this plot we have used the same bestfit parameters of the background, and on top of that ω b = 0.02284032, ωc = 0.1184566, τreio = 0.05183895, ln(10 10 As) = 3.039270, ns = 0.9778259, and, as for the size of the wave-vector, k = 0.1 Mpc −1 (to enhance the effect of the instability in the very short wavelengths, but still linear).

D. Evolution of the perturbations
After having studied in details the behavior of the background, we now focus on the dynamics of the perturbations. First of all, let us see what happens for the speed of propagation for the modes. In Fig. 6, we can see that the perturbation variables undergo an era of instability (because c 2 s < 0), whereas tensor modes acquire a velocity of propagation different from unity (at early times, i.e. for z > 743). The fact that there is an era of instability should make us worried, however such an era is relatively short 5 . Therefore, there is the chance that the matter perturbations 5 We find that in cosmic time, for the bestfit, this era lasts ∆t    Figure 9. Dynamics of the perturbation variables θ b and θc (related to the velocity perturbations for the baryons and the CDM components respectively). Also for these variables, after the transition-era ends, their evolution is not too different from GR.
still remain finite after that transition-era ends. This is exactly what happens. After all, we already know that the χ 2 is really good for this model. In fact, the χ 2 from Planck comes from correlations of quantities (for example TT correlations) which typically involve an integral over the line-of-sight. Therefore, we expect that, if perturbations follow again GR evolution after the transition-era, the results for any such integral will be smooth. Since linear perturbations remains under control, one may expect a similar behavior for higher order perturbations, although we have not studied here this complicated issue. Let us then study the evolution of the perturbation variables. As shown in Fig. 7, the photon energy density contrast δ γ remains finite and once the transition is over, its evolution returns to the one similar to GR. Similar considerations can be made for the other perturbation variables, as shown in Figs. 8, 9 and 10. To conclude the study of the best fit model, we want to emphasize again that the transition era is short-lived enough as to make the instability inefficient so that as soon as it ends, the dynamics is once again following the GR evolution. Therefore this model changes the χ 2 not because of the transition, but because of the different behavior which it shows at early times (i.e. for z > 743). It should be noticed that the numerics have been always stable as to be sensitive to small variations of a 1 (which on the bestfit is about +2 × 10 −3 , whereas a 1 = 0 corresponds to ΛCDM) and of the parameter β governing the duration of the transition era. This sensitivity has been achieved by reducing the standard tolerance parameters ( 10 −5 ) which determine the precision of the numerical solution of the  Boltzmann solver, CLASS, by about ten orders of magnitude.

VI. CONCLUSION
We have studied a "kink model" within the class of minimally modified gravity theories dubbed f (H) theories, which was shown in the present work to have a fit to several early times and late times data sets better than ΛCDM by ∆χ 2 = 16.6. We are impressed by several aspects of it. Most importantly, this model shows once for all that Planck data and ΛCDM do not fit so well with each other. In particular, the kink model performs better than ΛCDM especially on Planck 2018, but also on H 0 single-point data sets.
In the context of ΛCDM this aspect is ameliorated for Planck data alone notably if one introduces a non-zero curvature term, which however is not well motivated by inflation and moreover it is strongly suppressed by late time data (especially BAO). On the other hand, late time data do not strongly constrain the kink model, at most as much as flat-ΛCDM, and the kink model still keeps a better fit to Planck 2018 data.
On top of this consideration, the bestfit kink model is characterized by the presence of a transition era which happens at a redshift of z 743 6 . Such a transition era happens at intermediate redshifts and not at very high energies. This means that there is the possibility that such a transition may occur in other contexts with strong gravity (inside a star for example). Furthermore during this transition era, the scalar perturbations have an instability (c 2 s < 0) which, at least, as long as linear perturbation theory holds, does not seem to be catastrophic because of its short duration (∆z 0.1 around z 743) so that the matter perturbations (for which we have many constraints) do not show any divergent behavior. In particular, this transition era seems to be washed out when we look at observables which imply a line-of-sight integration.
The kink model then opens up an arena for modified gravity models, which so far are introduced mostly to address late time data only. This model in fact, addresses and resolves tension for ΛCDM in Planck data and at the same time shows a non-trivial window for phenomenology at intermediate redshifts. Future work is needed to address several points which remain to be explored, such as the existence, composition and evolution of compact objects for this model.