Fanni M Sélley
email: selley@lpsm.paris

Matteo Tanzi
email: matteo.tanzi@sns.it

Linear Response for a Family of Self-Consistent Transfer Operators

We study a system of globally coupled uniformly expanding circle maps in the thermodynamic limit. The state of the system is described by a probability density and its evolution is given by the action of a nonlinear operator called the self-consistent transfer operator. Self-consistency is understood in the following sense: if ϕ corresponds to the system's state, then the evolution of ϕ is given by the application of the transfer operator of the circle map F ϕ . This is a ϕ-dependent map describing the dynamics of a single unit in the finite system where the interaction coming from all the nodes is replaced with the average of the interaction term with respect to the density ϕ. Assuming some level of smoothness of the coupled circle maps and of the coupling, we prove that when the coupling strength is sufficiently small, the system has a unique smooth stable state. We then show that this stable state is continuously differentiable as a function of the coupling strength. Finally, we prove that the derivative satisfies a linear response formula which at zero simplifies to a formula reminiscent to the one obtained for the linear transfer operators of perturbed expanding circle maps.

Introduction

A fundamental question in the theory of dynamical systems is how the statistical properties of a system change when it is subjected to perturbation. A system exhibits linear response if the invariant measure depends smoothly on the perturbation, and an expression for the derivative of the invariant measure (in the strong or in the dual sense) is called a linear response formula.

The rigorous study of linear response dates back to Ruelle, who proved linear response of uniformly hyperbolic Axiom A systems [START_REF] Ruelle | Differentiation of SRB states[END_REF][START_REF] Ruelle | General linear response formula in statistical mechanics, and the fluctuation-dissipation theorem far from equilibrium[END_REF][START_REF] Ruelle | A review of linear response theory for general differentiable dynamical systems[END_REF]. Similar results can be proved in some non-uniformly hyperbolic [START_REF] Dolgopyat | On differentiability of SRB states for partially hyperbolic systems[END_REF][START_REF] Zhang | On the smooth dependence of SRB measures for partially hyperbolic systems[END_REF] or non-uniformly expanding cases such as intermittent maps [START_REF] Bahsoun | Linear response in the intermittent family: differentiation in a weighted C 0 -norm[END_REF][START_REF] Baladi | Linear response for intermittent maps[END_REF][START_REF] Korepanov | Linear response for intermittent maps with summable and nonsummable decay of correlations[END_REF] and piecewise expanding unimodal maps [START_REF] Baladi | Linear response formula for piecewise expanding unimodal maps[END_REF][START_REF] Baladi | Alternative proofs of linear response for piecewise expanding unimodal maps[END_REF]. But caution is required, as examples for the lack of linear response are also well known [START_REF] Aspenberg | On the fractional susceptibility function of piecewise expanding maps[END_REF][START_REF] Baladi | On the susceptibility function of piecewise expanding interval maps[END_REF][START_REF] Baladi | Whitney-Hölder continuity of the SRB measure for transversal families of smooth unimodal maps[END_REF][START_REF] Baladi | Linear and fractional response for the SRB measure of smooth hyperbolic attractors and discontinuous observables[END_REF] (see [START_REF] Baladi | Linear response, or else[END_REF] for a more comprehensive collection of references).

Deterministic systems subjected to random perturbations are also known to exhibit linear response in some cases [START_REF] Galatolo | A linear response for dynamical systems with additive noise[END_REF][START_REF] Galatolo | Quadratic response of random and deterministic dynamical systems[END_REF]. In addition to theoretical works, linear response theory is successfully applied in geophysics, in particular climate science [START_REF] Lucarini | Mathematical and physical ideas for climate science[END_REF][START_REF] Majda | Challenges in climate science and contemporary applied mathematics[END_REF].

In this paper we rigorously prove a linear response formula for all-to-all weakly coupled uniformly expanding maps in the thermodynamic limit. To the best of our knowledge, there are few results available on linear response of complex systems with multiple interacting units. A program focusing on the control of statistical properties of extended systems by exploiting linear response was outlined by Mackay [START_REF] Mackay | Management of complex dynamical systems[END_REF]. This question was further pursued in the case of a single dynamical unit defined by a uniformly expanding circle map [START_REF] Galatolo | Controlling the statistical properties of expanding maps[END_REF] and by a smooth map of a compact smooth Riemannian manifold [START_REF] Kloeckner | The linear request problem[END_REF]. In case of coupled map lattices, the smooth dependence of the SRB measure and the linear response formula for was studied by Jiang and de la Llave [START_REF] Jiang | Linear response function for coupled hyperbolic attractors[END_REF][START_REF] Jiang | Smooth dependence of thermodynamic limits of SRB-measures[END_REF]. Wormell and Gottwald [START_REF] Wormell | On the validity of linear response theory in high-dimensional deterministic dynamical systems[END_REF][START_REF] Wormell | Linear response for macroscopic observables in high-dimensional systems[END_REF] recently showcased numerical examples and some rigorous arguments for all-to-all coupled maps with mean-field interaction showing that it is possible for a high dimensional compound system to exhibit linear response, even if its units do not.

Our setup is based on the model studied by Fernandez [START_REF] Fernandez | Breaking of ergodicity in expanding systems of globally coupled piecewise affine circle maps[END_REF]. The finite size model is described by M coordinates (x 1 , ..., x M ) with x i ∈ T, whose evolution is given by

x i (t + 1) = f • Φ ε (x i (t); x 1 (t), ..., x M (t)) (1)
where Φ ε : T × T M → T,

Φ ε (x i ; x 1 , ..., x M ) = x i + ε M M j=1 h (x i , x j ) , (2)
is a smooth mean-field all-to-all coupling between the coordinates (ε is the coupling strength * ), and f : T → T is a sufficiently smooth map describing the uncoupled dynamics.

When taking the thermodynamic limit, M → ∞, one can replace the finite dimensional state space T M with the infinite-dimensional state space T N , and the coupling term in [START_REF] Bahsoun | Linear response in the intermittent family: differentiation in a weighted C 0 -norm[END_REF] with lim

M →∞ ε M M j=1 h (x i , x j ) .
The limit above might not exist. Let's assume that for the initial condition (x i (0)) ∞ i=1 ∈ T N there is a probability measure µ such that lim

M →∞ 1 M M i=1 δ x i (0) = µ (3)
where the limit is with respect to the weak convergence of measures † . This assump-tion ensures that for any continuous h,

lim M →∞ ε M M j=1
h(x i (0), x j (0)) = ε h(x i (0), y)dµ(y), and the time evolution of each coordinate can be written as:

x i (1) = f • Φ ε,µ (x i (0)) where Φ ε,µ (x) = x + ε h(x, y)dµ(y).

Calling f * and (Φ ε,µ ) * the push-forwards of f and Φ ε,µ , lim

M →∞ 1 M M i=1 δ x i (1) = f * (Φ ε,µ ) * (µ)
where the convergence is weak and follows from the continuity of the push forwards (which only requires the continuity of f and Φ ε,µ ). In case µ has a density ϕ with respect to a suitable reference measure (Lebesgue, for instance) it is sufficient to study the evolution of the density. By the above argument, this is given by a transfer operator that depends on the probability measure, thus ϕ itself (for a precise definition, see [START_REF] Baladi | Linear response formula for piecewise expanding unimodal maps[END_REF]). This nonlinear application is called a self-consistent transfer operator. This will be the object under study in this paper.

Fernandez [START_REF] Fernandez | Breaking of ergodicity in expanding systems of globally coupled piecewise affine circle maps[END_REF] studied the finite setting in the case of h(x, y) = g(x -y) with g : R → R defined to be a specific discontinuous, but piecewise affine function. Studying the thermodynamic limit of his setup in case of ε sufficiently small, Bálint and Sélley [START_REF] Sélley | Mean-field coupling of identical expanding circle maps[END_REF] first showed that when f (x) = 2x mod 1, the constant function is the only smooth fixed density. Expanding these results, in [START_REF] Bálint | Synchronization versus stability of the invariant distribution for a class of globally coupled maps[END_REF], they considered a general uniformly expanding smooth circle map f , and showed that the invariant density is unique and a Lipschitz continuous function of the coupling strength, with higher regularity being obstructed by the discontinuity of the function h.

In this paper we present an approach based on convex cones of functions to study self-consistent operators for general uniformly expanding f and general h satisfying certain smoothness conditions. We show that in the weak coupling regime (ε small), the self-consistent operator has a unique fixed smooth density that depends in a differentiable way on the coupling strength. Furthermore, we compute a linear response formula for the density.

We believe that assuming stronger smoothness conditions on h and f , higher smoothness of the invariant density in terms of the coupling could be also achieved and a quadratic response formula could be deduced, reminiscent of the formula obtained in [START_REF] Galatolo | Quadratic response of random and deterministic dynamical systems[END_REF]. Another, more difficult, open question is whether one could apply ideas similar to the one in this paper to treat the case where h is not smooth. This would require the definition of some suitable Banach spaces and cones playing analogous roles to the spaces and cones of smooth functions that we consider here.

Setup

Function spaces and norms Let S 1 denote the unit circle and denote the Lebesgue (Haar) measure on S 1 by λ. Throughout this paper, we consider the representation of S 1 as the unit interval [0, 1] with 0 and 1 identified. We denote by L 1 (S 1 , R) the real valued functions on S 1 for which

S 1 |f | dλ < ∞ holds. Let C 0 (S 1 , S 1
) denote the space of continuous mappings of S 1 , endowed with the norm

f C 0 = max x∈S 1 |f (x)|.
For integers k ≥ 1, let C k (S 1 , S 1 ) denote the space of k-times continuously differentiable mappings of S 1 , endowed with the norm

f C k = max x∈S 1 |f (x)| + k ℓ=1 max x∈S 1 |f (k) (x)|.
Uncoupled Dynamics Let f ∈ C k (S 1 , S 1 ) (with k > 1 to be specified later) be a uniformly expanding circle map, i.e. there exists an ω > 1 such that |f ′ (x)| > ω for all x ∈ S 1 , and denote by P its transfer operator with respect to the Lebesgue measure. That is, P :

L 1 (S 1 , R) → L 1 (S 1 , R) is defined as P ϕ(x) = y∈f -1 (x) ϕ(y) |f ′ (y)| , ϕ ∈ L 1 (S 1 , R). (4)
Mean-Field Coupling Given ε 0 > 0, t ∈ [-ε 0 , ε 0 ], h ∈ C k (S 1 × S 1 , R) (with k > 1 to be specified later), and ψ ∈ L 1 (S 1 , R + 0 ) ‡ with S 1 ψ = 1, define Φ t,ψ : S 1 → S 1 Φ t,ψ (x) = x + t S 1
h(x, y)ψ(y)dy.

We assume that ε 0 > 0 is sufficiently small so that Φ t,ψ is a diffeomorphism such that Φ ′ t,ψ > 0, and we denote by Q t,ψ its transfer operator. More precisely, similarly to (4), the operator

Q t,ψ : L 1 (S 1 , R) → L 1 (S 1 , R) is defined as Q t,ψ ϕ(x) = ϕ Φ ′ t,ψ • Φ -1 t,ψ (x), ϕ ∈ L 1 (S 1 , R). (5)
The evolution of the density of states ψ under the mean-field coupling is then given by the nonlinear application of Q t defined as Q t (ψ) := Q t,ψ ψ. In the following we use the short-hand notation A ψ (x) = S 1 h(x, y)ψ(y)dy, so we write Φ t,ψ (x) = x + tA ψ (x). Thus substituting ϕ = ψ in (5) gives the following expression for Q t :

Q t ψ(x) = ψ 1 + tA ′ ψ • Φ -1 t,ψ (x). ( 6
)
Coupled Dynamics Define F t,ψ := f • Φ t,ψ , and let L t,ψ := P Q t,ψ denote its transfer operator. Given t ∈ [-ε 0 , ε 0 ], define the self-consistent transfer operator

L t : L 1 (S 1 , R) → L 1 (S 1 , R) as L t := P Q t (7)
which describes the evolution of the all-to-all coupled system in the thermodynamic limit. ‡ R + 0 is the set containing the positive real numbers together with 0.

Smooth densities Throughout the paper we are going to work with the restrictions of P, Q t and L t to some appropriate spaces of smoothly differentiable functions. These restrictions are well-defined: on the one hand, ϕ

∈ C ℓ implies P ϕ ∈ C ℓ for ℓ < k supposing f ∈ C k (S 1 , S 1 ). On the other, if h ∈ C k ′ (S 1 × S 1 , R) then Φ t,ψ ∈ C k ′ (S 1 , S 1 ) for sufficiently small values of t, thus Q t ϕ ∈ C ℓ for ℓ < k ′ .
By continuity, f is an N -fold covering map of S 1 for some N ∈ N. We make a technical assumption on f ensuring that the expansion is large enough and the distortion is small enough, more precisely

N max |f ′′ | ω 3 + 1 ω 2 < 1. ( 8
)
Note that for example the map f (x) = kx mod 1 (where k is an integer) and moderate perturbations satisfy this condition. We remark that we could exclude this condition by working with higher iterates of the self-consistent transfer operator. However, this would complicate the calculations a great deal and we will assume [START_REF] Baladi | Alternative proofs of linear response for piecewise expanding unimodal maps[END_REF] for the sake of clarity of the presentation. From now on, K > 0 will denote an unspecified positive number that depends on h and bounds of its derivatives only, while K > 0 is another unspecified constant depending also on f , unless otherwise specified.

Summary of the results and outline of the paper

The first theorem of this paper claims that for t sufficiently small, the system has a unique equilibrium density ρ(t), and any sufficiently smooth initial state will evolve towards ρ(t).

Theorem 1.1. Assume f ∈ C 2 (S 1 , S 1 ) and h ∈ C 2 (S 1 × S 1 , R). Then, there is ε 0 > 0 sufficiently small such that for any t ∈ [-ε 0 , ε 0 ], L t has a unique fixed nonnegative probability density ρ(t) ∈ Lip(S 1 , R + ). Furthermore, there exist constants K > 0 and Λ < 1 such that for any non-negative ϕ ∈ Lip(S 1 , R)with ϕ = 1 it holds that L n t ϕ -ρ(t) C 0 ≤ KΛ n for all n ∈ N.

Remark 1.1. A very similar result was stated by [START_REF] Keller | An ergodic theoretic approach to mean field coupled maps[END_REF]Theorem 4]. However the coupling mechanism there is different in the sense that an integral average is computed from the current system statistic which feeds into the dynamics as a constant. Our coupling is more complicated: it varies from point to point due to the x-dependence of the function h. However, the results and the methods that can be used to prove them are essentially the same. We view this statement as a preparatory result for the next theorem, our main result.

The second theorem claims that taking f and h sufficiently smooth, the dependence of ρ(t) on t is C 1 in a strong sense.

Theorem 1.2. Assume f ∈ C 5 (S 1 , S 1 ) and h ∈ C 5 (S 1 × S 1 , R). Then (i) there is an ε 1 > 0, ε 1 ≤ ε 0 such that ρ(t) ∈ C 3 (S 1 , R) for all t ∈ [-ε 1 , ε 1 ]; (ii) there is an ε 2 > 0, ε 2 ≤ ε 1 such that t → ρ(t) belongs to C 1 ([-ε 2 , ε 2 ], C 1 (S 1 , R));
(iii) there is an ε 3 > 0, ε 3 ≤ ε 2 such that the following linear response formula holds for | t| < ε 3 :

∂ t ρ(t)| t= t = -(1 -P t + tP K t) -1 P K t(ρ( t)),
where P t = P Q t,ρ( t) and the operator K t is given by (10).

Remark 1.2. Both theorems depend intimately on the fact that when the coupling strength is sufficiently weak, the operators satisfy uniform Lasota-Yorke type inequalities implying uniform spectral properties. However, standard perturbation theory arguments are insufficient to obtain the results above due to the nonlinear character of the operators L t . Getting around the nonlinearity of the operators is where the main novelty of our approach lies.

We break down the proof of the above theorems into several propositions listed below and proved in the following sections.

In Section 2 we study existence, uniqueness and regularity of the fixed point for the self-consistent transfer operator L t [START_REF] Baladi | Linear response formula for piecewise expanding unimodal maps[END_REF]. To this end, we consider the cone of log-Lipschitz functions

V a := ϕ : S 1 → R + : ϕ(x) ϕ(y) ≤ exp(a|x -y|) ,
where |x -y| denotes the Euclidean distance between x and y on the unit circle, and show that for large enough a, L t keeps V a invariant. Most importantly, we can also show that L t is a contraction with respect to the Hilbert metric θ a (see [START_REF] Boyarsky | Laws of chaos: invariant measures and dynamical systems in one dimension[END_REF] below) on V a .

Remark 1.3. Cones of functions and the Hilbert metric are common tools to prove uniqueness of the invariant density for some linear transfer operators [START_REF] Liverani | Decay of correlations[END_REF]. The standard approach uses the fact that a linear operator mapping a cone with finite diameter strictly inside itself is a contraction in the Hilbert metric. Unlike the standard case, L t is a nonlinear operator, and one has to prove explicitly that L t is a contraction.

Proposition 1.1. There is a 0 > 0 sufficiently large such that for a > a 0 , there is ε 0 > 0 sufficiently small such that for t ∈ [-ε 0 , ε 0 ], L t keeps V a invariant. Furthermore, there are a > 0, and

ε 0 > 0 small such that L t is a contraction on (V a , θ a ) for any t ∈ [-ε 0 , ε 0 ].
The proof of Theorem 1.1 then follows by standard arguments showing that under repeated application of L t , all the elements in the cone V a converge exponentially fast to the fixed point ρ(t) (in the C 0 norm).

To show that the C 5 -smoothness of f and h implies that ρ(t) ∈ C 3 (S 1 , R), define the following set for the quadruple

C = (C 1 , C 2 , C 3 , C 4 ) ∈ (R + ) 4 : C C := ϕ ∈ C 4 (S 1 , R + ) : ϕ = 1, d i ϕ dx i C 0 < C i i = 1, ..., 4 .
We show that Proposition 1.2. There is

ε 1 > 0 and C 1 , C 2 , C 3 , C 4 > 0 such that for t ∈ [-ε 1 , ε 1 ] L t C C ⊂ C C .
The existence of a fixed point in C 3 is then implied by the fact that L t is continuous (Lemma 2.3 below), C C is relatively compact in C 3 , and the Schauder fixed-point theorem. Since one can find ϕ ∈ C C ∩ V a (for example take ϕ = 1), the fixed point ρ ∈ V a , whose existence was argued above, must belong to the closure of C C and therefore is in C 3 , that will prove point (i) of Theorem 1.2.

In Section 3 we show that, when f and h are C 5 , the curve of fixed points t → ρ(t) found above is

C 1 from [-ε 2 , ε 2 ] to C 1 (S 1 , R) for some ε 2 > 0.
Remark 1.4. In the case of a one-parameter family {f t } t of sufficiently smooth uniformly expanding maps, the strong differentiability in t follows from the spectral properties of the associated transfer operators and the smoothness of t → f t [START_REF] Baladi | Linear response, or else[END_REF].

Here, since the one-parameter family depends on the fixed density itself, i.e. ρ(t) is the fixed density of f t = F t,ρ(t) (self-consistent relation), the usual arguments do not apply since there is no a priori knowledge of the regularity of t → f t that has to be established first. This is achieved by showing that the set of curves

Ĉε 2 ,K 1 ,K 2 := γ : [-ε 2 , ε 2 ] → C C s.t. (t,x) →γ(t)(x) is C 2 sup t d dt γ(t) C 3 <K 1 , sup t d 2 dt 2 γ(t) C 2 <K 2
is invariant for suitably chosen K 1 and K 2 greater than zero under the mapping

(Lγ)(t) = L t (γ(t)) for all t ∈ [-ε 2 , ε 2 ]
naturally defined by the family of operators on a curve.

Proposition 1.3. Choose C 1 , C 2 , C 3 , C 4 > 0 so that Proposition 1.2 holds.
Then there is ε 2 > 0 sufficiently small and K 1 , K 2 > 0 such that In Section 4, we prove point (iii) of Theorem 1.2 by calculating a linear response formula. To ease the notation, we denote P t := L t,ρ(t) .

L Ĉε 2 ,K 1 ,K 2 ⊂ Ĉε 2 ,K 1 ,K 2 . Since Ĉε 2 ,K 1 ,K 2 is relatively compact in C 1 ([-ε 2 , ε 2 ], C 1 (S 1 , R)), for any γ ∈ Ĉ, there is an increasing sequence {n k } k∈N and γ ∈ C 1 ([-ε 2 , ε 2 ], C 2 (S 1 , R)) such that L n k γ → γ,
Proposition 1.4. We have the formula

∂ t ρ| t=0 = -(1 -P ) -1 P (ρ(0)A ρ(0) ) ′ , (9)
where ∂ t is understood as

∂ t : C 1 ([-ε 3 , ε 3 ], C 1 (S 1 , R)) → C 0 ([-ε 3 , ε 3 ], C 1 (S 1 , R)).
More generally, for any t ∈ (-ε 3 , ε 3 ) where ε 3 ≤ ε 2 as in Proposition 1.3,

∂ t ρ| t= t = -(1 -P t + tP K t) -1 P (K t(ρ( t)))
where P t = P Q tρ( t) and

K t (g) = ρ(t) Φ ′ t,ρ(t) A g • Φ -1 t,ρ(t) ′ , g ∈ C 1 (S 1 , R). (10)
The formula in ( 9) is reminiscent of the one obtained in [4, Theorem 2.2] for perturbations of expanding circle maps. This is not surprising, since for t approaching 0 the operator is very close to linear. At the end of the paper, Appendix A gathers some elementary formulae used throughout our calculations.

Invariant Densities for the Self-Consistent Transfer Operators

We now show that the fixed point found above is unique using an invariant cone argument and the Hilbert metric. We remind the reader of our defintion

V a := ϕ : S 1 → R + : ϕ(x) ϕ(y) ≤ exp(a|x -y|) .
and of the Hilbert metric

θ a (ϕ, ψ) = log 1 α a (ϕ, ψ)α a (ψ, ϕ)
where

α a (ψ, ϕ) = inf ϕ(x) ψ(x) , e a|x-y| ϕ(x) -ϕ(y) e a|x-y| ψ(x) -ψ(y) : x = y . (11)
It is a known important fact that P is a contraction with respect to this metric (see e.g. [START_REF] Liverani | Decay of correlations[END_REF][START_REF] Viana | Stochastic dynamics deterministic systems[END_REF]). In the following we show that for suitable values of t and a > 0, the operator L t is also a contraction. The lemmas below gather a few useful inequalities.

Lemma 2.1. There is ε 0 > 0 sufficiently small such that for every a > 0, every t ∈ (-ε 0 , ε 0 ), and ϕ, ψ ∈ C 0 with ϕ = ψ = 1,

Φ -1 t,ψ (x) -Φ -1 t,ϕ (x) ≤ K|t| ψ -ϕ C 0 (12)
Proof. Notice that

|Φ t,ψ (x) -Φ t,ϕ (x)| = |t| ∂ 1 h(x, y)(ψ(y) -ϕ(y)) ≤ |t|K ψ -ϕ C 0
By the mean value theorem, for any x

min |Φ ′ t,ψ | ≤ |Φ t,ψ (Φ -1 t,ψ (x)) -Φ t,ψ (Φ -1 t,ϕ (x))| |Φ -1 t,ϕ (x) -Φ -1 t,ψ (x)| = |Φ t,ϕ (Φ -1 t,ϕ (x)) -Φ t,ψ (Φ -1 t,ϕ (x))| |Φ -1 t,ϕ (x) -Φ -1 t,ψ (x)| ≤ K|t| ψ -ϕ C 0 |Φ -1 t,ϕ (x) -Φ -1 t,ψ (x)| .
Since by assumption min |Φ ′ t,ϕ | ≥ 1 -K|t| > 0 for any ϕ, (12) follows.

Lemma 2.2. There is ε 0 > 0 sufficiently small such that for every a > 0 there is K a > 0 such that for every t ∈ (-ε 0 , ε 0 ) and ϕ, ψ ∈ V a with ϕ = ψ = 1, the following holds:

Φ -1 t,ψ (x) -Φ -1 t,ϕ (x) ≤ K a |t| max 1 - ϕ ψ (13) Φ -1 t,ϕ (x) -Φ -1 t,ϕ (y) ≤ (1 + K|t|)|x -y| (14) max 1 - ϕ ψ ≥ K a log α(ϕ, ψ) (15)
Proof. The proof of ( 13) follows along the same lines of the proof of Lemma 2.1 after noticing that

∂ 1 h(x, y)(ψ(y) -ϕ(y)) ≤ max ψ • K • max 1 - ϕ ψ ≤ e a 2 K • max 1 - ϕ ψ .
Since ψ ∈ V a and has integral equal to one, max ψ ≤ e a 2 . To prove [START_REF] Galatolo | A linear response for dynamical systems with additive noise[END_REF], notice that for t sufficiently small

|(Φ -1 t,ϕ ) ′ | < 1 + Kt. To prove (15), pick x 0 ∈ T such that max |1 -ϕ ψ | = 1 -ϕ ψ (x 0 ). Then the inequality is implied by log α(ϕ, ψ) ≤ log ϕ ψ (x 0 ) = log 1 -1 - ϕ ψ (x 0 ) ≤ K -1 a 1 - ϕ ψ (x 0 )
as ϕ/ψ is lower bounded by e -a 2 . Analogous conclusion holds if max x∈T |1 -ϕ ψ (x)| = ϕ ψ (x 0 ) -1 for some x 0 ∈ T.

We are now ready to prove Proposition 1.1.

Proof of Proposition 1.1.

Step 1: There is λ < 1 such that L t V a ⊂ V λa First of all, we show that there is K such that Q t V a ⊂ V a+K(a+1)|t| for t sufficiently small and a > 0. In fact, we show more. Pick

ε 0 > 0 such that 1 -|t| max |∂ 1 h| > 0 for t ∈ [-ε 0 , ε 0 ]. For any ϕ, ψ ∈ V a we have ϕ • Φ -1 t,ψ (x) ϕ • Φ -1 t,ψ (y)
≤ e a(1+K|t|) |x -y| and

(1 + t∂ x A ψ ) • Φ -1 t,ψ (y) (1 + t∂ x A ψ ) • Φ -1 t,ψ (x) ≤ e K|t||x-y| ,
where the bound on the first expression follows from ϕ ∈ V a and ( 14), while that on the second expression follows from the smoothness of h (in particular K depends on the norms of the derivatives of h). Thus

Q t,ψ ϕ(x) Q t,ψ ϕ(y) = ϕ • Φ -1 t,ψ (x) ϕ • Φ -1 t,ψ (y) (1 + t∂ x A ψ ) • Φ -1 t,ψ (y) (1 + t∂ x A ψ ) • Φ -1 t,ψ (x)
≤ e (a+(a+1)K|t|)|x-y| .

Analogously, it is easy to check (by noticing that

f ′ (x) f ′ (y) ≤ e 1 ω max |f ′′ ||x-y| ) that if ϕ ∈ V a then P ϕ ∈ V a/ω+max |f ′′ |/ω 2 .
Picking

a 0 := max |f ′′ | ω(ω -1) one has that for any a > a 0 , a/ω + max |f ′′ |/ω 2 < a.
Therefore, for any a > a 0 there is ε 0 > 0 and a 1 , a 2 satisfying

a < a + (a + 1)K|ε 0 | =: a 1 < a 2 < ω a - max |f ′′ | ω 2 .
With this choice of a 1 and a 2 we have:

Q t V a ⊂ V a 1 ⊂ V a 2 for t ∈ (-ε 0 , ε 0 ), and P V a 2 ⊂ V a 2 ω + max |f ′′ | ω 2 ⊂ V a . More precisely, L t V a ⊂ V λa with λ < 1 depending on a. § .
From now on, a > a 0 is considered to be fixed.

Step 2: for t small, L t : V a → V a is a contraction with respect to θ a Let us start by showing how assuming

θ a 2 (Q t ϕ, Q t ψ) ≤ [1 + O(|t|)] θ a (ϕ, ψ) (16)
yields the proof of the claim in Step 2. Recall that with the choice of a and a 2 as in

Step 1, P : (V a 2 , θ a 2 ) → (V a , θ a ) is a contraction as implied by the linearity of P and the fact that V a ⊂ V a 2 . Call σ < 1 the contraction rate, i.e. θ a (P ϕ, P ψ) ≤ σθ a 2 (ϕ, ψ) for any ϕ, ψ ∈ V a 2 . This and ( 16) imply that θ a (L t ϕ,

L t ψ) ≤ σ [1 + O(t)] θ a (ϕ, ψ).
Picking |t| sufficiently small, equation ( 16) implies that L t is also a contraction. Now we proceed with the proof of [START_REF] Galatolo | Quadratic response of random and deterministic dynamical systems[END_REF]. Notice that by the triangle inequality

θ a 2 (Q t ϕ, Q t ψ) ≤ θ a 2 (Q t,ϕ ϕ, Q t,ϕ ψ) + θ a 2 (Q t,ϕ ψ, Q t,ψ ψ).
For the first term, the operator Q t,ϕ is linear and maps the cone V a to V a 1 . V a 1 is strictly inside the cone V a 2 . Since Q t,ϕ (V a ) is strictly contained in V a 2 and its diameter in (V a 2 , θ a 2 ) is finite [34, Proposition 2.5], the mapping Q t,ψ : (V a , θ a ) → (V a 2 , θ a 2 ) is a contraction [START_REF] Viana | Stochastic dynamics deterministic systems[END_REF]. Therefore,

θ a 2 (Q t,ϕ ϕ, Q t,ϕ ψ) ≤ θ a (ϕ, ψ).
For the second term, computations are lengthier and occupy the rest of the proof. Recall that

θ a 2 (Q t,ϕ ψ, Q t,ψ ψ) = log 1 α a 2 (ϕ, ψ)α a 2 (ψ, ϕ)
where α is the infimum of two expressions as defined in [START_REF] Boyarsky | Laws of chaos: invariant measures and dynamical systems in one dimension[END_REF]. For the first expression in the definition of α a 2 :

Q t,ϕ ψ(x) Q t,ψ ψ(x) = ψ • Φ -1 t,ϕ (x) ψ • Φ -1 t,ψ (x) • (1 + t∂ x A ψ ) • Φ -1 t,ψ (x) (1 + t∂ x A ϕ ) • Φ -1 t,ϕ (x) = = ψ • Φ -1 t,ϕ (x) ψ • Φ -1 t,ψ (x) • (1 + t∂ x A ψ ) • Φ -1 t,ψ (x) (1 + t∂ x A ψ ) • Φ -1 t,ϕ (x) • (1 + t∂ x A ψ ) • Φ -1 t,ϕ (x) (1 + t∂ x A ϕ ) • Φ -1 t,ϕ (x)
The bound on the first two factors follows from ψ, ϕ ∈ V a , (1 + tA ψ ) ∈ V K|t| (again K depends on bounds on h and its derivatives), and from inequality (13) that give

ψ • Φ -1 t,ϕ (x) ψ • Φ -1 t,ψ (x)
• (1 + t∂ x A ψ ) • Φ -1 t,ψ (x) (1 + t∂ x A ψ ) • Φ -1 t,ϕ (x) ≥ e -(a+K|t|)|Φ -1 t,ϕ (x)-Φ -1 t,ψ (x)| ≥ e -aKa|t| max 1-ϕ ψ .
For the third factor

(1 + t∂ x A ψ )(x) (1 + t∂ x A ϕ )(x) = 1 + t ∂ 1 h(x, y)ψ(y)dy 1 + t ∂ 1 h(x, y)ϕ(y)dy ≥ 1 + t ∂ 1 h(x, y)ϕ(y)dy + t min{ ψ ϕ } -1 |∂ 1 h(x, y)|ϕ(y)dy 1 + t ∂ 1 h(x, y)ϕ(y)dy ≥ 1 -K|t| max ψ ϕ -1 (17)
≥ α a (ϕ, ψ) K|t| .

Putting all the estimates together and using (15) one obtains

Q t,ϕ ψ(x) Q t,ψ ψ(x) ≥ α a (ϕ, ψ) Ka|t| .
where K a > 0 indicates a constant depending on h and on a only.

For the second expression in the definition of α we need to estimate

e a 2 |x-y| Q t,ϕ ψ(x) -Q t,ϕ ψ(y) e a 2 |x-y| Q t,ψ ψ(x) -Q t,ψ ψ(y) = e a 2 |x-y| ψ (1+tAϕ) • Φ -1 t,ϕ (x) - ψ (1+tAϕ) • Φ -1 t,ϕ (y) e a 2 |x-y| ψ (1+tA ψ ) • Φ -1 t,ψ (x) - ψ (1+tA ψ ) • Φ -1 t,ψ (y) = A • B
where

A := e a 2 |x-y| ψ (1+tAϕ) • Φ -1 t,ϕ (x) - ψ (1+tAϕ) • Φ -1 t,ϕ (y) e a 2 |x-y| ψ (1+tAϕ) • Φ -1 t,ψ (x) - ψ (1+tAϕ) • Φ -1 t,ψ (y)
B := e a 2 |x-y| ψ (1+tAϕ) • Φ -1 t,ψ (x) - ψ (1+tAϕ) • Φ -1 t,ψ (y) e a 2 |x-y| ψ (1+tA ψ ) • Φ -1 t,ψ (x) - ψ (1+tA ψ ) • Φ -1 t,ψ (y) A ≥ e a 2 |x-y| ψ (1+tAϕ) • Φ -1 t,ψ (x)e -a 1 |Φ -1 t,ψ (x)-Φ -1 t,ϕ (x)| - ψ (1+tAϕ) • Φ -1 t,ψ (y)e a 1 |Φ -1 t,ψ (y)-Φ -1 t,ϕ (y)| e a 2 |x-y| ψ (1+tAϕ) • Φ -1 t,ψ (x) - ψ (1+tAϕ) • Φ -1 t,ψ (y) ≥ e a 2 |x-y| ψ (1+tAϕ) • Φ -1 t,ψ (x)e -K ′ a |t| max 1-ϕ ψ - ψ (1+tAϕ) • Φ -1 t,ψ (y)e K ′ a |t| max 1-ϕ ψ e a 2 |x-y| ψ (1+tAϕ) • Φ -1 t,ψ (x) - ψ (1+tAϕ) • Φ -1 t,ψ (y) (18)
≥ e

K ′ a |t| max 1-ϕ ψ   1 + e -2K ′ |t| max 1-ϕ ψ -1 e a 2 |x-y| ψ (1+tAϕ) • Φ -1 t,ψ (x) - ψ (1+tAϕ) • Φ -1 t,ψ (y)   (19)
≥ e

K ′′ a |t| max 1-ϕ ψ ≥ α a (ϕ, ψ) Ka|t| (20)
where [START_REF] Jiang | Linear response function for coupled hyperbolic attractors[END_REF] follows from ( 13), [START_REF] Kloeckner | The linear request problem[END_REF] follows from the fact that the denominator in ( 19) is bounded since ϕ (1+tAϕ) • Φ -1 t,ψ belongs to V a 1 with a 1 < a 2 , and K a , K ′ a and K ′′ a are constants depending on a. The factor B can be bounded in a similar way:

B ≥ α a (ϕ, ψ) Ka|t| .
Putting together all the previous estimates one obtains

α a 2 (Q t,ϕ ψ, Q t,ψ ψ) ≥ α a (ϕ, ψ) 2Ka|t| ,
and analogously

α a 2 (Q t,ψ ψ, Q t,ϕ ψ) ≥ α a (ψ, ϕ) 2Ka|t| which implies θ a 2 (Q t,ϕ ψ, Q t,ψ ψ) = log 1 α a 2 (Q t,ϕ ψ, Q t,ψ ψ)α a 2 (Q t,ψ ψ, Q t,ϕ ψ) ≤ 2K a |t|θ a (ϕ, ψ).
We are now ready to conclude the proof of Theorem 1.1

Proof of Theorem 1.1. The proof can be finished by standard arguments. The existence of an invariant density ρ(t) in Lip(S 1 , R + ) follows for example by the argument presented in [34, Section 2.3] and the exponential convergence of L n t ϕ to ρ(t) in the C 0 norm, hence uniqueness of ρ(t) can be argued as in [START_REF] Viana | Stochastic dynamics deterministic systems[END_REF]Section 2.4].

From now on let's assume that f ∈ C 5 (S 1 , S 1 ) and h ∈ C 5 (S 1 × S 1 , R), and prove point (i) of Theorem 1.2, i.e. we show that in this case ρ(t) ∈ C 3 (S 1 , R + ). Recall that • C 0 denotes the supremum norm on C 0 (S 1 , R), and

C C := ϕ ∈ C 4 (S 1 , R + ) : ϕ = 1, d i ϕ dx i C 0 < C i i = 1, ..., 4 .
Proof of Proposition 1.2. Let t ∈ [-ε 0 , ε 0 ] and ϕ ∈ C 4 (S 1 , R + ). We remind the reader that F t,ϕ is defined as F t,ϕ = f • Φ t,ϕ . By (43)-(45) of the Appendix and some elementary bounds, we obtain

|F ′ t,ϕ | ≥ ω(1 -K 0 |t|), (21)
|F ′′ t,ϕ | ≤ c((1 + K 0 |t|) 2 + K 1 )|t|, (22)
|F ′′′ t,ϕ | ≤ c((1 + K 0 |t|) 3 + (1 + K 0 |t|)K 1 |t| + K 2 |t|), (23)
where

K i = max |h (i) | and c = max{max |f ′′′ |, 3 max |f ′′ |, max |f ′ |}.
We now choose the constants C i , i = 1, . . . , 4, starting with the choice of C 1 . We denote the inverse branches of F t,ϕ by F -1 t,ϕ,i , i = 1, . . . , N .

(L t ϕ) ′ C 0 = N i=1 ϕ F ′ t,ϕ • F -1 t,ϕ,i ′ C 0 = N i=1 ϕ ′ (F ′ t,ϕ ) 2 - ϕF ′′ t,ϕ (F ′ t,ϕ ) 3 • F -1 t,ϕ,i C 0 ≤ N ϕ ′ (F ′ t,ϕ ) 2 - ϕF ′′ t,ϕ (F ′ t,ϕ ) 3 C 0 ≤ N 1 ω 2 (1 -K 0 t) 2 ϕ ′ C 0 + max |f ′′ |(1 + K 0 t) 2 + cK 1 t ω 3 (1 -K 0 t) 3 ϕ C 0
by the bounds ( 21)-( 22). Since ϕ is continuous and has unit integral, we have

ϕ C 0 ≤ 1 + ϕ ′ C 0 giving us (L t ϕ) ′ C 0 ≤ N 1 ω 2 (1 -K 0 t) 2 + max |f ′′ |(1 + Kt) 2 + cK 1 t ω 3 (1 -K 0 t) 3 ϕ ′ C 0 + N (max |f ′′ |(1 + K 0 t) 2 + cK 1 t) ω 3 (1 -K 0 t) 3 =: σ 1 ϕ ′ C 0 + R 1 . ( 24
)
By the assumption ( 8) on f , we can choose t small enough (by potentially decreasing

ε 0 ) such that σ 1 < 1. Choose C 1 such that R 1 1 -σ 1 ≤ C 1 ,
and then

ϕ ′ C 0 ≤ C 1 implies (L t ϕ) ′ C 0 ≤ C 1 . We move on to the choice of C 2 . (L t ϕ) ′′ C 0 = N i=1 ϕ F ′ t,ϕ • F -1 t,ϕ,i ′′ C 0 ≤ N ϕ ′′ (F ′ t,ϕ ) 3 -3 ϕ ′ F ′′ t,ϕ (F ′ t,ϕ ) 4 - ϕF ′′′ t,ϕ (F ′ t,ϕ ) 4 + 3 ϕ(F ′′ t,ϕ ) 2 (F ′ t,ϕ ) 5 C 0 ≤ N 1 ω 3 (1 -K 0 t) 3 ϕ ′′ C 0 + N 3 F ′′ t,ϕ (F ′ t,ϕ ) 4 + F ′′′ t,ϕ (F ′ t,ϕ ) 4 -3 (F ′′ t,ϕ ) 2 (F ′ t,ϕ ) 5 C 0 ϕ ′ C 0 + N F ′′′ t,ϕ (F ′ t,ϕ ) 4 -3 (F ′′ t,ϕ ) 2 (F ′ t,ϕ ) 5 C 0 =: σ 2 ϕ ′′ C 0 + R (1) 2 ϕ ′ C 0 + R (2) 2 . ( 25
)
For sufficiently small t we have σ 2 < 1 by assumption [START_REF] Baladi | Alternative proofs of linear response for piecewise expanding unimodal maps[END_REF] and by ( 21)-( 23)

R (1) 2 , R (2)
2 are bounded uniformly in ϕ and in t if t is in some fixed neighborhood of zero, so choosing R

(1)

2 C 1 + R (2) 2 1 -σ 2 ≤ C 2
will be sufficient. Similar calculations can be done for the third and fourth derivative. We can conclude that

(L t ϕ) ′′′ C 0 ≤ σ 3 ϕ ′′′ C 0 + R (1) 3 ϕ ′′ C 0 + R (2) 3 ϕ ′ C 0 + R (3) 3 (L t ϕ) ′′′′ C 0 ≤ σ 4 ϕ ′′′′ C 0 + R (1) 4 ϕ ′′′ C 0 + R (2) 4 ϕ ′′ C 0 + R (3) 4 ϕ ′ C 0 + R (4) 4
where σ 3 , σ 4 < 1 for sufficiently small t and R (j) i only depend on the derivatives of f and h and on t (such that they are bounded uniformly in t if t is in some fixed neighborhood of zero). Choosing

R (1) 3 C 2 + R (2) 3 C 1 + R (3) 3 1 -σ 3 ≤ C 3 R (1) 4 C 3 + R (2) 4 C 2 + R (3) 4 C 1 + R (4) 4 1 -σ 4 ≤ C 4
will be sufficient for our purposes.

Lemma 2.3. If f ∈ C 5 (S 1 , S 1 ) and h ∈ C 5 (S 1 × S 1 , R), then for 1 ≤ k ≤ 3, Q t : C C → C C is continuous with respect to the C k distance.
Proof. Fix k ∈ {1, 2, 3}. For any ϕ, ψ ∈ C k (S 1 , R), using expression ( 6) for Q t one obtains

Q t ϕ -Q t ψ = ϕ -ψ 1 + tA ′ ϕ • Φ -1 t,ϕ (x)+ + ψ 1 + tA ′ ϕ - ψ 1 + tA ′ ψ • Φ -1 t,ϕ (x)+ + ψ 1 + tA ′ ψ • Φ -1 t,ϕ (x) - ψ 1 + tA ′ ψ • Φ -1 t,ψ (x).
For t sufficiently, small 1 1+tA ′ ϕ and Φ -1 t,ϕ belong to C 4 . Then it follows that

ϕ -ψ 1 + tA ′ ϕ • Φ -1 t,ϕ (x) C k ≤ M 1 1 + tA ′ ϕ • Φ -1 t,ϕ (x) C k (ϕ -ψ) • Φ -1 t,ϕ (x) C k ≤ M ′ 1 1 + tA ′ ϕ C k max ℓ=1,...,k Φ -1 t,ϕ ℓ C k 2 ϕ -ψ C k ≤ K ϕ -ψ C k ,
where the first two inequalities hold by ( 42) and (47) of the Appendix, respectively. The third inequality can be proved by using (52) of the Appendix and some standard bounds on the derivatives of Φ t,ϕ . To bound the second term, notice that

1 1 + tA ′ ϕ - 1 1 + tA ′ ψ = t A ′ ψ -A ′ ϕ (1 + tA ′ ϕ )(1 + tA ′ ψ )
and

A ′ ψ -A ′ ϕ C k ≤ K ψ -ϕ C 0 .
Then similarly to the calculation above,

ψ 1 + tA ′ ϕ - ψ 1 + tA ′ ψ • Φ -1 t,ϕ (x) C k ≤ M t ψ (1 + tA ′ ϕ )(1 + tA ′ ψ ) • Φ -1 t,ϕ (x) C k (ϕ -ψ) • Φ -1 t,ϕ (x) C k ≤ M ′ t ψ (1 + tA ′ ϕ )(1 + tA ′ ψ ) C k max ℓ=1,...,k Φ -1 t,ϕ ℓ C k 2 ϕ -ψ C k ≤ K C ϕ -ψ C k ,
where the last inequality also exploits that ψ ∈ C C . Using (38)-( 41) of the Appendix we can compute that there exists an explicit polynomial function p such that

ψ 1 + tA ′ ψ • Φ -1 t,ϕ (x) - ψ 1 + tA ′ ψ • Φ -1 t,ψ (x) C k ≤ p ψ 1 + tA ′ ψ C k+1 , Φ -1 t,ϕ C k , Φ -1 t,ψ C k ϕ -ψ C k . Since ψ 1+tA ′ ψ ∈ C 4
, simple calculations as before give us

ψ 1 + tA ′ ψ • Φ -1 t,ϕ (x) - ψ 1 + tA ′ ψ • Φ -1 t,ψ (x) C k ≤ K C ϕ -ψ C k .
By triangle inequality one can bound Q t ϕ -Q t ψ C k and the thesis follows.

We can finally prove that ρ(t

) ∈ C 3 (S 1 , R). Proof of Theorem 1.2 Point (i). C C ⊂ C 4 (S 1 , R) is convex and relatively compact in C 3 (S 1 , R).
The invariance of the set C C under the action of L t , as claimed in Proposition 1.2 and the continuity of L t restricted to C C as stated in Lemma 2.3 implies by the Schauder fixed point theorem that L t has a fixed density in C 3 (S 1 , R). By Theorem 1.1 this density must be the unique Lipschitz density ρ(t).

Smooth Dependence of the Invariant Densities on the Coupling Strength

Given t ∈ [-ε 1 , ε 1 ], the relation t → ρ(t) defines a curve of fixed densities, that is, for each t the density ρ(t) is such that

L t ρ(t) = ρ(t). Recall the definition of Ĉε,K 1 ,K 2 Ĉε,K 1 ,K 2 := γ : [-ε, ε] → C C s.t. (t,x) →γ(t)(x) is C 2 sup t d dt γ(t) C 3 <K 1 , sup t d 2 dt 2 γ(t) C 2 <K 2 and the operator L (Lγ)(t) = L t (γ(t)) for all t ∈ [-ε, ε] for γ ∈ Ĉε,K 1 ,K 2 .
The strategy consists of (i) showing that there are ε 2 , K 1 , and K 2 such that L keeps Ĉε 2 ,K 1 ,K 2 invariant, and (ii) using Proposition 1.1 to show that this implies that t → ρ(t) is C 1 .

Proof of Proposition 1.3. Let us start by noticing that the implicit function theorem implies that the function g : [-ε 1 , ε 1 ] × S 1 → S 1 , defined as g(t, y) := (Φ t,ϕ ) -1 (y), is differentiable and C 5 . Furthermore, its C 5 norm can be uniformly bounded with respect to ϕ

¶ . Let ε ≤ ε 1 and pick γ ∈ Ĉε,K 1 ,K 2 . d dt (Lγ)(t) = d ds L s s=t γ(t) + lim δ→0 L t γ(t + δ) -L t γ(t) δ . ( 26
)
¶ This is a consequence of the fact that the function to which we apply the implicit function theorem, F (t, y, x) = Φt,ϕ(x) -y, has C 5 norm uniformly bounded varying ϕ.

For the first term in the right hand side of ( 26)

d ds L s s=t ϕ = P lim δ→0 Q t+δ -Q t δ ϕ = P d ds ϕ ∂ x Φ s,ϕ • Φ -1 s,ϕ s=t = P ϕ ′ ∂ x Φ t,ϕ • Φ -1 t,ϕ ∂ t Φ -1 t,ϕ - ϕ (∂ x Φ t,ϕ ) 2 • Φ -1 t,ϕ ∂ t [∂ x Φ t,ϕ • Φ -1 t,ϕ ] = P ϕ ′ ∂ x Φ t,ϕ • Φ -1 t,ϕ ∂ t Φ -1 t,ϕ - ϕ (∂ x Φ t,ϕ ) 2 • Φ -1 t,ϕ [∂ t ∂ x Φ t,ϕ • Φ -1 t,ϕ + ∂ 2 x Φ t,ϕ • Φ -1 t,ϕ ∂ t Φ -1 t,ϕ ] . (27)
Therefore, putting ϕ = γ(t) and recalling that γ(t) ∈ C C , for t sufficiently small all the functions above have bounded norm in C 3 (S 1 , R), and therefore there is a constant K depending on f, h and C such that

d ds L s s=t γ(t) C 3 ≤ K. (28)
The second term in the RHS of (26) equals

Ξ(t) := P lim δ→0 Q t γ(t + δ) -Q t γ(t) δ = P lim δ→0 1 δ γ(t + δ) 1 + tA ′ γ(t+δ) • Φ -1 t,γ(t+δ) - γ(t) 1 + tA ′ γ(t) • Φ -1 t,γ(t) = P lim δ→0 γ(t + δ) -γ(t) δ(1 + A ′ γ(t+δ) )
• Φ -1 t,γ(t+δ)

+ P lim δ→0 1 δ γ(t) 1 1 + tA ′ γ(t+δ) - 1 1 + tA ′ γ(t) • Φ -1 t,γ(t+δ) + P lim δ→0 1 δ γ(t) 1 + tA ′ γ(t) • Φ -1 t,γ(t+δ) - γ(t) 1 + tA ′ γ(t) • Φ -1 t,γ(t) , giving Ξ(t) = = P Q t,γ(t) d dt γ(t) (29)
+ P γ(t) d ds

1 1 + tA ′ γ(s) s=t • Φ -1 t,γ(t) + ∂ x Q t,γ(t) γ(t) d ds Φ -1 t,γ(s) s=t (30)
For the first term of the sum, one can argue in a way similar to the proof of Proposition 1.2 and obtain the inequality

P Q t,γ(t) d dt γ(t) C 3 ≤ σ d dt γ(t) C 3 + K (31)
with σ < 1 and K > 0 (in general different from the K above). To see this, notice that in the referenced calculations we use uniform bounds for F (k) t,ϕ in ϕ, only exploiting its unit integral. Furthermore, since d dt γ(t) = 0, the bound

ϕ C 0 ≤ 1 + ϕ ′ C 0 simplifies to d dt γ(t) C 0 ≤ ( d dt γ(t)) ′ C 0 .
The regularity assumptions allow to deduce that the C 3 norm of the other terms is uniformly bounded. With computations similar to those in the proof of Proposition 1.2, one can deduce there are σ < 1 (by possibly decreasing ε 1 ) and K > 0 such that Ξ C 3 ≤ σ d dt γ C 3 + K. Putting this estimate together with the one in (28) we get

d dt (Lγ)(t) C 3 ≤ d ds L s s=t γ(t) C 3 + Ξ(t) C 3 ≤ σ d dt γ(t) C 3 + K.
We conclude that there is a K 1 sufficiently large such that if The first term corresponds to the derivative with respect to t of ( 27) when ϕ = γ(t), which in turn is a combination of the functions

∂ x γ(t) ∂ x Φ t,γ(t) • Φ -1 t,γ(t) , ∂ s Φ -1 s,γ(t) | s=t , γ(t) (∂ x Φ t,γ(t) ) 2 • Φ -1 t,γ(t) , ∂ s ∂ x Φ s,γ(t) | s=t • Φ -1 t,γ(t) , ∂ 2 x Φ t,γ(t) • Φ -1 t,γ(t) .
With the given regularity hypotheses all the derivatives with respect to t of the above terms are continuous in the variable x. In particular, for the first term notice that the assumed regularity of all the terms implies

d dt ∂ x γ(t) ∂ x Φ t,γ(t) = ∂ x d dt γ(t) ∂ x Φ t,γ(t)
.

We conclude that there exists a K > 0 such that d dt d ds L s s=t γ(t) C 2 < K. As for d dt Ξ(t), consider the decomposition of Ξ(t) given by ( 29)- [START_REF] Ruelle | General linear response formula in statistical mechanics, and the fluctuation-dissipation theorem far from equilibrium[END_REF]. The derivative of

P Q t,γ(t) d dt γ(t) is d ds P Q s,γ(s) d ds γ(s) s=t = P lim δ→0 (Q t+δ,γ(t+δ) -Q t,γ(t+δ) ) d dt γ(t + δ) δ + P lim δ→0 (Q t,γ(t+δ) -Q t,γ(t) ) d dt γ(t + δ) δ + P lim δ→0 Q t,γ(t) ( d dt γ(t + δ) -d dt γ(t)) δ .
The first term in the sum on the RHS can be treated with computations similar to those leading to [START_REF] Majda | Challenges in climate science and contemporary applied mathematics[END_REF], and equals

P lim δ→0 (Q t+δ,γ(t) -Q t,γ(t) ) d dt γ(t) δ = P ∂ x dγ(t) dt ∂ x Φ t,γ(t) • Φ -1 t,γ(t) ∂ s Φ -1 s,γ(t) | s=t - dγ(t) dt (∂ x Φ t,γ(t) ) 2 • Φ -1 t,γ(t) ∂ s ∂ x Φ s,γ(t) | s=t • Φ -1 t,γ(t) + ∂ 2 x Φ t,γ(t) • Φ -1 t,γ(t) ∂ s Φ -1 s,γ(t) | s=t
which can be verified to have uniformly bounded C 2 norm.

The second term can be expressed as

P lim δ→0 (Q t,γ(t+δ) -Q t,γ(t) ) d dt γ(t + δ) δ = d ds dγ dt (t) ∂ x Φ t,γ(s) • Φ -1
t,γ(s) s=t [START_REF] Ruelle | Structure and f-dependence of the acim for a unimodal map f of Misiurewicz type[END_REF] and noticing that

d ds 1 ∂ x Φ t,γ(s) = - 1 (∂ x Φ t,γ(s) ) 2 d ds ∂ x Φ t,γ(s) , d ds dγ(t) dt , and
∂ s Φ -1 t,γ(s)
have uniformly bounded C 2 norm, it follows from straightforward computations that the whole of (32) has bounded norm. ‖ For the third term we can find an inequality with a similar argument that leads to [START_REF] Ruelle | A review of linear response theory for general differentiable dynamical systems[END_REF]:

P Q t,γ(t) d 2 dt 2 γ(t) C 2 ≤ σ d 2 dt 2 γ(t) C 2 + K
For what concerns the derivative of [START_REF] Ruelle | General linear response formula in statistical mechanics, and the fluctuation-dissipation theorem far from equilibrium[END_REF], it involves the same terms listed above. Most of the factors appearing have been previously treated apart from d ds

1 1 + tA ′ γ(s) s=t and ∂ x Q t,γ(t) γ(t) .
For the first one is easy to see that all the derivatives in t of 1 + tA ′ γ(s) have bounded C 2 norm in x. For the second one, it follows from previous computations that Q t,γ(t) γ(t) is at least C 1 ((-ε, ε), S 1 ) and this implies

d dt ∂ x Q t,γ(t) γ(t) = ∂ x d dt Q t,γ(t) γ(t). d dt Q t,γ (
t) γ(t) can be shown to have uniformly bounded C 3 norm in the variable x with the same computations used to control [START_REF] Mackay | Management of complex dynamical systems[END_REF], and therefore

∂ x d dt Q t,γ(t) γ(t) C 2 ≤ K. Therefore, d 2 dt 2 (Lγ)(t) C 2 ≤ σ d 2 dt 2 γ(t) C 2
+ K so we can find K 2 sufficiently large so that d dt γ(t) C 3 < K 1 , with K 1 as above, and

d 2 dt 2 γ(t) C 2 < K 2 implies d 2 dt 2 (Lγ)(t) C 2 < K 2 .
‖ Notice that in the expression above the factor ∂x dγ dt appears, and it is for this kinds of factors that the choice of controlling two different norms for dγ dt and d 2 γ dt 2 becomes crucial.

Proof of Theorem 1.2 Point (ii). Given any γ ∈ Ĉε

2 ,K 1 ,K 2 , it follows from Proposi- tion 1.3 that {L n γ} n∈N 0 ∈ Ĉε 2 ,K 1 ,K 2 . Since Ĉε 2 ,K 1 ,K 2 ⊂ C 2 ([-ε 2 , ε 2 ], C 2 (S 1 , R)) is re- latively compact in C 1 ([-ε 2 , ε 2 ], C 1 (S 1 , R)) * * , there is a strictly increasing sequence n k and γ ∈ C 1 ([-ε 2 , ε 2 ], C 1 (S 1 , R)) such that L n k γ → γ in the C 1 topology.
In particular, for every t, L n k t γ(t) -γ(t) C 0 → 0 for k → ∞. However, from Proposition 1.1 we know that for every t, L n t γ(t)-ρ(t) C 0 → 0 for n → ∞ from which we deduce that γ(t) = ρ(t), and therefore t → ρ(t) belongs to C

1 ([-ε 2 , ε 2 ], C 1 (S 1 , R)).

Linear Response Formula

In this section we derive a linear response formula for the curve t → ρ(t) of fixed densities. We proceed to prove the proposition in an analogous way to [4, Theorem 2.2] and apply perturbation arguments to deal with the self-consistency.

Lemma 4.1. Considering {P t } t∈[-ε 2 ,ε 2 ] as bounded operators from C 2 (S 1 , R) to C 1 (S 1 , R), we claim that the map t → P t is differentiable at t ∈ (-ε 2 , ε 2 ) and ∂ t P t | t= t(ϕ) = -P ϕ Φ ′ t,ρ( t) (A ρ( t) + tA ∂tρ(t)| t= t ) • Φ -1 t,ρ( t) ′ .
Proof. Recall that P t = P Q t,ρ(t) which implies that

∂ t P t = P ∂ t Q t,ρ(t) . (33)
Recall also that for t

∈ [-ε 0 , ε 0 ], Φ t,ρ(t) is a diffeomorphism of S 1 and t → Φ t,ρ(t) is C 2 . Let ϕ ∈ C 2 . Q t,ρ(t) (ϕ) -Q t,ρ( t) (ϕ) = = ϕ Φ ′ t,ρ(t) • Φ -1 t,ρ(t) - ϕ Φ ′ t,ρ( t) • Φ -1 t,ρ( t) = -(t -t)(∂ t Φ t,ρ(t) | t=t * • Φ -1 t * ,ρ(t * ) ) ′ ϕ Φ ′ t * ,ρ(t * ) • Φ -1 t * ,ρ(t * ) -(t -t)∂ t Φ t,ρ(t) | t=t * • Φ -1 t * ,ρ(t * ) ϕ ′ (Φ ′ t * ,ρ(t * ) ) 2 • Φ -1 t * ,ρ(t * ) - ϕ • Φ ′′ t * ,ρ(t * ) (Φ ′ t * ,ρ(t * ) ) 3 • Φ -1 t * ,ρ(t * ) .
for some t * ∈ (t, t). Notice that

(∂ t Φ t,ρ(t) | t=t * • Φ -1 t * ,ρ(t * ) ) ′ ϕ Φ ′ t * ,ρ(t * ) • Φ -1 t * ,ρ(t * ) = ϕ (Φ ′ t * ,ρ(t * ) ) 2 • Φ -1 t * ,ρ(t * ) ∂ t Φ ′ t,ρ(t) | t=t * • Φ -1 t * ,ρ(t * ) * * Recall that C 2 (S 1 , R) is compactly embedded in C 1 (S 1 , R)
, and thus given a bounded ball B ′ ⊂ C 2 (S 1 , R) and any δ > 0 there are

F δ := {ϕi} N δ i=1 ⊂ C 1 (S 1 , R) such that the union of the δ/2 C 1 -balls centered at the functions in F δ covers B ′ . Fix a bounded ball B ⊂ C 2 ([-ε2, ε2], C 2 (S 1 , R)), there is B ′ ⊂ C 2 (S 1 , R) such that γ([-ε2, ε2]) ⊂ B ′ for any γ ∈ B. Now fix δ > 0 and pick points {ti} M i=1 partitioning [-ε2
, ε2] such that γ(s1)-γ(s2) C 0 < δ/2 whenever s1, s2 ∈ [ti, ti+1] for some i. Then, for any (i1, ..., iM ) ∈ {1, ..., N } M fix a C 1 function γ (i 1 ,...,i M ) : [-ε2, ε2] → C 1 (S 1 , R) such that γ (i 1 ,...,i M ) (tj) = ϕi j . It is then easy to check that for every γ ∈ B there is (i1, ..., iM ) ∈ {1, ..., N } M such that γ is at C 1 distance no more than δ from γ (i 1 ,...,i M ) .

implying that

∂ t Q t,ρ(t) | t= t(ϕ) = - ϕ (Φ ′ t,ρ( t) ) 2 ∂ t Φ ′ t,ρ(t) | t= t • Φ -1 t,ρ( t) - ϕ ′ (Φ ′ t,ρ( t) ) 2 - ϕ • Φ ′′ t,ρ( t) (Φ ′ t,ρ( t) ) 3 • ∂ t Φ t,ρ(t) | t= t • Φ -1 t,ρ( t) = - ϕ Φ ′ t,ρ( t) (∂ t Φ t,ρ(t) ) t= t • Φ -1 t,ρ( t) ′
and by ( 33)

∂ t P t | t= t(ϕ) = -P ϕ Φ ′ t,ρ( t) (∂ t Φ t,ρ(t) ) t= t • Φ -1 t,ρ( t) ′ = -P ϕ Φ ′ t,ρ( t) (A ρ( t) + tA ∂tρ(t)| t= t ) • Φ -1 t,ρ( t) ′ .
Notice that this lemma implies that

∂ t P t | t= tρ( t) = -P [K t(ρ( t)) + tK t(∂ t ρ| t= t)],
where K t is given by [START_REF] Bálint | Synchronization versus stability of the invariant distribution for a class of globally coupled maps[END_REF].

Proof of Proposition 1.4. Recall that P t is the transfer operator of f • Φ t,ρ(t) which is an at least C 3 transformation of S 1 . Furthermore, there is K > 0 such that |Φ ′ t,ρ(t) | > 1 -Kt and thus f • Φ t,ρ(t) is uniformly expanding for any t with |t| < ε for some sufficiently small ε. Standard results on uniformly expanding maps (see [START_REF] Boyarsky | Laws of chaos: invariant measures and dynamical systems in one dimension[END_REF]) imply that P t has a spectral gap in C 2 (S 1 , R). In particular, 1 is an isolated and simple eigenvalue and Rρ(t) is the corresponding eigenspace. Because of the uniform Lasota-Yorke type inequalities ( 24)-( 25), a uniform bound can be given on the spectral gap of all P t such that t is small enough. This implies that we can find a positively oriented curve γ on the complex plane around 1 such that no other element of the spectrum of any P t is contained inside of it, and the projection formula

ρ(t) = 1 2πi γ (z -P t ) -1 ϕ(z)dz (34)
holds for any ϕ ∈ C 2 such that S 1 ϕ = 1. One can prove this using the decomposition of L t,ρ(t) as stated in [START_REF] Boyarsky | Laws of chaos: invariant measures and dynamical systems in one dimension[END_REF]Theorem 7.1.1] (exploiting the fact that F t,ρ(t) is C 2 and uniformly expanding which properties by [START_REF] Krzyżewski | On invariant measures for expanding differentiable mappings[END_REF] imply that F t,ρ(t) mixing) and a residue computation. For z ∈ γ we have (z -P t) -1 -(z -P t ) -1 = (z -P t) -1 (P t -P t )(z -P t ) -1 [START_REF] Wormell | On the validity of linear response theory in high-dimensional deterministic dynamical systems[END_REF] where we view (z -P t ) -1 as an operator acting on C 2 (S 1 , R), (P t -P t ) as an operator from C 2 (S 1 , R) to C 1 (S 1 , R) and (z -P t) -1 as acting on C 1 (S

where the last step is a residue computation. Thus we can write (1 -P t + tP K t)∂ t ρ(t)| t= t = -P K t(ρ( t)).

It is a well-known fact that invertible operators form an open set in the space of bounded linear operators between normed spaces. More precisely, if an operator T is invertible, then any operator S such that T -S ≤ T -1 -1 is also invertible. Since 1-P t is invertible as an operator from C 2 (S 1 , R) to C 1 (S 1 , R), the invertibility of 1 -P t + tP K t follows if we are able to choose t so small that tP K t C 2 →C 1 ≤ (1 -P t) -1 -1

C 2 →C 1 ⇔ t ≤ 1 P K t C 2 →C 1 (1 -P t) -1 C 2 →C 1 .
For this it suffices to show that t → P K t C 2 →C 1 (1 -P t ) -1 C 2 →C 1 is uniformly bounded in a small interval around zero. The norm (1 -P t ) -1 C 2 →C 1 is finite for all sufficiently small t (because of the above mentioned spectral gap), so it can be uniformly bounded in an interval [-ε 2 , ε 2 ]. As for K t , we can compute that sup Thus we can choose ε 3 further decreasing ε 2 such that 1 -P t + tP K t is invertible for t ∈ [-ε 3 , ε 3 ] and ∂ t ρ(t)| t= t = -(1 -P t + tP K t) -1 P K t(ρ( t)).

x∈S 1 |K t ϕ(x)| = = sup x∈S 1 ρ(t) Φ ′ t,ρ(t) A ϕ • Φ -1 t,ρ(t) (x) ′ = sup x∈S 1 ρ ′ (t)A ϕ Φ ′ t,ρ(t) + ρ(t)A ′ ϕ Φ ′ t,ρ(t) -ρ(t)A ϕ Φ ′′ t,ρ(t) (Φ ′ t,ρ(t) ) 2 • Φ -1 t,ρ(t) (x) • (Φ -1 t,ρ(t) ) ′ (x) ≤ 1 (1 -tK) 3 K(1 + tK) sup

x∈S 1 |ρ

 1 ′ (t)| + (K(1 + tK) + tK 2 ) sup x∈S 1 |ρ(t)|and a similar formula for sup x∈S 1 |(K t ϕ) ′ (x)| depending smoothly on ρ, ρ ′ , ρ ′′ and t.Since (t, x) → ρ(t)(x) is C 2 , K t C 2 →C 1 is uniformly bounded for t ∈ [-ε 2 , ε 2 ].

  d dt γ(t) C 3 < K 1 for all t values sufficiently small, d dt (Lγ)(t) C 3 ≤ K 1 . Now, assuming that d dt γ C 3 < K 1 , we carry out similar computations in order to bound d 2 dt 2 (Lγ)(t) C 2 .

	d 2 dt 2 (Lγ)(t) =	d dt	d ds	L s	s=t	γ(t) +	d dt	Ξ(t)

  Letting t → t in[START_REF] Wormell | On the validity of linear response theory in high-dimensional deterministic dynamical systems[END_REF] by Lemma 4.1 we have that∂ t (z -P t ) -1 | t= t = (z -P t) -1 ∂ t P t | t= t(z -P t) -1 .(36)Substituting ϕ = ρ( t) in[START_REF] Viana | Stochastic dynamics deterministic systems[END_REF] and differentiating both sides we get∂ t ρ(t)| t= t = 1 2πi γ (z -P t) -1 ∂ t P t | t=0 (z -P t) -1 ρ( t)(z)dz P t) -1 P (K t(ρ( t)) + tK t(∂ t ρ(t)| t= t))(z) 1 z -1 dz = -(1 -P t) -1 P (K t(ρ( t)) + tK t(∂ t ρ(t)| t= t)),

	=	1 2πi γ	(z -P t) -1 ∂ t P t | t= tρ( t)(z)	1 z -1	dz
	= -	1 2πi γ	(z -

1

, R).

§ If Lt were a linear operator, this would be enough to conclude that Lt is a contraction with respect to the Hilbert metric. Since it is nonlinear, we need to show it directly by estimating the Hilbert distance.

(corresponding author). The research of F. M. Sélley was supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 787304).

A Appendix

In this appendix we recall some simple formulae on derivatives that have been used extensively throughout the document. First of all higher derivatives of a product: Assume Φ, Ψ ∈ C 4 , then

which implies that there is a constant M > 0 such that

Higher derivatives of a composition: Assume Φ, Ψ ∈ C 4 , then

which implies that there is a numerical constant M such that

We also give the formulae for the higher derivatives of the inverse function. Let Φ ∈ C 4 be invertible and such that min

and so there exists an M > 0 such that


Fanni M Sélley
email: selley@lpsm.paris

A SELF-CONSISTENT DYNAMICAL SYSTEM WITH MULTIPLE ABSOLUTELY CONTINUOUS INVARIANT MEASURES

Keywords: 2020 Mathematics Subject Classication. Primary: 37A05, 37A10, 37E05; Secondary: 37M25, 65P99 Self-consistent dynamics, phase transition, β-map, absolutely continuous invariant measure, discrete transfer operator

A self-consistent dynamical system with multiple absolutely continuous invariant measures

Fanni M Sélley

self-consistent in the sense that the discrete time dynamics is dierent in each step depending on current statistics. The general framework admits popular examples such as coupled map systems. Motivated by an example of [START_REF] Blank | Ergodic averaging with and without invariant measures[END_REF], we concentrate on a special case where the dynamics in each step is a β-map with some β ≥ 2. Included in the denition of β is a parameter ε > 0 controlling the strength of self-consistency. We show such a self-consistent system which has a unique absolutely continuous invariant measure (acim) for ε = 0, but at least two for any ε > 0. With a slight modication, we transform this system into one which produces a phase transition-like behavior: it has a unique acim for 0 < ε < ε * , and multiple for suciently large values of ε. We discuss the stability of the invariant measures by the help of computer simulations employing the numerical representation of the self-consistent transfer operator. [START_REF] Arnold | Random dynamical systems[END_REF]. Introduction. A self-consistent dynamical system is a discrete time dynamical system where the dynamics is not the same map in every time step, but computed by the same rule from some momentary statistical property of the system. Such systems arise in problems of both physical and mathematical motivation, but their rigorous mathematical treatment so far has been restricted to some special cases, mainly coupled map systems.

Self-consistent systems bear resemblance to the larger framework of systems that are governed by laws that vary over time. The uniqueness and stability of the invariant measure is thoroughly studied for examples including non-autonomous dynamical systems [START_REF] Boyarsky | Absolutely continuous invariant measures for nonautonomous dynamical systems[END_REF][START_REF] Ott | Young Memory loss for time-dependent dynamical systems[END_REF], random dynamical systems [START_REF] Arnold | Random dynamical systems[END_REF][START_REF] Baladi | Random correlations for small perturbations of expanding maps[END_REF][START_REF] Bogenschütz | Stochastic stability of invariant subspaces[END_REF][START_REF] Buzzi | Absolutely continuous SRB measures for random LasotaYorke maps[END_REF] and random perturbations of dynamical systems [START_REF] Blank | Random perturbations of chaotic dynamical systems: stability of the spectrum[END_REF][START_REF] Baladi | On the spectra of randomly perturbed expanding maps[END_REF]. However, a self-consistent system is not a special case of any of these examples, as the dynamics in each step is not chosen via an abstract rule or drawn randomly from a set of possibilities, but is computed in a deterministic way from the trajectory of an initial probability distribution on the phase space.

The introduction of self-consistent systems dates back to [START_REF] Kaneko | Globally coupled chaos violates the law of large numbers but not the central limit theorem[END_REF], who studied globally coupled interval maps. In a globally (or mean-eld) coupled map system the F. M. SÉLLEY dynamics is the composition of the individual dynamics of a single site and a coupling dynamics which is typically the identity perturbed by the mean-eld generated by the sites, hence self-consistency arises from coupling. The eect of the meaneld is usually multiplied by a nonnegative constant ε called the coupling strength, which controls to what extent the self-consistency distorts the uncoupled dynamics. The literature studying coupled map systems is quite extensive. As systems of coupled maps are just loosely connected to the present work, we refrain from giving a complete bibliography, as a starting point see [START_REF] Chazottes | Dynamics of coupled map lattices and of related spatially extended systems[END_REF][START_REF] Sélley | Asymptotic properties of mean eld coupled maps[END_REF] and the references therein. Typically the existence and uniqueness of the invariant measure is studied in terms of the coupling strength. Most available results prove the uniqueness of the SRB measure for small coupling strength [START_REF] Blank | Self-consistent mappings and systems of interacting particles[END_REF][START_REF] Jiang | Equilibrium Measures for Coupled Map Lattices: Existence, Uniqueness and Finite-Dimensional Approximations[END_REF][START_REF] Liverani | A spectral gap for a one-dimensional lattice of coupled piecewise expanding interval maps[END_REF], but in some specic models phase transition-like phenomena can also be observed [START_REF] Bardet | Stochastically stable globally coupled maps with bistable thermodynamic limit[END_REF]: unique invariant measure for small coupling strength, and multiple for stronger coupling.

The literature of self-consistent systems not arising from coupled map systems is particularly sparse (in fact the only example known to us is the one discussed below). In this paper our goal is to study such a system which is in some sense much simpler than a coupled map system, hence interesting phase transition-like phenomena can be shown by less involved methods than the ones used for example in [START_REF] Bardet | Stochastically stable globally coupled maps with bistable thermodynamic limit[END_REF]. As results of this type are particularly hard to obtain in the coupled map setting, our results, although obtained in a simplied self-consistent system, contribute to the few existing examples.

Our main point of reference is Section 5 of [START_REF] Blank | Ergodic averaging with and without invariant measures[END_REF], specically the two systems dened by Example 5.2 which we now recall. Let X = [0, 1] and E µ = 1 0 x dµ(x), where µ is a probability measure on X.

Let (a) T µ (x) = x • E µ , (b) T µ (x) =
x/E µ mod 1 (where 1/0 mod 1 is dened as 0).

The map T µ induces an action on the space of probability measures, and an invariant measure of such a system is a probability measure for which µ = (T µ ) * µ. As Blank noted, in case (a) the only invariant measures are the point masses supported on 0 and 1, as T µ is a contracting linear map in all nontrivial case. Case (b) is more interesting since now T µ is a particular piecewise expanding map, a beta map, rst studied by [START_REF] Parry | On the β-expansions of real numbers[END_REF][START_REF] Parry | Representations for real numbers[END_REF][START_REF] Rényi | Representations for real numbers and their ergodic properties[END_REF]. Blank pointed out, that the self-consistent system has innitely many mutually singular invariant measures, including the Lebesgue measure. We are going to show that this picture is not complete, as the existence of multiple Lebesgue-absolutely continuous invariant measures (acims) can be shown.

The stability of these invariant measures is a more delicate question. By stability we mean that the invariant density attracts all elements of some neighborhood in a suitable norm, hence these equilibrium states rightfully describe an asymptotic behavior of the system. Rigorous results in this direction are only available in case of smooth self-consistent dynamics [START_REF] Keller | An ergodic theoretic approach to mean eld coupled maps[END_REF][START_REF] Bálint | Synchronization versus stability of the invariant distribution for a class of globally coupled maps[END_REF] and the treatment of piecewise smooth dynamics (such as example (b) of Blank) would require a completely different approach. However, to obtain a rough picture of the phenomena to be expected, computer simulations can be very useful. Numerical approximations of transfer operators and invariant densities have been extensively studied in the last few decades, typically by the help of generalized Galerkin-type methods. The idea behind these discretization schemes is the construction of a sequence of nite rank operators approximating the transfer operator of the dynamical system. The most notable scheme is Ulam's method [33], a relatively crude but robust method. The convergence of the xed points of the nite rank operators to the invariant density was rst proved by [START_REF] Li | Finite approximation for the FrobeniusPerron operator. A solution to Ulam's conjecture[END_REF] in case of one-dimensional dynamics, and since then, many generalizations have followed. For a comprehensive study see [START_REF] Murray | Discrete approximation of invariant densities[END_REF] and the references within. Better approximation can be achieved by higher order Galerkintype methods [START_REF] Ding | High order approximation of the FrobeniusPerron operator[END_REF][START_REF] Ding | Finite approximations of FrobeniusPerron operators. A solution of Ulam's conjecture to multi-dimensional transformations[END_REF]. For a more extensive survey of the discretization of the Perron-Frobenius operator see [START_REF] Klus | On the numerical approximation of the PerronFrobenius and Koopman operator[END_REF]. The approximation of invariant densities of non-autonomous dynamical systems on the other hand have a much more limited literature, focusing mainly on the setting of random dynamical systems [START_REF] Froyland | Ulam's method for random interval maps[END_REF][START_REF] Froyland | Quas Stability and approximation of random invariant densities for LasotaYorke map cocycles[END_REF][START_REF] Froyland | Quenched stochastic stability for eventually expanding-on-average random interval map cocycles[END_REF]. The self-consistent case, to the best of our knowledge is an uncharted territory. To make the rst steps, we consider specic systems (motivated by example (b) of Blank) with piecewise linear dynamics. The advantage of such systems is that the transfer operator maps the space of piecewise constant functions to itself, hence no discretization scheme is needed to compute pushforward densities. However, the task is not completely trivial, as the chaotic nature of the dynamics causes computational errors to blow up rapidly.

The setting and our results are summarized in Section 2. In Section 3 we introduce an auxiliary function ψ ε providing the main tool for the proofs of our results. In Section 4 we study a self-consistent system which interpolates linearly between the doubling map and case (b) of Blank's example by a parameter ε. We show that the system has a unique acim only in the case of ε = 0 (giving the doubling map) and has multiple absolutely continuous invariant measures for any ε > 0 (in particular for ε = 1, giving Blank's example.) In Section 5 we study a modied version of this self-consistent system which indeed exhibits a phase transition like-behavior: it has a unique acim if ε is smaller than some ε * > 0 and multiple acims if ε is suciently large. In Section 6 we showcase some results of computer simulations intended to study the stability of invariant densities with respect to the iteration of the self-consistent transfer operator.

2. The results. Let X = [0, 1] and denote the space of probability measures on X by M(X). For a measure µ ∈ M(X), let

E µ = 1 0 x dµ(x). (1)
Given an initial probability measure µ ∈ M(X) and ε ≥ 0, dene the dynamics T ε µ : X → [0, 1) as

T ε µ (x) = 2 + εF 1 E µ -2 x mod 1, x ∈ X, (2)
where

F ∈ C 1 (R, R) is such F (0) = 0.
The parameter ε controls to what extent the measure µ inuences the dynamics. (If µ 0 = δ 0 the Dirac mass concentrated on zero, we dene T ε δ0 ≡ 0.) In particular for ε = 0, µ has no inuence at all and

T 0 µ (x) = 2x mod 1, x ∈ X (3)
for any µ ∈ M(X).

We are going to study the self-consistent system

(M(X), (T ε • ) * •). (4)
An invariant measure of the system (4) is a measure µ ∈ M(X) such that

(T ε µ ) * µ = µ.
It is easy to see that the system (4) has many invariant measures: for instance the Lebesgue measure λ, since

E λ = 1 0 x dx = 1 2 ⇒ T ε λ (x) = 2x mod 1 implying that (T ε λ ) * λ = λ. (5)
In addition to this, innitely many mutually singular invariant measures exist: consider a measure uniformly distributed on a periodic orbit of the doubling map which is symmetric about We can see that this orbit is symmetric about 1 2 : if y is in this orbit, then so is 1-y. This implies that if µ = 1 k+1 k j=0 δ 2 j x mod 1 , we have E µ = 1 2 and T ε µ (x) = 2x mod 1, by which (T ε µ ) * µ = µ. Our main question is if (4) has multiple invariant measures absolutely continuous with respect to the Lebesgue measure (acims).

For ε = 0 we have seen that irrespective of the measure µ, the dynamics is always the doubling map. So in this case we have a unique absolutely continuous invariant measure.

We rst show that by taking the identity as F (producing Blank's example for ε = 1) this property is immediately lost as we introduce self-consistency.

Theorem 2.1. Consider the self-consistent system [START_REF] Bálint | Synchronization versus stability of the invariant distribution for a class of globally coupled maps[END_REF] and suppose that F (x) = x. Then for any ε > 0, at least two acims exist: one is Lebesgue, and the other is equivalent to Lebesgue.

We then show that under some additional assumptions on F , the uniqueness of the acim persists for ε small enough. But not indenitely: we also show that for suciently strong self-consistency this is not the case, i.e. multiple acims exist if ε is large enough. Theorem 2.2. Consider the self-consistent system [START_REF] Bálint | Synchronization versus stability of the invariant distribution for a class of globally coupled maps[END_REF] and suppose that F (x) > 0 for all x = 0 and F

(x) = O 1 | log x| as x → 0+.
1. There exists an ε * 1 > 0 such that for 0 ≤ ε < ε * 1 the only acim is the Lebesgue measure. 2. There exists an ε * 2 ≥ ε * 1 such that for ε ≥ ε * 2 at least two acims exist: one is Lebesgue, and the other is equivalent to Lebesgue.

An example of the function F for which Theorem 2.2 holds is F (x) = x 2k , where k ≥ 1 is an integer.

To discuss the stability of these invariant densities we performed a series of computer simulations. Based on the results of these computations we make two conjectures.

Conjecture 1. In the setting of Theorem (2.1), [START_REF] Arnold | Random dynamical systems[END_REF]. In case of ε = 0, the uniform density is stable with respect to the BV -norm

• BV = var(•) + | • |. 2.
In case of ε > 0, the uniform density is not stable with respect to the L 1 -norm

• L 1 = | • |
but there exists a dierent invariant density in BV which is stable with respect to the L 1 -norm.

Part (1) of this conjecture is the stability of Lebesgue measure under the doubling map, which is well known, we included it in the conjecture to highlight the bifurcation phenomenon.

Conjecture 2. In the setting of Theorem (2.2), the uniform density is stable with respect to the BV -norm for all ε > 0. Other acims are unstable with respect to the BV -norm.

3. The auxiliary function ψ ε . The aim of this section is to give the denition and discuss some properties of an auxiliary function on which the proofs of Theorems (2.1) and (2.2) rely. Let

T β (x) = βx mod 1, x ∈ [0, 1] (6)
such that β > 1, and denote by µ β the unique acim of the system. By the classical results of [START_REF] Rényi | Representations for real numbers and their ergodic properties[END_REF] we in fact know that such a measure exists, and it is equivalent to the Lebesgue measure. Remember that by our notation (1)

E µ β = 1 0 x dµ β (x). Dene ψ ε : (1, ∞) → R as ψ ε (β) = 2 + εF 1 E µ β -2 . (7)
Suppose there exists a β such that ψ ε ( β) = β. Notice that in this case µ β is an invariant measure of the self-consistent system (4). Indeed,

β = 2 + εF 1 E µ β -2 ⇒ T µ β (x) = βx mod 1,
and this implies that

(T µ β ) * µ β = µ β .
This shows that every xed point of ψ ε gives rise to an absolutely continuous invariant measure of (4). Moreover, if F ≥ 0, the absolutely continuous invariant measures of (4) and the xed points of ψ ε are in onetoone correspondence. So it suces to study number of xed points of the function ψ ε to prove our Theorems 2.1 and 2.2. We now state a lemma implying regularity properties of ψ ε .

F. M. SÉLLEY Lemma 3.1. 1. Let 2 ≤ β < 3.
There exists a constant K > 0 such that

|E µ β -E µ2 | ≤ K|β -2|(1 + | log |β -2||).
2. Let 2 < a 0 ≤ a 1 . There exists a constant K(a 0 , a 1 ) > 0 such that for all

β, β ∈ [a 0 , a 1 ], |E µ β -E µ β | ≤ K(a 0 , a 1 )|β -β |(1 + | log |β -β ||).
Proof. As proved in [START_REF] Parry | Representations for real numbers[END_REF], the unnormalized invariant density of T β can be given by the formula

h β (x) = ∞ n=0 1 β n 1 [0,T n β (1)) (x). (8)
Thus

1 0 h β (x) dx = ∞ n=0 1 0 1 β n 1 [0,T n β (1)) (x) dx = ∞ n=0 1 β n T n β (1) (9)
1 0 xh β (x) dx = ∞ n=0 1 0 1 β n x1 [0,T n β (1)) (x) dx = ∞ n=0 1 2β n (T n β (1)) 2 . ( 10
)
Notice in particular that since T 2 has range [0, 1), we have

T 0 2 (1) = 1 and T n 2 (1) = 0 for all n ≥ 1, hence h 2 (x) = 1. This implies that 1 0 h 2 (x) dx = 1 and 1 0 xh 2 (x) dx = 1 2 .
By denition and the above observation, we have

E µ β = 1 1 0 h β (x) dx 1 0 xh β (x) dx and E µ2 = 1 2 .
We can deduce the following uniform bounds on E µ β :

1 4 ≤ 1 2 • β -1 β ≤ E µ β ≤ 1 2 , (11)
where the last inequality is always strict if β is not an integer. Indeed,

1 2 ∞ n=0 1 β n (T n β (1)) 2 ∞ n=0 1
β n T n β (1) ≤ 1 2 , since (T n β (1)) 2 ≤ T n β (1
) for all n. In particular, if β is an integer, T 0 β (1) = 1 and T n β (1) = 0 for all n ≥ 1, hence the equality (and this is the only case that equality can occur). We get the lower bound from the fact 1 ≤ h β ≤ β β-1 (including just the rst term in the sum (8) versus all terms).

In particular observe that

E µ β ≤ E µ2 (12) so |E µ β -E µ2 | = E µ2 -E µ β .
Let us write h β (x) = 1 + g β (x) where

g β (x) = ∞ n=1 1 β n 1 [0,T n β (1)) (x).
Then

E µ β = 1 0 x(1 + g β (x))dx 1 0 1 + g β (x)dx = 1 2 + 1 0 xg β (x)dx 1 + 1 0 g β (x)dx , so E µ β ≥ 1 2 1 + 1 0 g β (x)dx ≥ 1 2 1 - 1 0 g β (x)dx ⇒ E µ2 -E µ β ≤ 1 2 1 0 g β (x)dx, so it is enough to show 1 0 g β (x)dx = ∞ n=1 1 β n T n β (1) ≤ N η + 1 β N ( 13
)
for η = β -2 and N = log η log 1/β to prove the claim of part (a). To show [START_REF] Chazottes | Dynamics of coupled map lattices and of related spatially extended systems[END_REF], we rst note that

∞ k=N +1 1 β k = 1 β -1 1 β N (14) Now notice that if T n-1 β (1) = β n-2 (β -2) < 1
β , the images T n β (1) fall under the rst branch. By the choice of N , this is exactly the case for all k = 1, . . . , N . So

1 β n T n β (1) = β n-1 (β -2) β n ≤ 1 β (β -2) for 1 ≤ n ≤ N implying N k=1 1 β k ≤ N β (β -2) (15)
Putting together ( 14) and ( 15) we get [START_REF] Chazottes | Dynamics of coupled map lattices and of related spatially extended systems[END_REF]. Part 2 can be proved in an analogous way to [23, Proposition 2] as a consequence of [START_REF] Keller | Continuity properties of transport coecients in simple maps[END_REF]Corollary 1]. This lemma has the following important corollary: Corollary 1. 1. For 2 ≤ β < 3, there exists a constant K > 0 such that

1 E µ β -2 ≤ K|β -2|(1 + | log |β -2||). 2. ψ ε (β) is continuous for all β ≥ 2.
Indeed, by part (1) of Lemma 3.1 and ( 11)

1 E µ β - 1 E µ2 = E µ2 -E µ β E µ β E µ2 ≤ 8K|β -2|(1 + | log |β -2||).
We note that part (2) of Lemma 3.1 serves only the purpose to draw the conclusion of part (2) of the above corollary. We now outline our main idea behind the proofs of Theorems 2.1 and 2.2. First observe that ψ ε (2) = 2, since Lebesgue measure is an invariant measure of the doubling map. But since Lebesgue is invariant for any β-map where β is an integer, we have

ψ ε (k) = 2 for all k ≥ 2, k ∈ N. So if for some β ∈ (k, k + 1) we have ψ ε ( β) > β, we can conclude that ψ ε has a xed point β * ∈ (k, k + 1)
. This implies that µ β * is an invariant measure of the self-consistent system (4) that is equivalent, but not equal to Lebesgue.

In the following sections we are going to prove Theorem 2.2 part (1) by showing that no such β exists for the stated values of ε, while we prove Theorem 2.1 and Theorem 2.2 part (2) by showing the existence of a β such that ψ ε ( β) > β in the settings considered. F. M. SÉLLEY 4. Proof of Theorem 2.1. In this section we consider the special case when F is the identity. Now [START_REF] Blank | Collective phenomena in lattices of weakly interacting maps[END_REF] takes the form

ψ ε (β) = 2 + ε 1 E µ β -2 .
We are going to show that for any ε > 0 the map ψ ε has another xed point in addition to β = 2. This implies Theorem 2.1, as discussed in Section 3. According to our arguments in Section 3, it is more than enough to prove the following proposition:

Proposition 1. Let ε > 0. For any δ > 0 there exists a β > 2, |β -2| < δ such that

ψ ε (β) > β. ( 16
)
This is the consequence of the following lemma, which claims that the log-Lipschitz continuity of β → E µ β at β = 2 stated in Lemma 3.1 cannot be improved to Lipschitz continuity: Lemma 4.1. There exists a sequence

β k → 2+ such that |E µ β k -E µ2 | > c β k |β k -2| for some c β k > 0 such that lim β k →2+ c β k = ∞.
To see that Lemma 4.1 readily implies [START_REF] Froyland | Finite approximation of Sinai-Bowen-Ruelle measures for Anosov systems in two dimensions[END_REF], note rst that we can discard the absolute values from the statement of this lemma by [START_REF] Boyarsky | Absolutely continuous invariant measures for nonautonomous dynamical systems[END_REF]. As we previously showed that E µ β is bounded away from 0, we also have

1 E µ β k - 1 E µ2 > cβ k (β k -2)
for some

cβ k > 0 such that lim β k →2+ cβ k = ∞. But then for cβ k > 1 ε we have ψ ε (β k ) = 2 + ε 1 E µ β k -2 > 2 + εc β k (β k -2) > β k .
By choosing k 0 so large such that cβ k 0 > 1 ε holds for all k ≥ k 0 , we obtain that

ψ ε (β k ) > k for all k ≥ k 0 .
The proof of Lemma 4.1. We are going to construct the sequence β k explicitly. Let β k be such that the rst k images of 1 fall under the rst branch of T β k and the k + 1-th image of 1 is 0. This means that

T β k (1) = β k -2 T 2 β k (1) = β k (β k -2) T 3 β k (1) = β 2 k (β k -2) . . . T k β k (1) = β k-1 k (β k -2) T k+1 β k (1) = β k k (β k -2) = 1 ≡ 0. (17)
To obtain β k , one simply has to nd the unique positive solution of [START_REF] Froyland | Ulam's method for random interval maps[END_REF]. It is easy to see that the thus dened β k → 2+ as k → ∞.

x T β3 (x)

β 3 -2 β 3 (β 3 -2)
β 2 3 (β 3 -2) = 1 β 3 = 1 β 3 3 (β 3 -2) 1 Figure 1. The choice of β k , k = 3 is pictured.
Straightforward calculations using ( 9) and [START_REF] Bogenschütz | Stochastic stability of invariant subspaces[END_REF] give that

1 0 h β k (x) dx = ∞ n=0 1 β n k T n β k (1) = 1 + k • 1 β k (β k -2) 1 0 xh β k (x) dx = ∞ n=0 1 2β n k (T n β k (1)) 2 = 1 2 1 + β k k -1 β k -1 • 1 β k (β k -2) 2
Now as

E β k = 1 1 0 h β k (x) dx 1 0 xh β k (x) dx,
and in particular E 2 = 1 2 , we obtain that

|E β k -E 2 | = 1 2 • 1 + β k k -1 β k -1 • 1 β k (β k -2) 2 1 + k • 1 β k (β k -2) - 1 2 = |β k -2| • β k k -1 β k -1 (β k -2) -k 2β k + 2k(β k -2) .
Since

β k k -1 β k -1 (β k -2)-k 2β k +2k(β k -2) → ∞ we can choose c β k = β k k -1 β k -1 (β k -2) -k 2β k + 2k(β k -2) .
5. Proof of Theorem 2.2. Throughout this section we assume that F > 0 for all x = 0 and F has the property

F (x) = O 1 | log x| as x → 0+.
5.1. Weak self-consistency: unique acim. In this section we prove the following proposition: Proposition 2. There exists an ε * 1 > 0 such that for 0 ≤ ε < ε * 1 the function ψ ε has a unique xed point.

As discussed, this implies part (1) of Theorem 2.2.

Proof. As we previously noted ψ ε (2) = 2. We are going to show that no other xed point exists provided that ε is small enough.

First notice that ψ ε ≥ 2 by ( 12), so no xed point exists which is smaller that 2.

We now show that despite the irregularity of β → E µ β in β = 2 (as stated in Lemma 4.1), ψ ε is Lipschitz continuous in β = 2 as a result of the derivative of F vanishing in 0 at an appropriate rate.

Lemma 5.1. There exists a δ = δ(F ) > 0 and an M = M (F ) > 0 such that

|ψ ε (β) -ψ ε (2)| ≤ εM |β -2| (18)
for all β ∈ (2, 2 + δ).

Proof. Let G : (0, h) → (0, ∞) (for some small h > 0) be the inverse function of the

map x → Kx(1 + | log x|) dened for small x > 0, where K is from Corollary 1 part (1), that is, 1 E µ β -2 ≤ K(β -2)(1 + | log(β -2)|)
for all β > 2 suciently close to 2. Taking arbitrary α > 1, since x

1 α > Kx(1 + | log x|)
for all suciently small x > 0, we have that x α < G(x) and so

1 α| log x| < 1 | log G(x)| = KG (x)
for all suciently small x > 0. Therefore, since by assumption

F (x) is O 1 | log x| as x → 0+ we have that F (x) < M G (x)
for some constant M = M (F ) and for suciently small x > 0. Since F (0) = 0, it follows that F (x) < M G(x) for suciently small x > 0 and so

ψ ε (β) -2 = εF 1 E µ β -2 ≤ M εG 1 E µ β -2 ≤ εM (β -2).
Now we show that for small enough ε, we have ψ ε (β) < β for 2 < β < 2 + δ. By [START_REF] Froyland | Quas Stability and approximation of random invariant densities for LasotaYorke map cocycles[END_REF],

ψ ε (β) -2 ≤ εM (β -2) ψ ε (β) ≤ εM (β -2) + 2 and εM (β -2) + 2 < β whenever ε < 1 M .
We now show

ψ ε (β) = 2 + εF 1 E µ β -2 < β for β ∈ [2 + δ, ∞). (19)
By [START_REF] Buzzi | Absolutely continuous SRB measures for random LasotaYorke maps[END_REF] we have

F 1 E µ β -2 ≤ F 2 β -1 . Since for β ≥ 2 + δ 2 + εF 1 E µ β -2 ≤ 2 + εF 2 1 + δ < 2 + δ ≤ β when ε < δ F ( 2 1+δ )
, [START_REF] Froyland | Quenched stochastic stability for eventually expanding-on-average random interval map cocycles[END_REF] is proved. This implies the statement of the proposition with

ε * 1 = min    1 M , δ F 2 1+δ    .

5.2.

Strong self-consistency: multiple acims. In this section we prove Theorem (2.2), part [START_REF] Baladi | On the spectra of randomly perturbed expanding maps[END_REF]. For this it suces to nd a single β such that ψ ε ( β) > β holds for suciently large ε.

We in fact show that we can achieve ψ ε ( β) > β for arbitrary β ∈ (k, k + 1), provided that ε is large enough in terms of β.

Proposition 3. Let β ∈ (k, k + 1). There exists an ε * 2 = ε * 2 ( β) such that ψ ε ( β) > β holds for all ε > ε * 2 ( β).
Proof. We would like to have

ψ ε ( β) = 2 + εF 1 E µ β -2 > β This holds if ε > β -2 F 1 Eµ β -2
, so by the choice of

ε * 2 = β -2 F 1 Eµ β -2
the proposition is proved.

This proposition has the following corollary:

Corollary 2. The self-consistent system (4) can have an arbitrarily large nite number of invariant measures equivalent to the Lebesgue measure, provided that ε is large enough.

To see this, let ≥ 2 be an arbitrary integer. Choose βk ∈ (k, k + 1), k = 2, . . . , . Let ε * = max k∈{2,..., } ε * 2 ( βk ). Then Proposition 3 implies that for ε > ε * , the function ψ ε has a xed point on each of the intervals (k, k +1), k = 2, . . . , implying a total number of at least acims. 6.1. Illustrations of β → ψ ε (β). To illustrate the results of the previous sections, we present some computer approximations of the curve β → ψ ε (β) for some appropriate functions F . Since T β (x) = βx mod 1 is ergodic (proved in [START_REF] Rényi | Representations for real numbers and their ergodic properties[END_REF]), we can approximate E µ β by computing the reciprocal of ergodic averages. This means we can approximate the graph of ψ ε by

(b) N = 10 7 , ∆ = 10 -5 , β1 = 2.03, ε = 0.2
ψ ε (β) ≈ P x,N (β) = 2 + εF N N -1 n=0 T n β (x) - 2
for (Lebesgue) almost every x ∈ [0, 1] and N large. We rst consider the setting of Theorem 2.1: the case when F is the identity. We illustrate on Figure 2 that no matter how small ε is, the curve approximating β → ψ ε (β) always grows above the line x = y for β values suciently close to 2. This shows that β → 1 Eµ β -2 cannot be Lipschitz at 2, otherwise multiplying it with suciently small ε would force the curve β → ψ ε (β) = 2 + ε 1 Eµ β -2 to always stay below x = y.

We now consider the setting of Theorem 2.2. We rst study the setting of part (1), that is when ε is suciently small. Now

β → F 1 Eµ β -2 is Lipschitz at 2 as a result of F (x) = O 1 | log x| and the log-Lipschitz continuity of β → 1 Eµ β -2.
So suciently small ε will not let the curve β → ψ ε (β) rise above x = y. The results of our computer simulations are pictured on Figure 3, showing clearly that the curve approximating β → ψ ε (β) has a single intersection with the diagonal x = y at β = 2.

The setting of Theorem 2.2 part (2) is studied on Figure 4 rst for the special case F (x) = x 2 . We can clearly see, as suggested by Corollary 2, that if ε is larger and larger, the curve approximating β → ψ ε (β) has intersections with the line x = y on more and more intervals between two consecutive integers.

Similar plots can be made for F (x) = x 4 and F (x) = x 6 . On Figure 5 we can see that for suciently large ε the curve approximating β → ψ ε (β) has intersections with the line x = y, indicating multiple invariant densities. 6.2. Stability of the invariant densities. Although the existence of a unique or multiple invariant measures is an interesting phenomenon on its own, it is natural to further inquire about their stability. Let dµ = f dλ and let P ε f denote the transfer operator associated to the dynamics T ε µ . This operator maps the density h of a measure ν to the density of the measure (T ε µ ) * ν. Consider an initial measure dµ 0 = f 0 dλ. The associated dynamics is T ε µ0 and the associated transfer operator is P ε f0 . Let f 1 = P ε f0 f 0 be the pushforward density and dµ 1 = f 1 dλ be the pushforward measure. Continuing this further we obtain the nth step pushforward density as

f n+1 = P ε fn f n , n = 0, 1, . . .
As an ease of notation, consider the self-consistent transfer operator F ε dened as

F ε (f n ) = P ε fn f n .
An invariant density of the self-consistent system is a xed point of this operator, and we can study its stability, that is, if densities suciently close to it in some metric converge to it. A function space well suited to this problem is the space BV of functions of bounded variation, and the metric is the one given by the bounded variation norm.

(a) ε = 2.5, β1 = 3
(b) ε = 400, F (x) = x 6
In case of ε = 0, producing the doubling map as the dynamics in each step, it is well known that the Lebesgue measure is stable in the sense that all measures with a density of bounded variation converge to it exponentially fast. It is a natural question to ask if this still holds in our self-consistent system when Lebesgue is the unique absolutely continuous invariant measure, so in the setting of Theorem 2.2 part (1). On the other hand, in the setting of Theorem 2.1 and Theorem 2.2 part (2) we have proved that the Lebesgue measure is a unique invariant measure of the system for suciently small values of ε, but for larger values multiple invariant measures exits. It would be interesting to see what kind of bifurcation occurs at the critical value of the coupling: is Lebesgue stable for small ε, and does it stay that way when multiple invariant measures arise, or does it lose its stability? Also, are the new invariant measures stable or unstable? For example, it would be interesting to show a similar behavior to the pitchfork bifurcation observed in the system of coupled fractional linear maps of [START_REF] Bardet | Stochastically stable globally coupled maps with bistable thermodynamic limit[END_REF]: they show that the stable, unique invariant measure loses stability if the coupling strength is suciently increased, and two new stable invariant measures arise.

To study these questions we present the results of some computer simulations. We rst note that there exists an explicit expression for the transfer operator associated to the β map T β (x) = βx mod 1:

P β f (x) = y∈T -1 β (x) f (y) |T β (y)| , f ∈ L 1 ,
or more explicitly one can write

P β f (x) =    1 β β k=0 f x+k β if 0 ≤ x < β -β β β -1 k=0 f x+k β if β -β ≤ x ≤ 1. ( 20
)
It is clear that if f is a a nite linear combination of indicators of intervals (a step function), P β f is also a step function. As functions of this kind are easy to store and manipulate by computer programs, we will restrict to working with densities of this kind. Note that this restriction means that we do not need to apply a discretization scheme to compute pushforward densities. Let the step function f be represented by the vectors x f ∈ [0, 1] N such that x f (1) = 0 and x f (N ) = 1 and y f = (R + 0 ) N -1 . The vector x f contains the jumps of the step function in increasing order and y f contains the respective heights of the steps. We dene the total variation of f as

var(f ) = N -2 i=1 |y f (i + 1) -y f (i)|.
Our initial densities will be generated in the following way: x f (2), . . . , x f (N -1) are N -2 numbers drawn from the uniform random distribution on (0, 1) and then ordered increasingly. The values ỹf (1), . . . , ỹf (N -1) are also drawn from the uniform random distribution on (0, 1). We dene

y f (i) = 1 N -1 j=1 ỹf (j)(x f (j + 1) -x f (j)) • ỹf (i), i = 1, . . . , N -1,
as this is an easy way to generate a fairly general step function of integral 1.

F. M. SÉLLEY

Let the `expectation' associated to the density f be dened as

E f = N -1 i=1 y f (i) 2 (x f (i + 1) 2 -x f (i) 2 ).
The action of the self-consistent transfer operator F ε (f ) can be computed in the following way: take an initial density represented by the pair of vectors (x f0 , y f0 ) ∈

[0, 1] N0 × (R + 0 ) N0-1 . Compute β 0 = 2 + εF 1 E f 0
-2 dening the dynamics T β0 . Compute the new vector of jumps x f1 as the vector containing the values T β0 (x f0 (i)), i = 1, . . . , N 0 and 1 in increasing order.

To compute the new heights of the steps, ρ > 0 such that ρ < min |x f1 (i + 1) -x f1 (i)| and compute y f1 (i) = P β0 f 0 (x f1 (i + 1) -ρ) for i = 1, . . . , N 0 by using the formula [START_REF] Jiang | Equilibrium Measures for Coupled Map Lattices: Existence, Uniqueness and Finite-Dimensional Approximations[END_REF]. As f 0 is piecewise constant and exactly stored, the evaluation of a f 0 at prescribed places is not an issue. Set N 1 = N 0 + 1 and repeat with (x f1 , y f1 ) ∈ [0, 1] N1 × (R + 0 ) N1-1 . Our procedure is the following: we generate K 1 × K 2 of the above described step functions in the following way: for each k = 1, . . . , K 1 we generate a random integer between 1 and M (denote this by m(k)) and generate K 2 step functions with m(k) inner jumps (this means not counting 0 and 1). This gives us a fairly general pool of initial densities.

We then compute the T long trajectories of the densities with respect to the self consistent transfer operator F ε . We are going to use the notation

f i t = F ε f i 0
where the lower index refers to time and the upper index to which one of the initial densities we are considering, so t = 1, . . . , T and i = 1, . . . , K 1 × K 2 . However, our experience is that computational errors grow rapidly and seriously skew our results in the long run. So in each iteration we normalize f i t by the numerical integral

I(f i t ) = N j=1 y f i t (j)(x f i t (j + 1) -x f i t (j)), dening
f i t = 1 I(f i t ) f i t .
This assures us that the Perron-Frobenius operator is indeed applied to a density function. This density f i t might not be the actual density F t ε (f i 0 ), but one close to it. It can be thought of as F t ε (g i 0 ) for some dierent initial density g i 0 by anticipating a type of shadowing property.

We are going to study two mean quantities of the densities. Dene the mean slope of the densities at time t as

β(t) = 1 K 1 × K 2 K1×K2 i=1 2 + εF 1 E f i t - 2
and the mean total variation as

var(t) = 1 K 1 × K 2 K1×K2 i=1 var( f i t ).
Finally, we dene our notions of stability of an invariant density f * (ε). We are going to study two types of stability.

1. Stability in the BV-norm:

var f * (ε) (t) = 1 K 1 × K 2 K1×K2 i=1 var( f i t -f * (ε)) → 0 as t → ∞ 2. Stability in the L 1 -norm: int f * (ε) (t) = 1 K 1 × K 2 K1×K2 i=1 [0,1] | f i t -f * (ε)|dλ → 0 as t → ∞
If an explicit expression for f * (ε) not available, then in practice f * (ε) is f T (ε) for some T considerably larger than T and some xed initial density, as we assume this is a good approximation of an invariant density.

Note that L 1 -stability is weaker than BV -stability. We further remark that var(t) → 0 as t → ∞ implies stability of the uniform invariant density in the BV -sense.

We rst consider the setting of Theorem 2.1. In this case we studied the values ε = 0.1 and ε = 0.2. From Figure 2 we can read that there exists an invariant density f * (ε) for which 2 + εF

1 E f * (0.1) -2 ≈ 2.0006 and 2 + εF 1 E f * (0.2) -2 ≈ 2.0181.
Running our simulations we can see from Table 1 that in both cases var does not converge to zero, so the uniform density is not stable in BV -sense. To convince ourselves more thoroughly that the constant density is indeed not stable, we made computations with pools of initial densities very close to the constant one in the sense that var(f i 0 ) < 10 -4 . In this case convergence is slower, but it is clearly not to the constant density, see Table 2.

On the other hand, we can also read form Table 1 that β converges to 2.0006 and 2.0181 for ε = 0.1 and 0.2 respectively, and this suggests that the nontrivial invariant densities are stable. Table 3 provides further evidence pointing to the stability of the nontrivial invariant densities in the L 1 -sense. Note that convergence of int f * (ε) to zero is not something to be expected since f T (ε) is just an approximation of f * (ε). However, as T is larger, int f * (ε) converges to smaller values which supports our hypothesis.

So it seems likely that the constant density looses its stability (in the L 1 -sense) as ε becomes larger than zero, and a new stable invariant density arises (in the L 1 -sense).

To back our assumption that var and β correctly describes the behavior of a typical density, we plot all the total variation and β for the densities on Figures 6 and7 respectively, for each time instance. We can see that they all converge to a single value, hence the averaging does not give an average of dierent asymptotic behaviors but shows us the true one.

The asymptotic densities obtained from iterating an appropriate initial density for both cases ε = 0.1 and ε = 0.2 are pictured on Figure 8.

Finally we discuss the convergence of our method depending on the value of ε.

By convergence we mean that the quantities var and β settle at a value var * and β * in the sense that |var * -var(t)| < 10 -4 and |β * -β(t)| < 10 -4 for t = T 0 , . . . , T where we choose T 0 , T 1 to be some xed large numbers. In Figure 9 we plotted the mean total variation and mean slope of the last 50 iterates of our density pool for a range of ε values, that is, we chose T 0 = 150 and T 1 = 200. If this produces a considerable range of values for a single ε (an interval having length larger then 10 -4 above an ε value), then the method does not converge. So we can read from this gure that our method is converges until ε ≈ 0.4. Now we move on to consider the setting of Theorem 2.2. We studied the cases F (x) = x 2 , F (x) = x 4 and F (x) = x 6 for a few values of ε for which simulations similar to the ones discussed for the case of F (x) = x clearly suggest unique or multiple absolutely continuous invariant measures. In the rst columns of Tables 4 we see computations in the cases when the constant density is the unique invariant one, see Figure 3. We can see that in all cases var(t) decreases rapidly, suggesting the stability of the invariant constant density in BV -sense.

In the second columns of Tables 4, we considered situations where multiple invariant densities exist, see Figure 4 (a) and Figure 5. In case of F (x) = x 2 , our method does not converge in the sense that the values of var do not settle at a value with precision 10 -4 if we consider any subinterval t = T 0 , . . . , 100. We did other experiments with T 1 = 500 and T 1 = 1000, T 0 = T 1 -100 but obtained similar results. However, a slow decrease of var is observable, so we do not exclude the possibility that we could obtain convergence for T 1 = 10 k for some k large, but our limited resources prohibit us from carrying out such computations in a reasonable amount of time.

On the other hand, for F (x) = x 4 and F (x) = x 6 we see fast convergence of var to zero implying that the constant density is likely to be stable one in BV -sense. On Figure 10 we discuss the convergence of our method for a range of ε similarly as we have done for F (x) = x (plotted on Figure 9). We can see that for F (x) = x 2 our method becomes erratic for larger ε values than ε ≈ 1.6. But our method seems to converge nicely in the F (x) = x 4 and F (x) = x 6 case for all values of ε considered.

7. Concluding remarks. To answer the question regarding stability rigorously, the careful study of the self-consistent transfer operator is necessary. When ε = 0, F n ε = P n T , where P T is the transfer operator of the doubling map. The stability result regarding the doubling map can be proved by elementary means: one can show by explicit calculations that P n T f -1 BV ≤ C 2 n f -1 BV for all f ∈ BV ([0, 1], R) such that 1 0 f = 1 However, when ε > 0 we have to deal with the self-consistency. In this case F n ε = P ε f0 P ε f1 . . . P ε fn-1

for some densities f 0 , . . . f n-1 , so the problem does not simplify to the study of a single linear operator. Provided that ε is small enough, it is natural to expect that P ε f0 P ε f1 . . . P ε fn-1 is `close' to P n T in some sense, hence acts similarly. In the coupled map systems of [START_REF] Keller | An ergodic theoretic approach to mean eld coupled maps[END_REF] and [START_REF] Bálint | Synchronization versus stability of the invariant distribution for a class of globally coupled maps[END_REF] giving rise to self-consistent dynamics this is precisely the strategy to prove stability in BV . However, the (C 3 or C 2 and Lipschitz second derivative) smoothness of the stepwise dynamics is an essential part of their proof. In the setting of this paper, the stepwise dynamics T ε µ is discontinuous, posing a major technical diculty, so it can also be the case that dierent tools are needed to study the asymptotic behavior of the operator F ε . In correspondence to the numerical stability of the nontrivial invariant densities of the case F (x) = x, we also believe that in full generality only stability in the L 1 -sense is to be expected.

Another question that arises observing the Figures 2 and4 is if the intersection of the numerical approximation of β → ψ ε (β) and the line x = y approximates a single intersection of β → ψ ε (β) and x = y or innitely many accumulating ones. As the regularity of ψ ε is quite low, one can quite possibly imagine innitely many intersections reminiscent of the innitely many accumulating zeros of the trajectory of Brownian motion. This would be interesting, as it would give innitely many absolutely continuous invariant measures.

Figure 2 .

 2 Figure 2. Approximation of ψ ε with ergodic averages, F (x) = x. P x,N (β) is plotted for β ∈ [2, β 1 ] with gridsize ∆ and for x drawn uniform randomly from [0, 1]. The line x = y is plotted.

2

 2

Figure 3 .

 3 Figure 3. Approximation of ψ ε with ergodic averages. P x,N (β) is plotted for β ∈ [2, 3] with gridsize ∆ = 10 -4 , N = 10 6 and for x drawn uniform randomly from [0, 1], ε = 0.8. The line x = y is plotted.

  ε = 200, β1 = 5

Figure 4 .

 4 Figure 4. Approximation of ψ ε with ergodic averages, F (x) = x 2 . P x,N (β) is plotted for β ∈ [2, β 1 ] with gridsize ∆ = 10 -4 , N = 10 6 and for x drawn uniform randomly from [0, 1]. The line x = y is plotted.

Figure 5 .

 5 Figure 5. Approximation of ψ ε with ergodic averages, P x,N (β) is plotted for β ∈ [2, 3] with gridsize ∆ = 10 -4 , N = 10 6 and for x drawn uniform randomly from [0, 1]. The line x = y is plotted.

F 2 Figure 6 .

 26 Figure 6. Total variation for each f i t as a function of time, F (x) = x, T = 100, K 1 × K 2 = 100.

Figure 7 . 2 Figure 8 .

 728 Figure 7. Associated slope to each f i t as a function of time, F (x) = x, T = 100, K 1 × K 2 = 100.

Table 4 .

 4 Computation of the mean total variation var. T = 100,M = K 1 = K 2 = 10. Left hand side: F (x) = x 2 , center: F (x) = x 4 , right hand side: F (x) = x 6 .

  (a) F (x) = x 2 , ∆ = 10 -3 , E = 2F (x) = x 4 , ∆ = 0.5, F (x) = x 6 , ∆ = 10 -2 , E = 12

Figure 10 .

 10 Figure 10. For each ε = k • ∆, k = 1, . . . , E/∆ the values var(t) are plotted for t = 150, . . . , 200, K 1 × K 2 = 100.

6 .

 6 Numerical results.

		2.001						2.03
		2.0009
		2.0008						2.025
		2.0007						2.02
		2.0006
	( )	2.0005					( )	2.015
		2.0004
		2.0003						2.01
		2.0002						2.005
		2.0001
		2						2
		2	2.0002	2.0004	2.0006	2.0008	2.001	2	2.005	2.01	2.015	2.02	2.025	2.03
	(a) N = 10 8 , ∆ = 10 -6 , β1 = 2.001,
	ε = 0.1

Table 1 .

 1 Computation of the mean total variation var and mean slope for F (x) = x. T = 100, M = K 1 = K 2 = 10. Left hand side: ε = 0.1, right hand side: ε = 0.2.

	1 ,

Table 2 .

 2 Computation of the mean total variation var and mean slope for F (x) = x, var(f i 0 ) < 10 -4 . T = 300, M = K 1 = K 2 = 10. Left hand side: ε = 0.1, right hand side: ε = 0.2.

	t	int f * (0.1) int f * (0.2)
	0	0.4096	0.5018
	10	0.0101	0.0321
	20	0.0049	0.0174
	30	0.0028	0.0097
	40	0.0017	0.0064
	50	0.0011	0.0040
	60	0.0007	0.0021
	70	0.0005	0.0007
	80	0.0004	0.0003
	90	0.0003	0.0002
	100 0.0002	0.0001
	110 0.0002	0.0001
	120 0.0002	0.0001
	130 0.0001	0.0001
	140 0.0001	0.0001
	150 0.0001	0.0001
	160 0.0001	0.0001
	170 0.0001	0.0001
	180 0.0001	0.0001
	190 0.0001	0.0001
	200 0.0001	0.0001

Table 3 .

 3 Computation of int f * (ε) and mean slope for

  Figure 9. For each ε = k • 10 -3 , k = 1, . . . , 10 3 the values var(t) and β(t) are plotted for t = 150, . . . , 200, F (x) = x, K 1 ×K 2 = 100.

		1											2.3
		0.95
													2.25
		0.9
		0.85											2.2
	var	0.8											2.15
		0.75
		0.7											2.1
		0.65
													2.05
		0.6
		0.55											2
		0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
	t		ε = 1	ε = 2.5			t	ε = 1	ε = 35			t		ε = 1		ε = 400
	0		6.4694 5.0341			0	5.3766 6.4806			0		7.7582 5.1444
	5		1.2483 1.1951			5	1.2370 1.1785			5		1.3852 1.0419
	10		0.9808 0.8849			10	0.9882 0.9788			10	0.5224 0.5990
	15		0.9914 0.8471			15	0.6884 0.9873			15	0.0210 0.1435
	20		0.9976 0.8447			20	0.0474 0.5994			20	0.0007 0.0026
	25		0.9994 0.8275			25	0.0015 0.1093			25	0.0000 0.0002
	30		0.9999 0.8418			30	0.0000 0.0347			30	0.0000 0.0000
	35		1.0000 0.8305			35	0.0000 0.0023			35	0.0000 0.0000
	40		1.0000 0.8468			40	0.0000 0.0001			40	0.0000 0.0000
	45		1.0000 0.8321			45	0.0000 0.0000			45	0.0000 0.0000
	50		0.9531 0.8358			50	0.0000 0.0000			50	0.0000 0.0000
	55		0.5357 0.8382			55	0.0000 0.0000			55	0.0000 0.0000
	60		0.1802 0.7651			60	0.0000 0.0000			60	0.0000 0.0000
	65		0.0128 0.6366			65	0.0000 0.0000			65	0			0.0000
	70		0.0004 0.5117			70	0.0000 0.0000			70	0			0
	75		0.0000 0.4535			75	0		0.0000			75	0			0
	80		0.0000 0.4269			80	0		0.0000			80	0			0
	85		0.0000 0.4255			85	0		0			85	0			0
	90		0.0000 0.4015			90	0		0			90	0			0
	95		0.0000 0.3966			95	0		0			95	0			0
	100 0.0000 0.3813			100 0		0			100 0			0

Acknowledgments. The research was supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 787304) and by the Hungarian National Foundation for Scientic Research (NKFIH OTKA) grant K123782. The author would like to express gratitude to Péter Bálint, Imre Péter Tóth and Péter Koltai for helpful discussions. The author also expresses gratitude to an anonymous referee for providing a substantial amount of help to correct the proof of Lemma 5.1.

The author was supported by ERC grant No 787304 and NKFIH OTKA grant K123782.