Didier Emmanuel Rousseau

Didier Felbacq

Ray Chaos in a Photonic Crystal

Keywords: 42, 65, Sf -First pacs description 78, 67, Pt -Second pacs description nn, mm, xx -Third pacs description

Ray chaos in a photonic crystal

Introduction and issue. -Photonic crystals, i.e. artificial periodic structures exhibiting a photonic band gaps [1], have been studied at very different angles since the pioneering works of S. John and E. Yablonovitch [2,3]. From the basic effect of acting like mirrors, the engineering of band structures and the recognition that they behave as metamaterials [4], they have shown a wide range of physical phenomena. However, they have been hardly considered in the high frequency domain: it is indeed at most, say, the first three bands which are considered for band structure engineering or effective properties. In this work, we consider the propagation of light in periodic dielectric structures when the wavelength is very small with respect to the scatterers size. In this limit, band diagrams become highly complicated because of a huge number of eigenmodes. This is the domain of wave chaos where puzzling phenomena such as random lasing [5] or localization [6,7] may be expected. In electromagnetism and optics, most of the studies dealt with chaotic cavities [5,[START_REF] Stöckmann | [END_REF][14][15][16][17][18][19] but only few of them with open chaotic systems such as optical fibers [20] or photonic crystals [21,22]. In the case of a two dimensional photonic crystal, chaos is revealed by the help of the statistical properties of the energy level spacing [21,22]. They are computed from the eigenfrequencies of the Helmholtz equation solved with proper boundary conditions. These approaches might indicate that light dynamics in a 2D photonics-crystal is "quasiintegrable" [22] meaning that the light dynamics exhibits regular and chaotic trajectories characterized through dif-ferent level-spacing distributions. However, the dynamical properties of light are not a direct output of such an approach. Here we show that the dynamics of light and the transition from regular to chaotic trajectories can be understood with the help of light rays. We work then in the geometric approximation, the light rays representing the trajectories of the energy [20]. In order to be valid, our model requires that the wavelength of light be at least ten times smaller than the cylinder radius (see section I of the Supplementary Materials for a detailed discussion [40]). Within such an approximation, we show that the transition from a regular motion to a chaotic one is imprinted in the dynamics of the transmitted rays, the dynamics of the reflecting rays being always chaotic. Within the rays approximation, we compute the Lyapunov exponent that quantify the exponential sensitivity to the initial conditions. From the knowledge of the dynamics, we derive the diffusive properties of light in the photonic crystal. We show a transition from a brownian diffusion to an "anomalous" regime.

The Photonic Billiard. -The photonic crystal under study is depicted in Fig. 1: it is a biperiodic set of dielectric cylinders, in which the trajectories of a ray are studied, by applying Snell-Descartes relations at the boundary of the cylinders. If the cylinders are perfectly reflecting instead of transparent then the system is nothing else than a periodic Lorentz gas. It has an extensive bibliography both in the physical and mathematical literatures (see e.g. [START_REF] Gaspard | Chaos, scattering and statistical mechanics[END_REF] and references hereby). It consists Fig. 1: a) A small portion of the photonic crystal. T is the period of the photonic crystal in both vertical and horizontal directions. R is the radius of the cylinders. b) Notations used in the text. "Ap" referred to the p th point defined as the intersection between the light ray and the cylinder surface at the p th step, "Bp" is the p th point at the cylinder surface after refraction inside the cylinder. α denotes the angle between the ray and the horizontal axis. θA p (resp. θB p ) is the angle between the normal at the point Ap (resp. Bp) and the horizontal axis.

θ 1 Ap (resp. θ 1
Bp ) is the angle between the normal at the point Ap (resp. Bp) and the incoming light ray (resp. outcoming light ray).

in an ensemble of noninteracting point particles moving freely with elastic collisions on fixed scatterers (the cylinders). The photonic billiard can be seen as the refractive extension of the periodic Lorentz gas. Upon this analogy a light ray represents a particle trajectory and, instead of having infinite walls, the repulsive potential is finite with a value given by the optical index of the cylinders. An important difference with the periodic Lorentz gas is that now rays carry energy that is split between the refracted and the reflected ray. For moderate values of the cylinder optical index, the reflection coefficient in energy (4% for one single glass/air interface) is small. In such a case, energy is predominantly carried through the rays that are transmitted through the cylinder. See the Supplementary Materials (I A) [40] for more details. Whereas there exists a huge literature about the periodic Lorentz gas with hardwall scatterers, we have not been able to find references on the "refractive periodic Lorentz gas" [START_REF]We use the terminology "photonic billiard" or "refractive Lorentz gas[END_REF]. The Lorentz gas is an unfolding of the Sinai Billiard [START_REF] Ya | [END_REF] for which the dynamics is always chaotic,whatever the photonic crystal parameters, making impossible long-range precision concerning a single trajectory. So, despite its deterministic character, this system can only be studied from a statistical point of view. Transport is diffusive and the square geometry is known [26,27] to exhibit superdiffusion processes (also called anomalous diffusion) i.e. diffusion for which the mean square displacement < r 2 > grows faster than linearly with time t like < r 2 >∝ t ln t [28]. In the past few years, superdiffusion gained an increasing interest in optics [31] with the observation of weak localization [32] in disordered samples designed to exhibit anomalous diffusion [33]. Our photonic billiard, or refractive Lorentz gas, shows a richer dynamics as compared to the Lorentz gas and the main result here is that the geometry studied in this paper also exhibits diffusion processes that can be tuned from nearly ballistic to brownian, including anomalous diffusion. The letter is organized as followed. First the dynamics of the refractive Lorentz gas is discussed through Poincaré surfaces of section. The structure exhibits soft chaos, i.e. regular paths (also called ballistic paths), as well as chaotic motion as a consequence of its quasi-integrable dynamics [START_REF] Ozorio De | Hamiltonian Systems Chaos and Quantization[END_REF]. Above a threshold value of the normalized period t = T /R where T is the crystal period and R the cylinder radius, the dynamics becomes completely chaotic. Chaos strength is quantified through the computation of Lyapunov exponents. Diffusive properties are computed, showing that the refractive Lorentz gas exhibits superdiffusion due to the existence of Levy flights. The rays dynamics. -The dynamical properties of rays in the crystal are now described. Numerical computations were performed with an optical index of the cylinders equal to n = 1.5, which leads to a coefficient of reflection of the order of 4 % per air/glass interface. In Fig: S1 of the supplementary materials(I B) [40], we show the propagation of a laser beam computed by solving rigorously the Maxwell equations into the photonic crystal. Because of the low optical-index contrast between the cylinder and the surrounding medium, light propagation is dominated by the transmitted rays. Thus, the field intensity distribution can be predicted from the distribution of the transmitted rays. So we focus only on rays which are refracted and thus transmitted through the cylinders. Within the assumption of a low index-contrast, a geometrical ray that undergoes a reflection can be considered as a new ray with new initial conditions. For the next steps, it will be preferentially transmitted through the cylinders and its motion will be described by the dynamical properties of the transmitted rays. The ray dynamics has a Hamiltonian formulation in terms of the eikonal equation: d ds n d γ ds = ∇n [START_REF] Landau | Mechanics. Third Edition[END_REF] where γ is the ray trajectory and s is the parametric distance along the ray (curvilinear abscissa). As a consequence, the phase space is four dimensional. Because n is constant inside the cylinder, the solutions to this equation, i.e. the ray trajectories, are piecewise linear. In order to characterize them, we compute the Poincaré surface of section. Poincaré surface of section are a convenient way to represent as a 2D-graph a higher dimensionality phasespace by representing a continuous evolution in a discrete way. They preserve informations about the dynamics [START_REF] Cvitanović | Discrete Time Dynamics[END_REF].

As an example, periodic orbits in the phase space appears as periodic points in the Poincaré surface of section [START_REF] Ozorio De | Hamiltonian Systems Chaos and Quantization[END_REF].

For the refractive Lorentz gas, the Poincaré surface of section is the set of points (θ Ap , θ A 1 p ) defined by the intersection of the ray with the boundary of a cylinder: θ Ap is the angle between the horizontal axis and the light ray when it crosses the boundary of the cylinder and θ A 1 p is the angle between the incident ray and the normal to the cylinder boundary (see Fig: 1). The Poincaré map is constructed in three steps :

1. An incident ray intersects a cylinder, defining a point A p , 2. The ray is refracted twice at the boundary, giving the point B p . The following set of equations describes these refraction steps.

θ 1 Bp = -θ 1 Ap ( 1
)
θ Bp = θ Ap + 2 arcsin   sin θ 1 Ap n   ( 2
)
where n is the optical index of the cylinders.

3. The ray propagates in air until it hits a new cylinder, this gives A p+1 and the process goes on. The couple (θ Ap+1 , θ 1 Ap+1 ) is found with the help of the equations:

f B (x) = tan(α) x -R cos(θ Bp ) + R sin(θ Bp ) (3) (x A (p+1) -x c p+1 ) 2 + (f B (x A (p+1) ) -y c p+1 ) 2 = R 2 (4)
where α is the angle between the horizontal axis and the ray as depicted in Fig: 1-b). It satisfies:

α = θ Bp + θ 1 Bp = θ A (p+1) + θ 1 A (p+1) . (x c
p+1 , y c p+1 ) are the cylinder coordinates at the (p + 1) th collision.

These equations are solved numerically in order to compute the ray trajectory in the whole crystal. Examples of Poincaré surfaces of section are given in Fig: 2 for different values of the normalized period t (t varying from 2.2 to 6). It will be shown in the following that this parameter governs the stability properties of the ray dynamics.

The most prominent point is the existence of islands of stability characterized by closed curves in the Poincaré surface of section, surrounded by a chaotic sea. Islands of stability characterize ballistic trajectories in the crystal, i.e. rays that oscillate around the axis of symmetry of the crystal. Cylinders act as cylindrical lenses: they focus rays in the vicinity of the axis of symmetry of the crystal. As a result, these trajectories are periodic or quasi-periodic in real space as explained in the supplementary materials (III) [40] and illustrated by the figures Fig: S6 and Fig: S7. The area covered by these islands of stability decreases when t increases and they reduce to single points above the threshold t thres ≥ 2n n-1 as shown in the supplementray materials IV [40]. In the example of Fig: 2, the threshold is t thres = 6 for n = 1.5. At small angles, the islands of stability are ellipsis. This corresponds to Gauss conditions: sin(θ) ∼ θ and cos(θ) ∼ 1. When the angles increase, the islands of stability are deformed because of non-linear effects, although rays still propagate along the crystal axis. At the edges of the transition between the islands of stability and the chaotic sea, Cantori [START_REF]small islands of quasi-periodic trajectories showing self-similar patterns[END_REF] can be observed (see Fig: 3) because of an increase of non-linear effects: this is a direct consequence of KAM theory [START_REF] Ozorio De | Hamiltonian Systems Chaos and Quantization[END_REF]. Again rays propagate along the crystal axis. When nonlinear effects still increase, rays leave crystal axis and the motion is completely chaotic, with sometimes segments of ballistic motion. In that case rays diffuse in the entire crystal. We conclude that the "quasi-integrable" character of the light dynamics in a photonic crystal [22] is a consequence of the dynamical properties of the transmitted rays since the dynamics of the refracted rays, which is similar to the dynamic of the Lorentz gas, is always chaotic [START_REF] Gaspard | Chaos, scattering and statistical mechanics[END_REF]. We want to point out that this conclusion is independent of the index-contrast between the cylinder and the surrounding.

The strength of chaos, i.e. the sensitivity to initial conditions, can be quantified through the computation of Lya- punov exponents. In Hamiltonian systems with two degrees of freedom, Lyapunov exponents come in pairs with a null sum because of area preservation in the phase-space [39] (this is Liouville theorem). The positive exponent denoted λ, and called in the following Lyapunov exponent, quantifies the exponential sensitivity to initial conditions. The Lyapunov exponents λ were computed per step, that is, each time a ray hits a cylinder. Because of islands of stability, the dynamics is clearly not ergodic and the Lyapunov exponent may depend on initial conditions. The Lyapunov exponent takes low value for ballistic paths, as expected (blue color in the lower inset of Fig: 4) and the value sharply increases outside the islands of stability. As shown in Fig: 2, the size of the islands of stability decreases as t increases.

Above this threshold, the Lyapunov exponent is almost constant and large, whatever the initial conditions are. In Fig: 4, the Lyapunov exponents averaged over the initial conditions are represented (plain line). The averaged exponents increase monotonically when increasing the normalized period t . No apparent discontinuity in the behavior of the Lyapunov exponent is seen near the threshold. Its value is smaller than the Lyapunov exponent for the reflective Lorentz gas (dashed line in Fig: 4) but their behavior with t is similar. As a matter of fact for large normalized period, both of them grow as λ ∼ ln(t ) [see [START_REF] Gaspard | Chaos, scattering and statistical mechanics[END_REF][START_REF] José | [END_REF] for the reflective Lorentz Gas and the supplementary materials (V) for the refractive Lorentz Gas]. This is due to the fact that for large t refractions at cylinders become uncorrelated and the distance between two refractions is given by the mean free path which is identical for both the reflective and the refractive periodic Lorentz gases. As an example, the averaged Lyapunov exponents is λ = 16.7 for t = 6 (see Fig: 4). This value means that, starting with an initial accuracy of ∆θ 1 0 = 10 -14 radian for the initial angle, after hitting only two cylinders the uncertainty on the angle becomes ∆θ 1 2 = 10 -14 × e 2×16.7 = 3.2. Thus, after hitting only two cylinders, the value of the angle θ 1 could be any value in [-π, π]. We have performed experiments as shown in the supplementary materials (VI) [40]. They clearly demonstrate this exponential sensitivity to the initial conditions that can lead to the diffusion of a laser beam into the photonic crystal.

Fig. 4: Lyapunov exponent versus the normalized period t of the crystal. The plain line is the Lyapunov exponent for the refractive Lorentz gas whereas the dashed blue curve corresponds to the reflective Lorentz gas. Lyapunov exponents are given per step unit. The inset shows the Lyapunov exponent for the refractive Lorentz gas as a function of t and of the initial angle between the ray and the horizontal axis (launching angle). It takes low values (blue color) for launching angles around the crystal axis of symmetry. For t ≥ t thres = 6 the Lyapunov exponent is independent of the launching angle. Diffusive properties of the photonic billiard. -From an experimental point of view, a light beam can be seen as an ensemble of rays with a spatial extension. The rays constituting the light beam have different directions of propagation. As a consequence, from an experimental point of view, an incoming light beam can be seen as an ensemble of rays with different initial conditions {θ A (0), θ 1 A (0)}. Because of the high values of the Lyapunov exponent, a light beam will undergo a diffusive propagation inside the photonic crystal. In order to characterize the diffusive properties of light in the photonic crystal, the propagation of a bundle of rays is computed. This simulation models the emission of a point source (e.g. a molecule or a quantum dot) inside the photonics crystal. As we will explain in the following our numerical results also provide informations about the propagation of a light beam inside the photonic crystal through the computation of the exponent of the diffusion law. In the following, we first show how rays are distributed into the photonic crystal by computing the density of rays. Then we compute the time evolution of the mean-square displacement. The simulation is performed as followed. The rays are initially launched at the center of the photonic crystal with random initial conditions (under a uniform distribution). After N max steps (i.e. each ray hints N max times a cylinder), the density of rays (i.e. the number of rays per unit area) is plotted. For low values of t (see Fig: 5-a) for t = 3) the dynamics is barely chaotic with stability islands covering a large part of the Poincaré surface of section (Fig: 2) and characterized by a lowvalue of the Lyapunov exponent (Fig: 4). In that case one can clearly notice a non-isotropic density of rays. Most of the rays concentrate in 8 spots localized on the principal axis of the photonics crystal. These rays propagate in the photonic crystal through regular trajectories. When t is above the threshold t ≥ t thres , the dynamics only admits chaotic motion (see Fig: 5-c) with t = 6) and the density of rays becomes isotropic with a maximum at the center of the photonic crystal. The density of rays follows a Gaussian distribution law, which indicates that rays undergo a brownian motion. For intermediate values of the normalized period t , the density of rays is characterized by spots on the crystal axis plus an isotropic distribution which is maximum at the photonic crystal center. For this intermediate situation, some rays follow regular paths which lead to the spots in the rays density and some others follow chaotic motion giving the isotropic distribution.

We now describe the computation of the evolution of the mean-square displacement with time. Hamiltonian systems are known to lead to superdiffusion [26]. It has been shown that in a two-dimensional periodic potential [26,27] the mean square-displacement varies as < r 2 n >∼ t α as a results of ballistic flights, random motion and the "sticky" barriers formed by Cantori [START_REF] Denisov | [END_REF]. From an experimental point of view, let us consider a layer of width L of a diffusive material. The transmission coefficient ,i.e. the fraction of light that is transmitted through the layer of material, can be related to the exponent of the diffusion law [31]. Indeed, it has been shown that the transmission coefficient is of the form T = 1 1+AL (3-α)/2 where L is the thickness of the layer of diffusive material, α the exponent of the diffusion law and A a constant. As a consequence, computing the diffusion law in the photonic crystal gives informations on the repartition of the light emitted by a point-source but also to the transmission of a light beam through a finite size photonic crystal. In order to find the exponent of the diffusion law α, we compute the mean square-displacement for an ensemble of 10 4 rays with initial conditions chosen randomly (with a uniform distribution). The simulation time is long enough to reach the stationary regime, following the prescription of ref. [38]. The exponent of the diffusion law is plotted on Fig: 6 as a function of the normalized period t . Its values are larger than 1 when the dynamics contains chaotic regions and islands of stability, which confirms the superdiffusive behavior. It reaches values close to 2 when the gap between the cylinders is small (t ∼ 2), i.e. when the dynamics is dominated by islands of stability and it decreases to 1 as the area of stability islands decreases. For t ≥ 6, that is to say when the dynamics is completely chaotic, diffusion is close to normal. Nevertheless, the exponent remains slightly larger than 1 because of ballistic paths exactly along the crystal principal axis and because of infinitehorizon trajectories (i.e. light rays that can move arbitrarily far without crossing a cylinder). Our results are in strong contrast with the periodic Lorentz gas. Indeed, the periodic Lorentz gas exhibits a weak form of superdiffusion [38] with a logarithmic correction to the normal diffusion [28,38]. Moreover the diffusion exponent does not depend on the period of the periodic Lorentz gas.

Conclusion. -In the high frequency limit, light propagation in a photonic crystal exhibits complex propagation patterns with regular paths and chaotic motion. For this limiting case, light transport is best described as a diffusive process that can be tuned from a superdiffusive regime to a nearly-brownian motion. An interesting feature of the system is its simplicity, allowing for experimental investigations of wave chaos footprint in optics.

Fig. 2 :

 2 Fig. 2: Poincaré surface of section for four values of the crystal photonic period T normalized by the circle radius R, t = T /R.They result from 100 initial conditions, each evolved for 200 steps. For small normalized period (i.e. low density) the Poincaré surface of section are dominated by ballistic paths (t = 2.2) whereas the dynamics is completely chaotic for t higher than t = 6 here. The area of islands of stability decreases when increasing t (examples are given for t = 3 and t = 5). For t = 5 islands around θA = π/4[π] no more exist because t √ 2 ≥ 2n n-1 as explain in the text.

1 Fig. 3 :

 13 Fig. 3: Details of Fig:2 for t = 3 showing Cantori surrounded by the chaotic sea.

Fig. 5 :

 5 Fig. 5: Rays density after 300 simulation steps for a) t = 3, b) t = 5 and c) t = 6. The simulation was running with 3.10 4 rays starting from (0,0) with random initial conditions under a uniform distribution.

Fig. 6 :

 6 Fig. 6: Exponent of the diffusion law defined by the law < (r -r0) 2 >= Dαt α versus the normalized period t of the crystal. The inset is a zoom on the interesting region where islands of stability exist t ∈ [2, 6].10 4 rays was used to compute the exponent of the diffusion law.


P Dreyfuss

Results for a turbulent system with unbounded viscosities: weak formulations, existence of solutions, boundedness, smoothness

Keywords: degenerate elliptic system, weak solution, regularity, oceanography. 1991 MSC: 35J70, 35D05, 35D10, 86A05

We consider a circulation system arising in turbulence modelling in fluid dynamics with unbounded eddy viscosities. Various notions of weak solution are considered and compared. We establish existence and regularity results. In particular we study the boundedness of weak solutions. We also establish an existence result for a classical solution.

Introduction

Let Ω be an open bounded set in R 3 , with a Lipschitz boundary. We consider the following turbulent circulation model : Email address: Pierre.Dreyfuss@iecn.u-nancy.fr (P. Dreyfuss). URL: www.iecn.u-nancy.fr/~dreyfuss (P. Dreyfuss).

(P)                    -div(ν(k)∇u) = f in Ω -div(a(k)∇k) = ν(k)|∇u|
We study Problem (P ) under the following main assumption:

(H 0 )              f ∈ L r (Ω), with r > 3 2
a, ν : R + → R + are continuous ∃ δ > 0 : a(s), ν(s) ≥ δ ∀s ∈ R + Problem (P) is a simplified scalar version of the RANS model arising in oceanography (see [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF][START_REF] Lederer | A RANS 3D model with unbounded eddy viscosity[END_REF][START_REF] Bernardi | A model for two coupled turbulent fluids Part III: Numerical approximation by finites elements[END_REF]): the function u is an idealisation of the mean velocity of the fluid and k is the turbulent kinetic energy. The mathematical analysis of (P) is a step towards better understanding the RANS model. Various studies were made in this direction. Some existence results were established in [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF][START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF].

In this paper we focus on the case where the viscosity functions a and ν are not a priori bounded. In fact (see [START_REF] Lederer | A RANS 3D model with unbounded eddy viscosity[END_REF][START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF]), in the relevant physical situation, we have

(H p )      a(s) = a 1 + a 2 √ s ν(s) = ν 1 + ν 2 √ s
We will establish an existence result for a weak solution for (P) under less restrictive assumptions than in [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF]. An important feature is that our assumptions are satisfied under (H p ), contrarily to the assumptions made in [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF]. Moreover we give additional regularity results for the weak solution we obtain. In particular, under (H 0 ) and the following additional assumption: a is proportional to ν, ∂Ω is of class C 2,α , f ∈ C 0,α (Ω) and ν ∈ C 1,α (R + ), we prove the existence of a classical solution for (P).

We also compare our results with the results presented in [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF].

Another feature of our work is to considere various notions of weak solution for Problem (P): W -solution, H-solution, distributional solution, renormalized solution, 'energy solution', classical solution. We give some relations between these notions.

Notions of weak solution for (P)

We can reformulate equation (P).2 by using the Kirchoff transform. Let

A(s) := s 0 a(t)dt.
Instead of (P).2, we can consider

(P ).2 ′ -∆K = ν • A -1 (K)|∇u| 2 on Ω, where K = A(k).
In fact, from every distributional solution K ∈ W 1 (Ω) of (P).2' we obtain a distributional solution k of (P).2 by setting k = A -1 (K). This property is related to the facts that A is invertible, A -1 (0) = 0 and |A -1 (s)| ≤ C.s (this can be seen by using the assumptions made on ν in (H 0 )).

The situation is more complicated for equation (P).1, where the a priori unbounded coefficient ν(k) appears in the principal part of the operator and cannot be removed.

Hence we have to restrict u to satisfy the energy condition

Ω ν(k)|∇u| 2 < ∞. (1)
Nevertheless we will see later on that various non equivalent notions of weak solution can be considered for (P).1.

We will introduce the notions of W-solution and H-solution. It is also possible to consider the notion of renormalized solution (see [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF] chap.5). In [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF] the authors defined another notion that they call energy solution.

We will give some relations between these notions in the Appendix I.

Remark now that under the restriction (1), the right hand side in (P).2 (or in (P).2') is only a priori in L 1 (Ω). Hence (see [START_REF] Boccardo | Non-linear elliptic or parabolic equations involving measure date[END_REF]) it is natural to seek k in the space ∩ p<3/2 W 1,p 0 (Ω). We want to find a function u vanishing on ∂Ω that satisfies the energy condition [START_REF] Bernardi | A model for two coupled turbulent fluids Part III: Numerical approximation by finites elements[END_REF]. This leads to considering the following spaces:

W k = v ∈ H 1 0 (Ω) : [v] k < ∞ H k = closure of C ∞ c (Ω) with respect to [.] k
where we used the notation

[v] k = Ω ν(k)|∇v| 2 1/2 .
For any measurable function k, the map [.] k defines a norm on W k . In the general situation H k and W k are not equal. Moreover W k is not necessarily complete and a function in H k does not always have a uniquely defined gradient (see [START_REF] Zhikov | Weighted Sobolev spaces[END_REF]). If we assume that ν(k) ∈ L 1 (Ω) then W k is complete and in fact H k ⊂ W k are Hilbert spaces (see [START_REF] Dreyfuss | Higher integrabiblity of the gradient in degenerate elliptic equations[END_REF][START_REF] Zhikov | Weighted Sobolev spaces[END_REF][START_REF] Cassano | On the local boundedness of certain solutions for a class of degenerate elliptic equations[END_REF]) when they are equipped with the scalar product

(v, w) = Ω ν(k)∇v∇w.
Consequently, we will consider the following two distinct notions of solution for (P).1:

u is called a H k -solution of (P).1 if u ∈ H k and Ω ν(k)∇u∇v = Ω f v ∀v ∈ H k u is called a W k -solution of (P).1 if u ∈ W k and Ω ν(k)∇u∇v = Ω f v ∀v ∈ W k
Finally, we define the following notions of weak solution for (P):

(u, k) is called a H-solution of (P) if k ∈ ∩ p<3/2 W 1,p 0 (Ω), u ∈ H k , k is a distributional solution of (P).2 and u is a H k -solution of (P).1 (u, k) is called a W-solution of (P) if k ∈ ∩ p<3/2 W 1,p 0 (Ω), u ∈ W k k is
a distributional solution of (P).2 and u is a W k -solution of (P).1

Main results

Let (H 1 ) and (H 2 ) denote the following conditions:

(H 1 ) ∃γ > 0 : a(s) ≥ γν(s) ∀s ∈ R + (H 2 ) ∃γ > 0 : a(s) = γν(s) ∀s ∈ R +
We will establish: Theorem 1 Assume that (H 0 ) and (H 1 ) hold. Then there exists at least one Wsolution (u, k) for (P) such that

u ∈ L ∞ (Ω) and Ω a(k)|∇k| 2 < ∞.
(

) 2
Corollary 1 Assume that in addition to (H 0 ) and (H 1 ) we have

∃ ν 0 > 0 : ν(s) ≤ ν 0 (1 + s 6 ), ∀s ∈ R + . (3)
Then the W-solution (u, k) given in Theorem 1 is a distributional solution of (P).

Theorem 2 Assume that (H 0 ) and (H 2 ) hold. Then the W-solution (u, k) given in Theorem 1 satisfies u, k ∈ C 0,α (Ω), for some α ∈ (0, 1). ( 4)

Moreover (u, k) is also a H-solution of (P) (and in fact a classical weak solution). If in addition to (H 0 ) and (H 2 ) we assume that ∂Ω is of class C 2,α , f ∈ C 0,α (Ω) and ν ∈ C 1,α (R + ) then u, k ∈ C 2,β (Ω), for some β ∈ (0, 1), and (u, k) is a classical solution of (P).

Discussion of the results

In Theorem 1 we give an existence result of a W-solution. We next give some regularity results: firstly the property (2) and secondly (in Theorem 2) the property (4). Finally, in Theorem 2 we give an existence result for a classical solution for (P).

The main previous studies of Problem (P) are presented in [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF] chap. 5 and in [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF].

In [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF] chap.5, the authors prove the existence of a renormalized solution for (P) under the assumptions (H 0 ) and (H 2 ). It seems that their proof also works under (H 0 ) and (H 1 ). Nevertheless the notion of renormalized solution is very weak. A renormalized solution (u, k)

for (P) is a distributional solution if ν(k) ∈ L ∞ (Ω), whereas a H-or a W-solution is a distributional solution if ν(k) ∈ L 1 (Ω) (see the Appendix I).
In [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF] the authors introduced a notion of solution that they call 'energy solution' (see the Appendix I). In fact an 'energy solution' is a W-solution which satisfies an additional property ensuring that H k = W k (the additionnal property imposed is sufficient but not necessary to have this equality). Under this point of view an 'energy solution' is slighty stronger than a W-solution. However, their existence result is obtained by assuming complicated conditions on the coefficients a and ν which are not exactly satisfied in the physically relevant situation (H p ), but only in the following approximate situation:

(H ′ p )             
for some ǫ > 0 we have:

a(s) = a 1 + a 2 √ s + ǫ ν(s) = ν 1 + ν 2 √ s + ǫ
On the contrary, our assumptions in Theorem 1 and Corollary 1 are very simple, and they are satisfied in (H p ).

Note also that we establish the regulartity property (2) which are not established in [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF] (or in [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF]).

In the Appendix I we also give a new existence result for an 'energy solution'.

In Theorem 2 we assume that (H 0 ) and (H 2 ) hold. These assumptions are fulfilled in the physical situation (H p ) if a 2 ν 1 = a 1 ν 2 . We then prove that u and k are Hölder continuous. In particular we give here a positive answer to a central question put in [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF] : k is bounded. Note that in this situation we clearly have

W k = H k .
We next establish the existence of a classical solution for Problem (P) by assuming some differentiability properties for a and ν. These properties are fulfilled in the situation (H ′ p ) if a 2 ν 1 = a 1 ν 2 . It seems that this result is completely new: the existence of a classical solution for (P) was not studied in any previous work.

Organization of the paper

In the sequel n will always denote an arbitrary integer greater or equal to one, and C (possibly with subscript) will denote a positive real that does not depend on n, but that can differ from one part to another. We always consider the space H 1 0 (Ω) equipped with the gradient norm. The condition (H 0 ) is always assumed.

• In section 2 we introduce an approximate sequence (u n , k n ) of solutions obtained by truncating the coefficients a and ν. We immediatly obtain the basic estimates :

Ω ν n (k n )|∇u n | 2 ≤ C ∀p < 3 2 : Ω a n (k n )∇k n p ≤ C
The point is that we establish the following fundamental estimates:

u n L ∞ (Ω) ≤ C Ω a n (k n )|∇k n | 2 ≤ C ( * )
The first estimate above is proved by developping further a technique due to Stampacchia.

The second is obtained under the assumption (H 1 ). The proof is based on the following idea: if (u, k) is a solution of (P), we formally have

1 ν(k)|∇u| 2 = -div(ν(k)∇u).u =f u + div(ν(k)u∇u). (5)
In other words one can hope that the second member in the second equation in (P) is more regular than it seems.

In fact, we prove that a similar relation to (5) holds for the approximate sequence. By using next that (u n ) is uniformly bounded in L ∞ (Ω), we obtain ( * ) which is the key estimate to prove Theorem 1.

• In section 3 we extract from (u n , k n ) a subsequence converging to some element denoted by (u, k). Under the assumptions (H 0 ) and (H 1 ), we directly obtain that

u ∈ H 1 0 (Ω) ∩ L ∞ (Ω), k ∈ H 1 0 (Ω).
We prove that moreover we have:

Ω ν(k)|∇u| 2 < ∞, Ω a(k)|∇k| 2 < ∞.
• In section 4 we pass to the limit in the approximating Problems. In a first step we prove that u is a W k -solution of (P).1. To do this, we use the test functions v = h q (k n )ϕ (where ϕ ∈ W k ∩ L ∞ (Ω) and (h q ) is a sequence of functions that cut off the large values), and we pass to the limits n → ∞, q → ∞.

We next prove that the energies of the approximating sequence converge to the en-

ergy Ω ν(k)|∇u| 2 .
Finally we can pass to the limit in the second equation in order to prove that k is a distributional solution of (P).2. We then obtain Theorem 1 and Corollary 1 follows.

• In section 5 we assume that (H 0 ) and (H 2 ) hold. In a first step we obtain the estimate k n L ∞ (Ω) ≤ C. Hence k ∈ L ∞ and by using the De Giorgi-Nash Theorem we prove the Hölder continuity of u and k. Next, by assuming additional regularity on ν, ∂Ω and f we can apply the Schauder's estimates and we prove Theorem 2.

• In the Appendix I we study some relations between the notions of W -solution, H-solution, distributional solution, renormalized solution and 'energy solution' for Problem (P). We continue the discussion begun in Subsection 1.3 and we also establish a new existence result for an 'energy solution' for Problem (P). In the Appendix II we recall some basic properties of Hölder continuous functions.

Approximating sequence and estimates

We assume that (H 0 ) holds and we set

ν n (s) = T n (ν(s)) (6) a n (s) = T n (a(s)), ( 7
)
where T n is the truncated function defined by T n (t) = min(n, t).

We consider the Problem of finding (u n , k n ) ∈ (H 1 0 (Ω)) 2 such that

(P n )      Ω ν n (k n )∇u n ∇v = Ω f v ∀v ∈ H 1 0 (Ω) Ω a n (k n )∇k n ∇ϕ = Ω T n ν n (k n )|∇u n | 2 ϕ ∀ϕ ∈ H 1 0 (Ω) For any n ≥ 1, Problem (P n ) is well posed because a n , ν n ∈ L ∞ (R) and a -1 n , ν -1 n ∈ L ∞ (R) by construction.
It is proved in [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF] that a solution (u n , k n ) exists for any n ≥ 1. Moreover, the following basic properties were established:

k n ≥ 0 (8) Ω ν n (k n )|∇u n | 2 ≤ C 1 (9) ∀p < 3 2 : Ω a n (k n )∇k n p ≤ C 2 (10)
We now establish Lemma 3 The sequence u n is uniformly bounded in the L ∞ (Ω)-norm, that is,

u n L ∞ (Ω) ≤ C 3 (11)
Before proving this lemma we point out that the assumption f ∈ L r (Ω), with r > 3 2 made in (H 0 ) implies that

f ∈ W -1,ρ (Ω), with ρ = 3r 3 -r > 3. ( 12
)
This last property is easy to prove by using the Sobolev injection Theorem.

Proof

We will obtain the estimate (11) by using the technique presented on p.108 in [START_REF] Stampacchia | Equations elliptiques du second ordre à coefficients discontinus[END_REF].

In order to prove that C 3 is independent of n we have to detail the technique of

Stampacchia. Let b n (u, v) := Ω ν n (k n )∇u∇v.
Recall that f satisfies (12) and then by using a classical result (see [START_REF] Brezis | Analyse fonctionnelle, théorie et applications[END_REF]) there exists g ∈ (L ρ (Ω)) 3 such that -div(g) = f and g

(L ρ (Ω)) 3 ≤ C f L r (Ω) . Hence the sequence u n satisfies b n (u n , v) = Ω g∇v ∀v ∈ H 1 0 (Ω). ( 13
)
For s ≥ 0, we define the measurable set A n (s) ⊂ Ω by setting

A n (s) = x ∈ Ω : |u n (x)| ≥ s .
We also introduce ϕ := max (|u n | -s, 0) sgn(u n ). ( 14) It is proved in [START_REF] Stampacchia | Equations elliptiques du second ordre à coefficients discontinus[END_REF] that ϕ ∈ H 1 0 (Ω) and

∇ϕ = ∇u n in A n (s) ∇ϕ = 0 in Ω \ A n (s)
By testing [START_REF] Cassano | On the local boundedness of certain solutions for a class of degenerate elliptic equations[END_REF] with v = ϕ, we obtain

b n (ϕ, ϕ) = b n (u n , ϕ) = An(s) g∇ϕ. (15)
Remark now that assumption ν(s) ≥ δ > 0 in (H 0 ) implies that ν n (k n ) ≥ min(δ, 1). Consequently the bilinear form b n is uniformly coercive on H 1 0 (Ω). By using this property together with the Hölder inequality, we obtain from (15):

ϕ 2 H 1 0 (Ω) ≤ C An(s) |g| 2 1/2 ϕ H 1 0 (Ω) .
Hence by using the Cauchy inequality together with the Hölder inequality we obtain

ϕ 2 H 1 0 (Ω) 2 ≤ C1 g 2 L ρ (Ω) A n (s) ρ-2 ρ .
(16) On the other hand, the Poincaré-Sobolev inequality gives An(s)

|ϕ| 6 1/3 ≤ C2 ϕ 2 H 1 0 (Ω) . ( 17
)
Let now t > s. It is clear that A n (t) ⊂ A n (s) and consequently

An(s) |ϕ| 6 1/3 ≥ An(t) |ϕ| 6 1/3 ≥ An(t) |t -s| 6 1/3 ≥ |t -s| 2 A n (t) 1/3 (18)
We set ψ n (s) := A n (s) , ∀s ≥ 0 For fixed n, ψ n is a decreasing function, and from the estimates ( 16)-( 18), we obtain

ψ n (t) ≤ C3 |ψ n (s)| β (t -s) -6 ∀t > s ≥ 0,
where we have used the notation β := 3(ρ-2) ρ > 1 and where C3 = C3 ( C1 , C2 , f L r ). Both quantity β and C3 do not depend on n. Hence by using Lemma 4.1 in [START_REF] Stampacchia | Equations elliptiques du second ordre à coefficients discontinus[END_REF] it follows:

ψ n (θ) = 0, where θ = 2 β/(β-1) C3 |Ω| β-1 1/6 < ∞ does not depend on n.

This property tells precisely that (11) holds true with C 3 = θ.

Notice that the bilinear form

(u, v) → Ω a n (k n )∇u∇v,
is also uniformly coercive on H 1 0 (Ω). Moreover, the sequence

h n := T n ν n (k n )|∇u n | 2
is imbedded in L ∞ (Ω). We can then apply again the technique of Stampacchia detailed in the proof of lemma 3, and obtain:

for n ≥ 1 : k n ∈ L ∞ (Ω) ( 19
)
Nevertheless the control we have on {h n } is obtained from ( 9), which gives a uniform bound in the L 1 -norm for the sequence. This is not enough to obtain a uniform estimate for {k n } in the L ∞ -norm. However we can establish:

Lemma 4 Assume that (H 0 ) and (H 1 ) hold. Then we have

a n (s) ≥ γ 1 ν n (s), γ 1 = min(1, γ) (20) Ω a n (k n )|∇k n | 2 ≤ C 5 (21)

Proof

The estimate (20) is easy to obtain. Its verification is left to the reader. Let (u n , k n ) be the chosen approximating sequence. We have from ( 11) and ( 19) that

∀n ≥ 1 : u n , k n ∈ H 1 0 (Ω) ∩ L ∞ (Ω)
It follows (see [START_REF] Brezis | Analyse fonctionnelle, théorie et applications[END_REF])

that v := u n .k n ∈ H 1 0 (Ω) ∩ L ∞ (Ω) is admissible for (P n ).1 and we get 2 Ω ν n (k n )|∇u n | 2 k n = Ω f u n k n - Ω ν n (k n )u n ∇u n ∇k n (22)
By testing (P n ).2 with ϕ = k n , we obtain:

Ω a n (k n )|∇k n | 2 = Ω T n ν n (k n )|∇u n | 2 k n ≤ Ω ν n (k n )|∇u n | 2 k n , (23)
by using the properties T n (s) ≤ s and (8). Hence, by combining ( 22) with (23) we have:

I := Ω a n (k n )|∇k n | 2 ≤ Ω |f u n k n | :=II + Ω ν n (k n )u n ∇u n ∇k n :=III (24)
We can estimate the term II as follows:

2 more generally:

ν n (k n )|∇u n | 2 = -div(ν n (k n )∇u n ).u n =f un +div(ν n (k n )u n ∇u n ) in D ′ (Ω) II ≤ C 3 Ω |f k n | Hölder Ineq. ≤ C 3 f L 3/2 k n L 3 Poincaré-Sobolev Ineq. ≤ C1 f L 3/2 Ω |∇k n | 2 1/2 ≤ C1 δ f L 3/2 Ω a n (k n )|∇k n | 2 1/2 Young Ineq. ≤ C1 δ 1 ǫ f 2 L 3/2 + ǫ Ω a n (k n )|∇k n | 2 for any ǫ > 0 given ≤ 1 3 Ω a n (k n )|∇k n | 2 + C2 f 2 L 3/2
where δ > 0 is the constant given in (H 0 ). The last inequality was obtained by choosing ǫ = δ/(3 C1 ), using the estimate [START_REF] Lederer | A RANS 3D model with unbounded eddy viscosity[END_REF] and by setting C2 = 3 C1 2 /δ 2 .

We next estimate the term III:

III = Ω u n ν n (k n )∇u n ν n (k n )∇k n ≤ C3 Ω ν n (k n )∇u n a n (k n )∇k n , C3 = C 3 γ -1/2 1 ≤ 1 3 Ω a n (k n )|∇k n | 2 + C4 Ω ν n (k n )|∇u n | 2 , C4 = C4 ( C3 )
where C 3 , γ 1 are the constants that appear in [START_REF] Lederer | A RANS 3D model with unbounded eddy viscosity[END_REF] and (20). The last inequality follows from the Young inequality.

Recall now the inequality (24) and use the estimates established for the terms II and III. We obtain:

1 3 Ω a n (k n )|∇k n | 2 ≤ C2 f 2 L 3/2 (Ω) + C4 Ω ν n (k n )|∇u n | 2 . (25)
By using (25) together with (9) we finally obtain (21).

3 Basic convergence results for (u n , k n )

The estimates established in the previous section allow us to extract a converging subsequence from (u n , k n ). We have Lemma 5

1. Assume that (H 0 ) holds. Then we can extract a subsequence (still denoted by

(u n , k n )) such that a n (k n )∇k n ⇀ a(k)∇k in L p (Ω), p < 3 2 (26)
k n → k a.e in Ω (27) u n ⇀ u in H 1 0 (Ω) (28) u n * ⇀ u in L ∞ (Ω) (29)
2. If in addition the condition (H 1 ) is fulfilled then we may assume that

k n ⇀ k in H 1 0 (Ω) (30) Proof 1.
The properties (26) and ( 27) are obtained from [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF]. The property (28) is obtained by using the estimate (9) together with the assumption ν(s) ≥ δ > 0 in (H 0 ). We establish (29) from the estimate [START_REF] Lederer | A RANS 3D model with unbounded eddy viscosity[END_REF].

2. By using Lemma 4 together with the assumption a(s) ≥ δ > 0 in (H 0 ) we obtain (30). Notice that the k appearing in (26), ( 27) and ( 30) is necessarily the same in the three situations.

We are able to prove additional regularity results for the element (u, k) introduced in Lemma 5. For technical reasons we introduce the sequence {h q } q∈N of real functions defined in [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF] p. 185. It satisfies:

|h q (s)| ≤ 1 ∀(q, s) ∈ N × R (31) h q (s) = 0 when |s| > 2q (32) |h ′ q (s)| ≤ 1 q ∀q ∈ N,
a(k)|∇k| 2 < ∞ (36) Proof 1.
We take over the arguments presented in [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF] p. 192. For q ≥ 1, we set

η n,q := h q (k n )ν n (k n ) 1/2 ∇u n
Let now q be fixed. The sequence

h q (k n )ν n (k n ) 1/2 n≥1
is uniformly bounded in L ∞ (Ω). Consequently, {η n,q } n≥1 is bounded in (L 2 (Ω)) 3 and we can extract a subsequence weakly convergent to some η q ∈ (L 2 (Ω)) 3 .

On the other hand, we have

h q (k n )ν n (k n ) 1/2 → h q (k)ν n (k) 1/2
a.e in Ω ∇u n ⇀ ∇u in L 2 (Ω), and thus η q = h q (k)ν(k) 1/2 ∇u.

We now use a classical property of the weak convergence in L 2 (Ω):

η q L 2 (Ω) ≤ lim inf n→∞ η n,q L 2 (Ω) ≤ lim inf n→∞ Ω ν n (k n )|∇u n | 2 1/2 ≤ C 1/2 1 ,
where C 1 is a constant independent of q given in [START_REF] Ladyzenskaja | ceva: 'Equations aux Dérivées Partielles de Type Elliptiques[END_REF]. By using properties (34) and (31) we can see that

η 2 q → q→∞ ν(k)|∇u| 2 a.e. in Ω η 2 q ≤ ν(k)|∇u| 2
Hence by the Fatou Lemma we finally obtain:

Ω ν(k)|∇u| 2 ≤ lim inf q→∞ η q 2 L 2 ≤ C 1 2.
If the additional assumption (H 1 ) holds, then we have the estimate (21) and the previous reasoning allows us to obtain (36)

4 The proof of theorem 1

In the previous section we have proved that under (H 0 ) we can extract a converging subsequence of (u n , k n ). If moreover (H 1 ) holds then the limit (u, k) obtained satisfies:

u ∈ W k ∩ L ∞ (Ω) (37) k ∈ H 1 0 (Ω) (and in fact k ∈ W k ) (38)
4.1 Passing to the limit in (P n ).1

We recall that the space W k was defined by

W k = v ∈ H 1 0 (Ω) : [v] k < ∞
We now establish:

Lemma 7 Assume that (H 0 ) and (H 1 ) hold. Then the element (u, k) given in Lemma 5 satisfies (37), (38) and:

Ω ν(k)∇u∇v = Ω f v ∀v ∈ W k (39) Proof Let n ≥ 1, q ∈ N and ϕ ∈ W k ∩ L ∞ (Ω).
We consider the function v := h q (k n )ϕ. By recalling the properties (31)-(34) of h q , we can verify that h

q (k n ) ∈ H 1 0 (Ω) ∩ L ∞ (Ω). Consequently v ∈ H 1 0 (Ω) ∩ L ∞ (Ω)
. By testing (P n ).1 with v, we obtain:

I := Ω ν n (k n )h q (k n )∇u n ∇ϕ + Ω h ′ q (k n )ν n (k n )∇u n ∇k n ϕ :=II = Ω f h q k n ϕ :=III (40)
In a first step we fix q and we study the behaviour of terms I, II and III when n tends to infinity. By using the property (32) we see that

|ν n (k n )h q (k n )| ≤ max s∈[0,2q] ν(s) := C q ,
and by using (32) together with (27) we obtain

ν n (k n )h q (k n ) → n→∞ ν(k)h q (k) a.e in Ω.
Consequently

ν n (k n )h q (k n )∇ϕ → n→∞ ν(k)h q (k)∇ϕ in (L 2 (Ω)) 2 ,
and by also employing (28) we get:

I → n→∞ Ω ν(k)h q (k)∇u∇ϕ (41)
We now estimate II. From (33) we obtain:

II ≤ 1 q {q≤kn≤2q} ν n (k n )∇u n ∇k n ϕ ≤ ϕ L ∞ C q Ω ν n (k n )|∇u n | 2 1/2 Ω a n (k n )|∇k n | 2 1/2 ≤ C q , (42)
where the second inequality is obtained by using (20).

For the last term we get

III → n→∞ Ω f h q (k)ϕ (43)
By using the estimates (41)-(43) together with (40) we obtain that for any fixed ϕ ∈ W k ∩ L ∞ (Ω) the following holds Ω ν(k)h q (k)∇u∇ϕ

:=J 1 = Ω f h q (k)ϕ :=J 2 +O( 1 q ). ( 44
)
We next remark that the integrand in J 1 converges for a.e. x ∈ Ω to ν(k)∇u∇ϕ when q tends to infinity. Moreover by using (31) together with the fact that ϕ ∈ W k we can see that the integrand in J 1 is dominated by |ν(k)∇u∇ϕ| ∈ L 1 (Ω). Consequently, by the Dominated Convergence Theorem we get

J 1 → q→∞ Ω ν(k)∇u∇ϕ
Similarly we can see that

J 2 → q→∞ Ω f ϕ
At this stage we have proved that

Ω ν(k)∇u∇ϕ = Ω f ϕ ∀ϕ ∈ W k ∩ L ∞ (Ω), (45)
and it remains to show that the condition ϕ ∈ L ∞ (Ω) is not necessary. Let ϕ ∈ W k and i ∈ N. We consider ϕ i ∈ W k ∩ L ∞ (Ω) given by ϕ i = T i (ϕ). By using some basic properties of T i (see [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF]), we see that |ϕ i | ≤ |ϕ|, |∇ϕ i | ≤ |∇ϕ|, ϕ i → ϕ a.e , and ∇ϕ i → ∇ϕ a.e in Ω. Consequently, if we take ϕ i as test function in (45), we can pass to the limit i → ∞ and we obtain (39).

In Lemma 7 we have showed that u is a W k -solution of (P).1. In order to prove Theorem 1 we have to prove that k is a distributional solution of (P).2. We need first to establish:

Lemma 8 Assume that (H 0 ) and (H 1 ) hold. Then, in addition to the results presented in Lemma 5, we may assume:

ν n (k n )|∇u n | 2 → n→∞ ν(k)|∇u| 2 in L 1 (Ω) (46)
Proof We test (P n ).1 with the function u n . By using (28) we obtain:

Ω ν n (k n )|∇u n | 2 → n→∞ Ω f u = Ω ν(k)|∇u| 2 , ( 47
)
where the latter equality is obtained by testing (39) with u.

We set η n := ν n (k n )∇u n and η := ν(k)∇u. The relation (47) tells us that

η n L 2 (Ω) → n→∞ η L 2 (Ω) (48)
We can next take over the arguments presented in [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF] Lemma 5.3.4 in order to obtain:

η n ⇀ n→∞ η in (L 2 (Ω)) 2 (49)
Finally properties (49) and (48) imply that the convergence is strong in (49), and (46) follows.

  2 in Ω u = 0 on ∂Ω k = 0 on ∂Ω Here f, a and ν are given, and the functions u, k : Ω → R are the unknowns. ⋆ This work was partially supported by the Swiss National Science Foundation under Grant Number 200020-100051/1.

Assume that (H 0 ) and (H 1 ) hold. In Lemma 5 we have extracted a subsequence (u n , k n ) which converges in a certain sense to an element (u, k). This element has the properties (37)-(38). Next we have established (39). Let now ϕ ∈ C ∞ c (Ω). By using (26) we get:

We next remark that the property (46) ensures that

Recall that the sequence (u n , k n ) satisfies (P n ).2. Then relation (50) together with (51) allows to take the limit in (P n ).2. We get:

Thus (P).2 is fulfilled in the distributional sense.

At this point we have obtained (37), ( 38), (39) and (52). The proof of Theorem 1 is complete.

Assume now that the condition (3) in Corollary 1 is fulfilled. By using (38) together with the Sobolev Injection Theorem we get k ∈ L 6 (Ω) and thus ν(k) ∈ L 1 (Ω). Then we can conclude the proof of Corollary 1 by using Proposition 9 in the Appendix I: (u, k) is a distributional solution of (P).

The proof of Theorem 2

We assume in this section that (H 0 ) and (H 2 ) hold.

In this situation all the results presented in section 2 and section 3 are valid. For technical reasons we slightly modify the definition of a n by setting

where γ > 0 is the constant appearing in (H 2 ) and ν n is defined as before.

We will now consider Problems (P n ) modified by the new definition (53) of a n . Nevertheless the modification is very sligth, and all the results presented in the previous section can be recovered easly. The verifications are left to the reader.

We now prove that we have the new estimate:

In order to prove this result we set

and we remark that (P n ).2 leads to

and note that the sequence f u n is uniformly bounded in L r (Ω) with r > 3/2. These properties are sufficient (see the proof of Lemma 3) to get the estimate

where C does not depend on n.

The estimate (54) is finally obtained by using Lemma 3 together with (56). Consequently, in addition to the properties in Lemma 5 we may assume that

We will now prove that u, k ∈ C 0,α (Ω) for some α ∈ (0, 1). (58)

Let λ := ν(k). We have λ, λ -1 ∈ L ∞ (Ω) and

Recall also that f have the property [START_REF] Stampacchia | Equations elliptiques du second ordre à coefficients discontinus[END_REF]. Hence we can apply the De Giorgi-Nash Theorem (see for instance [START_REF] Dautray | Analyse mathématique et calcul numérique pour les sciences et techniques, volume 2, l'opérateur de Laplace[END_REF] Prop. 6 p.683 or [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF] Th. 8.22 and Th. 8.29). We obtain that u ∈ C 0,α 1 (Ω) for some α 1 ∈ (0, 1). We next set χ := k+(γ/2)u 2 . Then χ ∈ H 1 0 (Ω) and we have

By using the fact that u ∈ L ∞ (Ω) in (60), we can again apply the De Giorgi-Nash Theorem to get χ ∈ C 0,α 2 (Ω) for some α 2 ∈ (0, 1). Hence also k is Hölder continuous, and (58) follows.

Let α ∈ (0, 1) be a generic parameter that can differ from one part to another. We assume now that ∂Ω is of class C 2,α , f ∈ C 0,α (Ω) and ν ∈ C 1,α (R + ). We will prove the second part of Theorem 2 by iterating the Schauder estimates.

We have λ = ν(k) ∈ C 0,α (Ω) (see the Appendix II) and then, by applying the Schauder estimate (see [START_REF] Chen | Second Order Elliptic Equations and Elliptic Systems[END_REF] Theorem 2.7 p. 154) on (59) we get u ∈ C 1,α (Ω). Similarily, from equation ( 60) we obtain χ ∈ C 1,α (Ω) and thus k ∈ C 1,α (Ω). Hence (see Appendix II) λ ∈ C 1,α (Ω). By iterating again the Schauder estimates (see now Theorem 2.8 p.154 in [START_REF] Chen | Second Order Elliptic Equations and Elliptic Systems[END_REF]) we obtain that u and k are in C 2,α (Ω).

Finally we see that (u, k) is a classical solution of (P). Theorem 2 is proven.

Appendix I: Some relations between the notions of weak solution

We give here some relations between the various notions of weak solution: Wsolution, H-solution, distributional solution, renormalized solution, 'energy solution'.

Comparison with renormalized solution

We have:

then any W -or H-solution of Problem (P) is a also a distributional solution of (P). Proof 1. Let (u, k) be a W -solution of (P). Then the conditions (5.2.1)-(5.2.5) in [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF] chap.5 are satisfied. We have to prove that (5.2.6) holds.

By testing (39) with v we obtain the relation (5.2.6).a in [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF].

We remark that v is also admissible in (52). This allows us to obtain the condition (5.2.6).b in [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF]. Consequently (u, k) is a renormalized solution of (P).

If we consider a H-solution (u, k) of (P) we can take over the previous argument because the function v is now in

In consequence a W k -or a H k -solution of (P).1 is also a distributional solution of this equation. Hence (h, k) is a distributional solution of (P).

Remarks

1. The first point in Proposition 9 tells that the notions of H-or W -solution are stronger that the notion of renormalized solution. This fact is coherent with the second point established in Proposition 9: a H-or W -solution is a distributional solution if ν(k) ∈ L 1 (Ω) whereas a renormalized solution is only a priori a distributional solution if ν(k) ∈ L ∞ (Ω) (see [START_REF] Lewandowski | Analyse mathématique et océanographie[END_REF] p.185). 2. if we have k ∈ H 1 0 (Ω) and if ν satisfies the growth condition (3) then ν(k) ∈ L 1 (Ω).

Comparison with 'energy solution'

We have seen that when ν(k) ∈ L 1 (Ω) then any W -(or H-) solution is a distributional solution. Moreover the notion of W -solution coincides with the notion of H-solution iff W k = H k (see [START_REF] Zhikov | Weighted Sobolev spaces[END_REF]).

Some sufficient conditions to have this last equality were established in [START_REF] Zhikov | Weighted Sobolev spaces[END_REF] and in [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF], but necessary and sufficient conditions are not known.

Let us consider the following condition:

(R)

It was schown in [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF] that the first condition in (R) together with the property

In [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF] the authors introduced the notion of 'energy solution'. They impose (H 0 ) as the basic assumption. Then an 'energy solution' (u, k) for (P) is in fact a W -solution which satisfies (R). This implies that W k = H k . The energy solution is also a Hsolution, and moreover a distributional solution (because the first assumption in (R) implies that ν(k) ∈ L 1 (Ω)).

We see then that the notion of 'energy solution' (in the sense of [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF]) has the advantage of unifying various notions by putting us in the situation where ν(k) ∈ H 1 (Ω). The disadvantage is that we have to impose more complicated conditions on the coefficients a and ν, in order to obtain a solution. In particular in [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF] Theorem 2.1, the authors prove the existence of an 'energy solution' under the assumptions (H 0 ) and (H 3 ) (see below).

∃ C > 0 and γ > 1/2 such that:

This condition is not verified in the physical situation (H p ), but only in the approximate situation (H ′ p ).

In Theorem 1 we obtain a W -solution under much simpler conditions which are satisfied by (H p ). This solution is a distributional solution under an additionnal simple assumption (see Corollary 1) which is again satisfied in (H p ).

Note also that in the first part of Theorem 2 we prove that under the assumptions (H 0 ) and (H 2 ) (which are satisfied in (H p ) if a1ν 2 = a 2 ν1), the functions u and k are Hölder continuous. In particular ν(k) ∈ L ∞ which implies that W k = H k , and the notions of H-solution, W -solution, distributional solution and renormalized solution coincide in this case.

In order to conclude this Appendix we give a last existence result. Let (H 4 ) be the following condition:

We have:

Proposition 10 Assume that (H 0 ), (H 1 ) and (H 4 ) hold. Then the W -solution given in Theorem 1 is an 'energy solution' (in the sense of [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF]).

Proof

We have assumed that (H 0 ), (H 1 ) hold and consequently all the results presented in the sections 2, 3 and 4 can be recovered.

Let (u, k) be the W -solution given by Theorem 1. By using [START_REF] Boccardo | Non-linear elliptic or parabolic equations involving measure date[END_REF] we see that the second condition in (R) is satisfied. Nevertheless we cannot directly conclude that ν(k) ∈ H 1 (Ω), but we can obtain a new estimate for the approximating sequence (u n , k n ). More precisely, we have:

In fact, by using the property that

It follows that:

Hence, by using (21) we obtain

Moreover ν n (k n ) = ν(0) on ∂Ω and thus we obtain (63) by using a Poincaré inequality.

Remark

The hypotheses made in Proposition 10 are verified under assumption (H ′ p ). In the hypotheses, we require only very weak growth condition at infinity for ν. For instance (contrarily to the result presented in [START_REF] Gallouet | On a turbulent system with unbounded eddy viscosities[END_REF]) the Proposition 10 works if we have:

Appendix II: Hölder continuity and composition

Let Λ ⊂ R d and α ∈ (0, 1). We recall that the space C 0,α (Λ) of Hölder continuous (with exponent α) functions on Λ is defined by:

More generally, for any integer k, the space C k,α (Λ) is the space of those f ∈ C k (Λ) whose kth derivative is in C 0,α (Λ).

A first elementary result tells that the product of two Hölder continuous functions is an Hölder continuous function. More precisely we have (see relation (4.7) in [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]):

In Section 5 we used a function defined as a composition of two Hölder continuous functions. We needed the following result:

Lemma 12 Let Ω be a compact in R d and α ∈ (0, 1). We consider the following three conditions:

We have:

In this situation we clearly have λ(k) ∈ C 0 (Ω) and

We want to prove that sup

x,x 0 ∈Ω I(x, x 0 ) < ∞.

(63)

• Assume that (A) holds. Then in addition of (62) we have: In this situation we can estimate I(x, x 0 ) as follows:

Hence (63) is again satisfyed.

2. Assume that (C) holds and let µ := λ(k). Clearly µ ∈ C 1 (Ω) and ∇µ = λ ′ (k)∇k. We remark that λ ′ ∈ C 0,α (R) and k ∈ C 1,α (Ω). We can then apply the first point of this lemma to obtain: λ ′ (k) ∈ C 0,α (Ω). Moreover ∇k ∈ (C 0,α (Ω)) d . Hence the product λ ′ (k)∇k is Hölder continuous.