domain: info.info-mo math sde

text:
Jérôme LUQUEL

Advisors : Hélène BARUCQ, Julien DIAZ & Henri CALANDRA Mail : jerome.luquel@inria.fr Equipe-Projet Magique-3D INRIA -Bordeaux Sud-Ouest Laboratoire de Mathématiques et leurs Applications, UPPA

Imaging of complex media with elastic wave equations

Motivation and context

The drilling is expensive, the petroleum industry uses methods able to give images of the area before to do the drilling. From all the famous techniques used by oil companies, the Reverse Time Migration (RTM) is very efficient.

Seismic survey on earth

Seismic survey on sea Acquisition methods Speed model

• The seismic acquisition is the artificial generation and the recording of seismic datas.

• The seismic energy is measured by receivers who register seismic signals as waves.

• The RTM works with velocity models of subsurface and allows to do simulations of seismic acquisitions.

Modelling of the problem

Equations

We use acoustic and elastic wave equations to model the problem

{ ρ(x)∂ t u(x, t) + ∇p(x, t) = 0, 1 µ(x) ∂ t p(x, t) + ∇ • (u(x, t)) = f (x, t) { ρ(x)∂ t u(x, t) = ∇ • σ(x, t), ∂ t σ(x, t) = C(x)ϵ( -→ v ( -→ x , t)) + f (x, t)
where

(x, t) ∈ Ω × [0, T ]

Imaging conditions

We note u the forward field which is the solution of the previously equations where f (x, t) is the source of the experience, and d the backward field which is the solution when f (x, T -t) are the data of the receivers. We correlate the forward and backward wavefield to image heterogeneities of the subsurface of the subsurface.

I(x) = ∫ T 0 u(x, t) • d(x, t)dt
Others formula exist, we can separate the waves P and S :

I P P (x) = ∫ T 0 [∇ • u(x, t)][∇ • d(x, t)]dt, I SS (x) = ∫ T 0 [∇ × u(x, t)] • [∇ × d(x, t)]dt,
The goal is to build new imaging conditions using physical parameters of an elastic media.

Results

This first result is obtained using the first order wave acoustic equation, the velocity model is the marmousimodel described above and we use a discontinuous Galerkin method to solve equations.

Storage strategy

The Griewank algorithm

Several strategies exist to calculate I, the differences are the computation and the storage of u n .

• Compute u n for each n = N -1, ..., 0. This means O(N 2 ) total steps, which is too expensive for big size problems. In other words, as soon as we need u n we launch a new simulation, this means too many calculus. • Record the entire state u n , n = 0, ..., N -1 and access when we need. The cost is only O(N ) steps but we need a big memory. • Record only a given number of propagations, for example all the k iterations and build the field u n from 0 to N -1 by interpolation. The advantage of this approach is to reduce the implementation cost but it products only an approximate solution.

We have to find a compromise between computation time, precision and storage cost.

• Choose checkpoints {n i : i = 0, ..., N C } between 0 and N and give the buffers {b j : j = 1, ..., N B } for the states which must be saved and the states which must be initialized. • During the forward : store u n at the checkpoint {n = n i j : j = 1, ..., N B } in the buffers. We do a first forward et we record 3 checkpoints (0,6 and 11) in buf1, buf2 et buf3 respectively. We begin the backward at 14. At this step, we need I(u 14 , d 15 ), as in buf3 we have stored u 11 , all that we need to do is a forward starting from n = 11 until 14 (we consider in this case u 11 as a Cauchy condition) Now we are at 13. At this step, we need I(u 13 , d 14 ), as in buf3 we have stored u 11 , all that we need to do is a forward starting from n = 11 until 13 (we consider in this case u 11 as a Cauchy condition) And so on, the table opposite summarise our example until the time 0.

-The columns bufk are the checkpoints recorded in the buffer k -The column recomp is the steps of the forward that they will be recomposed -The numbers in bold are Cauchy conditions for the next forward -The column n is the u n used to produce I(u n , d n+1 )

step
4 3 0 1 3 -3 2 0 1 3 2 2 1 0 1 3 -1 0 0 1 3 -0
Adjoint state method and RTM The adjoint state methode is a general method to calculate the gradient of a cost function depending of state variables which are solutions of a direct problem. J. Tromp (http://www.princeton.edu/geosciences/tromp/) and his colleagues define three kernels in this method :

K ρ (x) = ∫ T 0 ρ(x)∂ t s † (x, T -t) • ∂ t s(x, t)dt, K µ (x) = - ∫ T 0 2µ(x)D † (x, T -t) : D(x, t)dt, K κ (x) = - ∫ T 0 κ(x)[∇ • s † (x, T -t)][∇ • s(x, t)]dt. s † is the adjoint wavefield, D = 1 2 [∇s + (∇s) T ] -1 3 (∇ • s)I et D † is his adjoint.
We focus on performances of imaging method when we use the kernel :

K ′ ρ = K ρ + K µ + K κ ,
This kernel is recognized to reduce numerical artefacts in the case of a global inversion of the earth. But, in seismic migration, the backward wavefield speed d is connected to the displacement of the field s by d(t) = ∂ t s(t) and the forward speed wavefield u is connected to the adjoint wavefield by u(t) = -∂ t s † (T -t) In this thesis, we will extend such imaging methods used in seismology to the RTM.

•

  During thre backward : Calculate u n serveral times from the last checkpoint, using the value stored in the buffer b i N B to initialize the evolution. • When the backward step has reached the last checkpoint : n = n i N B , use the buffer b n storing the corresponding value. Example We consider the following example : N = 15, N B = 3, N C = 6. Let I an imaging condition calculing I(u n , d n+1 ) where u n is the forward solution and d n is the backward solution.


Nadim Zgheib
email: nzgheib@gmail.com

Thomas Bonometti
email: thomas.bonometti@imft.fr

S Balachandar

Tim N Colonius

Zgheib

Keywords: Density currents, Gravity currents, Buoyancy-driven flows, Box model, Spectral methods

We present the results from laboratory experiments and fully resolved simulations pertaining to finite-release turbulent density flows with a non-axisymmetric initial shape. First, we demonstrate that the effects of the initial shape influence the current's evolution well into the long-time phase which would corresponds to the inertial self-similar phase in the case of planar or axisymmetric configurations. Then, we identify the physical mechanisms responsible for this dependence and propose a new model capable of capturing the dynamics of such releases. Finally, we show that this dependence on the initial configuration is robust for various types of gravity currents over a wide range of parameters such as Reynolds number, density ratio, and aspect ratio.

Fig. 1 Field observations of a non-circular gravity current. Aerial views of (top) a building of square cross section and (bottom) a T-shaped building being demolished. The structure of the debris cloud is highly dependent on the building's shape. Observe that material very close to the center of the building moves farther out than material located at the building's extremities. The dynamics of the cloud is reasonably well captured by our proposed model, where the front location of the cloud is marked at equal instants of time as blue lines (color figure online) [START_REF] Faillettaz | Two-threshold model for scaling laws of noninteracting snow avalanches[END_REF][START_REF] Hopfinger | Snow avalanche motion and related phenomena[END_REF] and storms [START_REF] Hall | Wind shear observations in thunderstorm density currents[END_REF]. These flows are driven by a difference in density either stemming from temperature, salinity or suspended sediments.

Predicting the short-time as well as long-time evolution of these material fronts is of crucial interest [START_REF] Simpson | Gravity currents in the laboratory, atmosphere and oceans[END_REF][START_REF] Huppert | Gravity currents: a personal perspective[END_REF]. Most real gravity currents are generated by a sudden or continuous release from an arbitrary, often non-circular, source. However, nearly all the prior theoretical, experimental and numerical investigations are based on planar or axisymmetric configurations [START_REF] Von Karman | The engineer grapples with nonlinear problems[END_REF][START_REF] Benjamin | Density currents and related phenomena[END_REF][START_REF] Hallworth | Axisymmetric gravity currents in a rotating system: experimental and numerical investigations[END_REF][START_REF] Huq | The role of aspect ratio on entrainment rates of instantaneous, axisymmetric finite volume releases of density fluid[END_REF]. The underlying (implicit) assumption is that, after a relatively short transitional phase, the material front becomes either planar or axisymmetric. Here, we present results from laboratory and numerical experiments that were performed with non-axisymmetric finite initial releases. The most striking feature was that the effect of the initial non-circular shape of the release persists for the whole duration of the observation. This is a unique, but robust, behavior of propagating material fronts, which is quite distinct from propagating informational fronts, such as sound waves and shock fronts, which are well known to quickly become independent of the source shape.

In the present paper, we identify the physical mechanisms responsible for this peculiar behavior and present a novel model which accounts for the shape of the initial release and predicts the non-axisymmetric propagation of the front of the gravity current. This model when applied to the problem of building demolition captures the time evolution of the observed debris cloud (Fig. 1).

Non-circular spreading of density currents

When a patch of heavy (resp. light) fluid spreads in a lighter (resp. heavier) ambient, it generally goes through successive stages [START_REF] Huppert | The slumping of gravity currents[END_REF][START_REF] Cantero | On the front velocity of gravity currents[END_REF], namely an acceleration phase at the end of which the current's front velocity reaches its maximum value, a slumping phase with constant or nearly constant front velocity, and eventually a phase denoted here as long-time phase where the buoyancy driving force is balanced either by inertia, viscosity or surface tension, respectively. Note that the long-time phase is often referred to as self-similar phase, because it is possible to find an exact long-time similarity solution of the simplified equations of motion (e.g., the shallowwater equations) in some specific configurations as in the planar or axisymmetric case. When a self-similar solution exists, the time evolution of the front height h N and front velocity u N usually follows the scaling law h N ∼ t α and u N ∼ t β with α and β being some constants. For instance, in the case of axisymmetric gravity currents, the one-layer shallow-water equations give α = -1 and β = -1/2. The reader is referred to the work [START_REF] Zemach | Gravity currents in non-rectangular cross-section channels: analytical and numerical solutions of the one-layer shallow-water model for high-Reynolds-number propagation[END_REF][START_REF] Ungarish | An Introduction to Gravity Currents and Intrusions[END_REF] for more details about the existence of a similarity solution in various configurations.

While the slumping and self-similar phases have been extensively studied in the past [START_REF] Hoult | Oil spreading in the sea[END_REF][START_REF] Benjamin | Density currents and related phenomena[END_REF][START_REF] Huppert | The slumping of gravity currents[END_REF][START_REF] Cantero | On the front velocity of gravity currents[END_REF][START_REF] Rottman | Density currents produced by instantaneous releases of a heavy fluid in a rectangular channel[END_REF][START_REF] Klemp | On the dynamics of density currents in a channel[END_REF], the initial acceleration phase has received very little attention [START_REF] Cantero | On the front velocity of gravity currents[END_REF][START_REF] Cantero | Turbulent structures in planar gravity currents and their influence on the flow dynamics[END_REF] as initial conditions were believed to be unimportant in the prediction of the long-term dynamics. Here, we argue, by means of experiments and fully resolved simulations, that a non-axisymmetric finite-volume release does not reach an axisymmetric shape nor preserves its initial shape for a significant time. In the following, we define the "long-time" regime as the regime for which the current has entered a phase for which the local front height and speed obey a scaling law of the form h N ∼ t α and u N ∼ t β , as found for the self-similar (inertial) phase.

The laboratory experiments consist of swiftly releasing a column of dense salty water (ρ c ≈ 1,100 kg/m 3 ) of height h 0 into a lighter ambient fluid of height H inside a 1.2 × 1.2× 0.4-m tank, thanks to a hollow cylinder of equivalent radius R 0 . The ambient fluid is either clear water (ρ a ≈ 1,000 kg/m 3 ) or air (ρ a ≈ 1.2 kg/m 3 ). Unless stated otherwise, we set h 0 = H (full-depth release). Three different cross-sectional shapes are considered: (a) a circular section, (b) a rounded rectangular section, i.e., a rectangle where the two shorter edges are replaced by semicircles, and (c) a plus-shape section for which concave corners are present (see e.g., Fig. 2). Fluorescent dye is added to the fluid inside the cylinder. The front location and the current's height are measured thanks to a mirror placed beneath the tank, which allows for a plan (bottom) view of the front evolution, while the side view of the current provides information about the height's evolution.

Several experiments have been performed for a wide range of Reynolds number Re = HU/ν(U = [g(ρ cρ a )H/ρ a ] 1/2 being the velocity scale, g the gravitational acceleration, and ν the kinematic viscosity of water), initial height aspect ratio H/R 0 , and density ratio ρ c /ρ a . As seen in Table 1, the initial aspect ratio based on the local distance from the center of mass is in the range 0.5-8.5, covering both configurations of O(10 -1 ) and O(10 1 ) aspect ratio. In all cases, the flow was fully turbulent. The commonly accepted distance of propagation for which transition from the slumping phase to the inertial self-similar phase occurs is about 2 (resp. 5-9) in the cylindrical (resp. planar) configuration [START_REF] Rottman | Density currents produced by instantaneous releases of a heavy fluid in a rectangular channel[END_REF][START_REF] Cantero | High-resolution simulations of cylindrical density currents[END_REF]. Here, the distance of propagation was in the range 7-25; hence, the current is likely to enter a regime which resembles the self-similar phase in all cases so the "long-time" regime is reached, as it will be confirmed later (Figs. 3,4).

H/R max H/R 0 H/R min 0.5 1-2-4 8.5 Density ratio ρ c /ρ a 0.93-1.1-10 3 Distance of propagation L/R max L/R 0 L/R min 7 12 25
h 0 (H ) is the initial height of the current (ambient), R 0 , R min , R max are the initial equivalent, minimum, maximum radius of the cross-sectional area of the cylinders, L is the maximum distance of propagation

As mentioned above, the behavior of material fronts is in stark contrast to propagation of information (or wave) fronts. For example, an outward propagating sound wave front quickly becomes spherical irrespective of the shape of the source [START_REF] Lighthill | Waves in Fluids[END_REF] and the front radius is linearly related to time through sound speed as r = c 0 t (note that sound intensity distribution depends on the source detail). Similarly, a blast wave from a finite source quickly becomes spherically independent of the source shape. Here, again the blast radius can be expressed by the power law [START_REF] Sachdev | Shock Waves and Explosions[END_REF] r ∼ t 2/5 . Although the blast front velocity decreases with time, it remains the same along the entire blast front. Clearly, the propagation of a material front, e.g., gravity currents, cannot be modeled as that of an information front.

The temporal evolution of some axisymmetric and non-axisymmetric gravity currents is presented in Figs. 2, 3 and 4. When the release is non-circular, the current's evolution depends on the initial shape within the time of the observation. For instance, in the case of the gravity current of initial rounded-rectangle cross section, the local front velocity and height at the tips of the initial major (slow) and minor (fast) axes (see the marks S and F in Fig. 7) show large differences in magnitude early after material redistribution (see Figs. 3,4, t > 2). Over time, the magnitude of the difference somewhat decreases, but persists up to the end of the experiment/simulation, indicating that the current has not become axisymmetric nor reverted back to its original shape within the limited time of observation. For instance, the patch of heavy fluid of initial rounded-rectangle shape is observed to flip axes, while that of initial plus shape turns into a square.

Initially, the height of the patch is uniform across the surface, and hence, the pressure force is the same along the current's front. As a consequence, just after release, the current rapidly accelerates outward with increasing speed, but initially the current's speed u N and height h N are independent of the initial shape and are uniform around the circumference of the front (see Figs. 3,4 for t <1). Upon the release of the heavy fluid inward propagating perturbations initiate at the front, which reflect back at the center of the patch and eventually catch up to the front [START_REF] Ungarish | An Introduction to Gravity Currents and Intrusions[END_REF]. During the acceleration phase, the heavy material is redistributed within the patch following the direction perpendicular to the initial front. At the end of the redistribution phase, the height h N is not uniform along the front of the current (Fig. 3, t ≈ 2). Since the front velocity u N scales as √ h N (see e.g., Ref. [START_REF] Huppert | Gravity currents: a personal perspective[END_REF][START_REF] Huppert | The slumping of gravity currents[END_REF]), the speed of propagation is non-uniform along the current's front. This non-uniform circumferential distribution of the heavy fluid within the current, as dictated by the initial shape, remains fixed over time leading to a non-axisymmetric spreading of the current during the slumping and the subsequent spreading phase of the current.

For comparison, the time-dependent solution of the axisymmetric one-layer shallow-water equations (equations 6.24-6.25 in Ref. [START_REF] Ungarish | An Introduction to Gravity Currents and Intrusions[END_REF]) obtained with a finite-difference method similar to that described in appendix 2 of [START_REF] Ungarish | An Introduction to Gravity Currents and Intrusions[END_REF] is plotted in Figs. 3b and4b together with the fully resolved simulation and experimental results. Figure 3b shows that the temporal evolution of the local front height of the gravity current of initial rounded-rectangle cross section roughly follows a slope -1 in log-log representation at times larger than t ≈ 5 approximately. The same trend is observed for the time-dependent solution of the one-layer axisymmetric shallow-water equations at a somewhat earlier time t ≈ 3. Similarly, the local front speed in Fig. 4b is observed to roughly follow a slope -1/2, as clearly seen for the experimental local "fast" front speed. Note that some fluctuations are present in the evolution of the local "slow" front speed, making the comparison more difficult. Comparing the temporal evolution of the local front speed of the non-axisymmetric gravity currents with the time-dependent solution of the one-layer axisymmetric shallow-water equations indicates that the non-axisymmetric gravity currents have entered a long-time phase which resembles the self-similar regime of the axisymmetric configuration, in that the local height and speed roughly follow a law of the type h N ∼ t α and u N ∼ t β , with α and β being some constants. Note, however, that determining the precise value of α and β, and the time at which the long-time regime starts is difficult with the present set of experimental and numerical data because of the significant velocity fluctuations observed in Fig. 4b and the somewhat limited range of parameters investigated here. Larger-size experiments and/or simulations would help to clarify this point.

The dependence of material front propagation on initial condition of release is robust in the sense that the peculiar behavior observed in Fig. 2 is not restricted to the presently discussed configurations (roundedrectangle and plus-shaped cross sections; Boussinesq currents of density ratio close to unity). We have conducted many more laboratory and numerical experiments of different non-axisymmetric geometries, density ratios ρ c /ρ a (including Boussinesq currents ρ c /ρ a ≈ 1, heavy currents of dam-break flow type ρ c /ρ a = O( 103 ), and light currents ρ c /ρ a < 1), different aspect ratio releases (radius to height ratio ranging between 0.5 and 7), and Reynolds numbers (see Table 1). We found that provided the Reynolds number is large enough, say Re ≥ O( 103 ), the initial shape of the current always influences the current's continued propa- gation as well as its preferred asymptotic non-axisymmetric shape for a significant amount of time. Finally, we performed simulations of both full-depth and partial-depth releases (Fig. 5) and found a qualitative similar dynamics. This indicates that the presently observed behavior is independent of the initial depth ratio. Note that we performed simulations for only two values of h 0 /H . Other experiments or simulations for a wider range of depth ratios would help in quantitatively determining the influence of this parameter on the dynamics of non-cylindrical gravity currents.

A new model for the prediction of the propagation of non-circular density flows

This finding has several theoretical implications. The classical box model, which is classically used for predicting the evolution of gravity currents [START_REF] Huppert | The slumping of gravity currents[END_REF], despite its simplicity, is able to reproduce the dynamics of axisymmetric and planar releases. However, straightforward application of the box model fails for non-axisymmetric releases (Fig. 7). According to this model, the height remains uniform along the entire spreading patch, so the speed of propagation remains uniform along the current's front during all the phases of spreading (recall that the local speed of propagation evolves as the square root of the local height of the current). Using the classical box model, an initially non-axisymmetric current inevitably becomes axisymmetric. Similarly, theories based on slumping and self-similar phases also fail to predict the sensitive dependence on the initial shape and the preferential propagation of non-axisymmetric gravity currents for the same reasons.

Fully resolved simulations support the experimental findings (Figs. [START_REF] Gröbelbauer | The propagation of intrusion fronts of high density ratio[END_REF][START_REF] Doyle | Power laws, highly optimized tolerance, and generalized source coding[END_REF][START_REF] Meiburg | Turbidity currents and their deposits[END_REF][START_REF] Faillettaz | Two-threshold model for scaling laws of noninteracting snow avalanches[END_REF]. The simulations are performed using a spectral code [START_REF] Cantero | High-resolution simulations of cylindrical density currents[END_REF] to solve the Navier-Stokes equations using the Boussinesq approximation ρ c /ρ a ≈ 1. The numerical domain consists of a Cartesian parallelepiped (length × width × height = 15 × 15 × 1), with a spatial resolution of 880 × 880 × 179 (140 million degrees of freedom). Boundary conditions are no-slip at the bottom wall, free-slip at the top, and periodic at the sidewalls. The propagation of the current front is visualized via iso-contours of the vertically averaged dimensionless density field ρ = (ρ *ρ a )/(ρ cρ a ) = 0.001 (ρ * is the dimensional local density). Quantitative agreement is found between experiments and fully resolved simulations. The undulations seen in Fig. 6 for both simulations and experiments at later times are due to the lobe and cleft instability [START_REF] Simpson | Effect of the lower boundary on the head of a gravity current[END_REF].

Even though such simulations are able to reproduce the peculiar dynamics of non-axisymmetric gravity currents, they are unlikely to be used for rapid prediction, as needed in operational models especially those dealing with the high-Reynolds numbers gravity currents.

Here, we propose an extended box model based on partitioning of the initial release using geometric rays that are perpendicular to the front (Fig. 7). Once the various subvolumes are obtained, the local fronts are advanced normal to themselves as in the box model. This initial partitioning is the key aspect of the present model, since it allows for non-uniform height and speed along the patch's advancing front, during all the phases of spreading. This allows the model to capture the non-axisymmetric propagation of the front. To be explicit, the formulation of the extended box model makes use of a Benjamin-type boundary condition (1) relating the outward normal front velocity u to the front height h, kinematic relations (2-3) for the advancement of the

x y -3 -2 -1 0 1 2 3 -4 -2 0 2 4 x -3 -2 -1 0 1 2 3 -4 -2 0 2 4 (a) (b)
Fig. 6 Temporal evolution of a non-axisymmetric material front (fully resolved simulation vs. experiment). The initial crosssectional geometry is a rounded rectangle. a Numerical simulations with Re = 8,950. Time separation between contours is t = 0.35 and the final time is t f = 12.6; b laboratory experiments. Here, Re = 28,000, time separation between contours is t = 0.26, and the final time is t f = 12.73 front position (x, y) and horizontal area per arc length σ , respectively, and mass conservation (4). This results in a system of coupled nonlinear PDEs for the unknowns u, x, y, σ and h, viz

u = Fr √ h (1) ∂{x, y} ∂t = u {∂ y/∂s, -∂ x/∂s} [(∂ x/∂s) 2 + (∂ y/∂s) 2 ] 1/2
(2)

∂σ ∂t = u (3) ∂σ h ∂t = 0 (4)
where, here, Fr is the Huppert-Simpson Froude number1 (see [START_REF] Huppert | The slumping of gravity currents[END_REF]). Note that since the flow is incompressible and entrainment is neglected, the area per arc length σ and the current height h is such that the total volume V of the current is given by

V = s σ hds (5)
The independent variables s and t denote the curvilinear coordinate along the front and time, respectively. The solution to (1)-( 4) is far easier and faster than the direct numerical simulations displayed in Fig. 6a. As shown in Figs. 1 and7, the solution of the extended box model is capable of capturing the propagation of currents with arbitrary initial forms.

At this point, we have shown that (i) fully resolved simulations and (ii) the proposed extension of the box model were able to reproduce the observed propagation of non-axisymmetric gravity currents contrary to the classical box model. It is of major interest, however, to assess the capability of approaches based on the shallow-water equations to reproduce such a dynamics. This would require the development of a numerical approach for solving the two-dimensional one-layer (or better yet two-layer) shallow-water equations as done for example in [START_REF] Rocca | Experimental and numerical simulation of three-dimensional gravity currents on smooth and rough bottom[END_REF]. The development of such an approach is, however, a non-trivial undertaking, which is beyond the scope of the present work. Comparing the present results with a one-and/or two-layer shallowwater approach would allow clarifying the capabilities of shallow-water approaches within the hierarchy of available models for the description of non-axisymmetric gravity currents.

Summary and discussion

We have presented results from laboratory experiments and numerical simulations of the propagation of turbulent material fronts stemming from non-axisymmetric finite initial releases. It is found that the effect of the initial non-circular shape of the release persists for the whole duration of the observation. The duration of the experiments was such that the gravity currents have crossed a distance of 12-25 times the initial radius, if one considers the equivalent or minimum radius, respectively. This allowed to cover the acceleration phase, the (quasi-)slumping phase and a regime for which the dynamics resembles the self-similar inertial phase predicted by the one-layer axisymmetric shallow-water equations. Finally, we presented a novel model, extending the classical box model, which accounts for the shape of the initial release and predicts the non-axisymmetric propagation of the front of the gravity current.

It is important to note that the local speed of propagation of a material front generated by the release of a patch of arbitrary shape can vary significantly, thus leading to local "fast fronts" and "slow fronts." In Fig. 4, the fastest front is twice as fast as the slowest front during all the observed phases of spreading. Such long-lasting speed variations between the different sections of the front may result in dramatically different front locations that depend on the shape of the initial release. In the context of massive oil spills such as that the Deepwater Horizon oil spill stemming from the explosion of a sea-floor oil gusher in the Gulf of Mexico in April 20, 2010, the flow is likely to remain turbulent for longer times contrary to more moderate oil spills generated by tankers running aground, for which the gravity current is likely to be dominated by viscous and subsequently capillary effects after an hour approximately [START_REF] Hoult | Oil spreading in the sea[END_REF]. In the former case, the error made in the estimation of the propagation of the oil spill front, which depends on the nature of initial release, unavoidably decreases the predictive capability of the precise location of impact along the coast. Clearly, many other factors, such as currents, cross flows, bottom topology, further influence the propagation of non-axisymmetric gravity currents. In light of the present findings, suitable additions to the extended box model described here can help improve the prediction of such gravity currents of arbitrary shapes.

Finally, it must be stressed that the present experiments/simulations were done for a limited range of parameters and on a limited spatial domain. The initial shape of the non-cylindrical gravity currents was found to influence the dynamics during the whole but limited duration of observation. For very large domains, however, non-cylindrical finite-release gravity currents are likely to enter, after some time, a regime where viscous effects are predominant. In that case, viscous diffusion of momentum may homogenize the front height and velocity so that the viscous current may become axisymmetric. More experiments on larger domains are needed to clarify whether the presently observed behavior holds for longer distances of propagation and corresponding times than those attained here.

Fig. 2

 2 Fig. 2 Temporal evolution of the experimental collapse of a column of salt water with different cross sections at the center of a tank containing fresh water (H/R 0 = 2, Re = 2.8 × 10 4 , ρ c /ρ a = 1.1). The schemes on the left depict the three-dimensional initial shape of the heavy material in the experiment, namely cylinders of a circular, b rounded-rectangle and c plus-shaped cross sections. Gravity is oriented opposite to z axis. Time is scaled by T = H /U (see text for definition)

Fig. 3 Fig. 4

 34 Fig. 3 Temporal evolution of the local front height in the vertical mid-plane of x direction (dashed line) and y direction (solid line) of the current of initial rounded-rectangle cross section. a Results obtained from the fully resolved simulation with H/R 0 = 2 and Re = 8,950. b Same as (a) in log-log representation. For comparison, the time-dependent solution of the axisymmetric one-layer shallow-water equations using Huppert and Simpson's [15] front condition in the full-depth configuration is shown in (b) (dash-dot line). The dotted lines indicate a slope of -1, as predicted by the self-similar solution of the axisymmetric one-layer shallow-water equations (see e.g., Ref. [23, p. 122])

Fig. 5

 5 Fig. 5 Temporal evolution of the local front's height (a) and speed (b) in the vertical mid-plane of x direction (dashed lines) and y direction (solid lines) of the current of initial rounded-rectangle cross section for two initial depth ratios, namely a partial-depth release h 0 = H/2 (blue lines with symbols) and a full-depth release h 0 = H (black lines without symbols). These results are obtained from fully resolved simulations with H/R 0 = 2 and Re = 8,950 (color figure online)

Fig. 7

 7 Fig. 7 Experiments versus box model/extended box model. Front location with a time separation between contours of t = 1.3, and a final time of t f = 13. In the classical box model, as the current propagates, the height is intrinsically averaged over the entire patch of fluid making the propagation inevitably axisymmetric. In the extended box model, the volume of release is initially divided into multiple subvolumes. The size of each subvolume depends on the inwardly propagating geometric rays starting from and perpendicular to the current's front. The extended box model is in quantitative agreement with experiments, contrary to the classical box model. The marks S and F refer to the tips of the initial major (slow) and minor (fast) axes from which the front velocity and height are computed in Figs. 3 and 4

Table 1

 1 Parameters used in the experiments and/or simulations

	Reynolds number	Re
	Initial depth ratio	[7 × 10 3 , 10 6 ] H/ h 0
		1-2
	Initial aspect ratio

Note that another model of Froude number function could be used without losing generality, provided this function is valid for the whole range of height ratio of nose to ambient h/H , this ratio being here in the range 0 ≤ h/H ≤ 1.

Acknowledgments

This work was supported by the National Science Foundation under Partnership for International Research and Education (PIRE) in Multiphase Flows at The University of Florida and by the French Embassy through the Chateaubriand Fellowship at the Institut de Mecanique des Fluides de Toulouse (IMFT), where the experiments were conducted.


Olivier Faugeras

James Maclaurin
email: james.maclaurin@inria.fr

A Large Deviation Principle and an Expression of the Rate Function for a Discrete Stationary Gaussian Process

Keywords: Stationary Gaussian process, Large Deviation Principle, Spectral Density

We prove a Large Deviation Principle for a stationary Gaussian process over | b , indexed by ) d (for some positive integers d and b), with positive definite spectral density and provide an expression of the corresponding rate function in terms of the mean of the process and its spectral density. This result is useful in applications where such an expression is needed.

Introduction

In this paper we prove a Large Deviation Principle (LDP) for a spatially stationary, indexed by ) d , Gaussian Process over | b , and obtain an expression for the rate function.

Our work in mathematical neuroscience involves the search for asymptotic descriptions of large ensembles of neurons [START_REF] Baladron | Mean-field description and propagation of chaos in networks of Hodgkin-Huxley and Fitzhugh-Nagumo neurons[END_REF][START_REF] Faugeras | Asymptotic description of neural networks with correlated synaptic weights[END_REF][START_REF] Faugeras | A constructive mean field analysis of multi population neural networks with random synaptic weights and stochastic inputs[END_REF]. Since there are many sources of noise in the brains of mammals [START_REF] Rolls | The noisy brain: stochastic dynamics as a principle of brain function[END_REF], the mathematician interested in modelling certain aspects of brain functioning is often led to consider spatial Gaussian processes that (s)he uses to model these noise sources. This motivates us to use Large Deviations techniques. Being also interested in formulating predictions that can be experimentally tested in neuroscience laboratories, we strive to obtain analytical results i.e. effective results from which for example numerical simulations can be developed. This is why we determine a more tractable expression for the rate function in this article.

Our result concerns the Large Deviations of ergodic phenomena, the literature of which we now briefly survey. Donsker and Varadhan obtained a Large Deviations estimate for the law governing the empirical process generated by a Markov Process [START_REF] Donsker | Asymptotic evaluation of certain markov process expectations for large time, iv[END_REF]. They then determined a Large Deviations Principle for a )-indexed Stationary Gaussian Process, obtaining a particularly elegant expression for the rate function using spectral theory. [START_REF] Bryc | Large deviations and strong mixing[END_REF][START_REF] Chiyonobu | The large deviation principle for hypermixing processes[END_REF][START_REF] Deuschel | Microcanonical distributions for lattice gases[END_REF] obtain a Large Deviations estimate for the empirical measure generated by ergodic processes satisfying certain mixing conditions. [START_REF] Baxter | An approximation condition for large deviations and some applications[END_REF] obtain a variety of results for the Large Deviations of ergodic phenomena, including one for the Large Deviations of )-indexed | b -valued stationary Gaussian processes. [START_REF] Steinberg | On tests for normality[END_REF] have proved an LDP for a stationary ) d -indexed Gaussian Process over |, and [START_REF] Bryc | On large deviations of empirical measures for stationary gaussian processes[END_REF] obtain an LDP for an |-indexed, |-valued stationary Gaussian process.

In the work we are developing [START_REF] Faugeras | Asymptotic description of stochastic neural networks[END_REF], we need Large Deviation results for spatially ergodic Ornstein-Uhlenbeck processes. This requires Theorem 1 of this paper.

In the first section we make some preliminary definitions and state the chief theorem, Theorem 1, for zero mean processes. In the second section we prove the theorem. In appendix A, we state and prove several identities involving the relative entropy which are necessary for the proof of Theorem 1. In appendix B we prove a general result for the Large Deviations of exponential approximations. We prove Corollary 2 which extends the result in Theorem 1 to non-zero mean processes in appendix C.

Preliminary Definitions

For some topological space Ω equipped with its Borelian σ-algebra B(Ω), we denote the set of all probability measures by M(Ω). We equip M(Ω) with the topology of weak convergence. Our process is indexed by ) d . For j ∈ ) d , we write j = (j(1), . . . , j(d)). For some positive integer n > 0, we let

V n = {j ∈ ) d : |j(δ)| ≤ n for all 1 ≤ δ ≤ d}. Let T = | b ,
for some positive integer b. We equip T with the Euclidean topology, T ) d with the cylindrical topology, and denote the Borelian σ-algebra generated by this topology by B(T ) d ). For some µ ∈ M(T ) d ) governing a process (X j ) j∈) d , we let µ Vn ∈ M(T Vn ) denote the marginal governing (X j ) j∈Vn . For some j ∈ ) d , let the shift operator S j : T ) d → T ) d be S(ω) k = ω j+k . We let M s (T ) d ) be the set of all stationary probability measures µ on (T

) d , B(T ) d )) such that for all j ∈ ) d , µ • (S j ) -1 = µ. We use Herglotz' Theorem to characterise the law Q ∈ M s (T ) d )
governing a stationary process (X j ) in the following manner. We assume that E[

X j ] = 0 and E[X 0 † X k ] = 1 (2π) d [-π,π] d exp(i k, θ ) G(θ)dθ. (1)
Our convention throughout this paper is to denote the transpose of X as † X and the spectral density with a tilde.

•, • is the standard inner product on | d . Here G(θ) is a continuous function [-π, π] d → C b×b
, where we consider [-π, π] d to have the topology of a torus. In addition G(-θ) = † G(θ) = Ḡ (x indicates the complex conjugate of x). We assume that for all θ ∈ [-π, π] d , det G(θ) > Gmin for some Gmin > 0, from which it follows that for each θ, G(θ) is Hermitian positive definite. If x ∈ C b , then † xX j is a stationary sequence with spectral density † x Gx. We employ the operator norm over C b×b . Let p n : T ) d → T ) d be such that p n (ω) k = ω k mod Vn . Here, and throughout the paper, we take k mod V n to be the element l ∈ V n such that for all 1 ≤ γ ≤ d, l(γ) = k(γ) mod (2n + 1). Define the process-level empirical measure μn :

T ) d → M s T ) d as μn (ω) = 1 (2n + 1) d k∈Vn δ S k pn(ω) . (2)
Let Π n ∈ M M s (T ) d ) be the image law of Q under μn . We note that in this definition we need not have chosen V n to have an odd side length (2n+1): this choice is for notational convenience, and these results could easily be reproduced in the case that V n has side length n.

In the context of mathematical neuroscience, (X j ) could correspond to a model of interacting neurons on a lattice (d = 1, 2 or 3), as in [START_REF] Faugeras | Asymptotic description of stochastic neural networks[END_REF][START_REF] Faugeras | Large deviations of a spatially-ergodic neural network with learning[END_REF]. We note that the Large Deviation Principle of this paper may be used to obtain an LDP for other processes using standard methods, such as the Contraction Principle or Lemma 13. Definition 1. Let (Ω, H) be a measurable space, and µ, ν probability measures.

I (2) (µ||ν) = sup f ∈B {E µ [f ] -log E ν [exp(f )]} ,
where B is the set of all bounded measurable functions. If Ω is Polish and H = B(Ω), then we only need to take the supremum over the set of all continuous bounded functions.

Let (Y j ) be a stationary Gaussian Process on

T such that E[Y j ] = 0, E[Y j † Y k ] = 0 and E[Y j † Y j ] = Id b . Each Y j
is governed by a law P , and we write the governing law in M s (T ) d ) as P ) d . It is clear that the governing law over V n may be written as P ⊗Vn (that is the product measure of P , indexed over V n ). Definition 2. Let E 2 be the subset of M s (T ) d ) defined by

E 2 = {µ ∈ M s (T ) d ) | E µ [ ω 0 2 ] < ∞}.
We define the process-level entropy to be, for µ ∈ M s (T ) d )

I (3) (µ) = lim n→∞ 1 (2n + 1) d I (2) µ Vn ||P ⊗Vn . It is a consequence of Lemma 11 that if µ / ∈ E 2 , then I (3) (µ) = ∞.
For further discussion of this rate function, and a proof that I (3) is well-defined, see [START_REF] Deuschel | Microcanonical distributions for lattice gases[END_REF]. If furthermore the set {x : I(x) ≤ α} is compact for all α ≥ 0, we say that I is a good rate function.

Definition 4. For µ ∈ E 2 , we denote the C b × C b -valued spectral measure on ([-π, π] d , B([-π, π] d
)) (which exists due to Herglotz's Theorem) by μ. We have

1 (2π) d [-π,π] d exp(i k, θ )dμ(θ) = E µ ω 0 † ω k . For θ ∈ [-π, π] d , let H(θ) = G(θ) -1
2 be the Hermitian positive definite square root of G(θ) -1 and

H(θ) = j∈) d H j exp (-i j, θ ) . (3)
The b × b matrices H j are the coefficients of the absolutely convergent Fourier series (due to Wiener's Theorem) of G-1/2 . Define β : T ) d → T ) d as follows

(β(ω)) k = j∈) d † H j ω k-j .
The Theorem below is the chief result of this paper.

Theorem 1. (Π n ) satisfies a strong LDP with good rate function

I (3) (µ • β -1
). Here

I (3) µ • β -1 =    I (3) (µ) -Γ(µ) if µ ∈ E 2 ∞ otherwise, (4)
where

Γ(µ) =    Γ 1 + Γ 2 (µ) if µ ∈ E 2 0 otherwise.
Here

Γ 1 = - 1 2(2π) d [-π,π] d log det G(θ)dθ Γ 2 (µ) = 1 2(2π) d [-π,π] d tr Id b -G(θ) -1 dμ(θ) .
Finally, the rate function uniquely vanishes at µ = Q.

Corollary 2. Suppose that the underlying process Q is defined as previously, except with mean E Q [ω j ] = c for all j ∈ ) d and some constant c ∈ | b . If we denote the image law of the empirical measure by Π n c , then (Π n c ) satisfies a strong LDP with good rate function (for µ ∈ E 2 )

I (3) (µ) -Γ(µ) + 1 2 † c G(0) -1 c - † c G(0) -1 m µ . (5)
Here

m µ = E µ(ω) [ω j ] for all j ∈ ) d . If µ / ∈ E 2 ,
then the rate function is infinite. The rate function has a unique minimum, i.e.

I (3) (µ) -Γ(µ) + 1 2 † c G(0) -1 c - † c G(0) -1 m µ = 0 if and only if µ = Q.
We prove this in appendix C.

Proof of Theorem 1

In this proof we essentially adapt the methods of [START_REF] Baxter | An approximation condition for large deviations and some applications[END_REF][START_REF] Donsker | Large deviations for stationary Gaussian processes[END_REF]. We introduce the following metric over T ) d . For j ∈ ) d , let

λ j = 3 -d d δ=1 2 -|j(δ)| . Define the metric d λ as follows, d λ (x, y) = j∈) d λ j min x j -y j , 1 ,
where the above is the Euclidean norm. Let d λ,M be the induced Prohorov metric over M(T ) d ). These metrics are compatible with their respective topologies.

For θ ∈ [-π, π] d , let F (θ) = G(θ) 1/2 be
the Hermitian positive square root of G and

F (θ) = j∈) d F j exp (-i j, θ ) . (6)
The b × b matrices F j are the coefficients of the absolutely convergent Fourier series of the positive square root. Define τ :

T ) d → T ) d and τ (M ) : T ) d → T ) d as follows, (τ (ω)) k = j∈) d † F j ω k-j , τ (n) (ω) k = j∈Vn † F j ω k-j d δ=1 1 - |j(δ)| 2n + 1 . (7)
We note that τ = β -1 (on a suitable domain where the series are convergent) and that τ (n) is a continuous map, but τ is not continuous (in general).

We note (using Lemma 6) that P ) d • τ -1 has spectral density G. We define

F(n) (θ) = j∈Vn F j exp (-i j, θ ) d δ=1 1 - |j(δ)| 2n + 1 . (8)
We write ε

(n) = sup θ∈[-π,π] d F (θ) -F(n) (θ) 2 . By Fejer's Theorem, ε(n) → 0 as n → ∞. We define G(n) (θ) = F(n) (θ) 2 , noting that this is the spectral density of P ) d •τ -1 (n) . Let Γ (n) (µ) = Γ 1,(n) + Γ 2,(n) (µ), where for µ ∈ E 2 , Γ 1,(n) = - 1 2(2π) d [-π,π] d log det G(n) (θ)dθ, (9)
Γ 2,(n) (µ) = 1 2(2π) d [-π,π] d tr Id b -G-1 (n) (θ) dμ(θ) . ( 10
) If µ / ∈ E 2 , let Γ 1,(n) = Γ 2,(n) (µ) = 0. Let R n ∈ M(M s (T ) d ))
be the law governing μn (Y ), where we recall that the stationary process (Y j ) is defined just below definition 1.

Let Π n (m) ∈ M(M s (T ) d )) be the law governing μn (τ (m) (Y )).
Lemma 3. (R n ) satisfies a Large Deviation Principle with good rate function

I (3) (µ). If µ / ∈ E 2 , then I (3) (µ) = ∞.
Proof. The first statement is proved in [START_REF] Deuschel | Microcanonical distributions for lattice gases[END_REF]. The last statement follows from Lemma 11 below.

Lemma 4. (Π n (m) ) satisfies a strong LDP with good rate function given by, for µ ∈ E 2

I (3) (m) (µ) = inf ν∈E 2 :µ=ν•τ -1 (m) I (3) (ν) = I (3) (µ) -Γ (m) (µ). ( 11
)
If µ / ∈ E 2 , then I (3) (m) (µ) = ∞.
Proof. The sequence of laws governing μn (Y ) • τ -1 (m) (as n → ∞, with m fixed) satisfies a strong LDP with good rate function as a consequence of an application of the contraction principle to Lemma 3 (since τ (m) is continuous). Now, through the same reasoning as in [10, Lemma 2.1,Theorem 2.3], it follows from this that (Π n (m) ) satisfies a strong LDP with the same rate function as that of μn (Y ) • τ -1 (m) . The last identification in [START_REF] Faugeras | Asymptotic description of neural networks with correlated synaptic weights[END_REF] follows from Lemmas 6 and 9 in Appendix A. We only need to take the infimum over E 2 because by Lemma 3,

I (3) (ν) is infinite for ν / ∈ E 2 . Lemma 5. If 0 < b < 1 2ε(m) then for all n > 0 1 (2n + 1) d log E P ) d exp b k∈Vn τ (m) (ω) k -τ (ω) k 2 ≤ - 1 2 log (1 -2bε(m)) .
The proof is almost identical to that in [10, Lemma 2.4]. We are now ready to prove Theorem 1.

Proof. We apply Lemma 13 in the Appendix to the above result. We substitute Y (m) = τ (m) (ω) and W = τ (ω). Taking m → ∞ in the equation in Lemma 5, we find (38) is satisfied if we stipulate that κ = 0. After noting the LDP governing Π n (m) in ( 11), we may thus conclude that (Π n ) satisfies a strong LDP with good rate function

lim δ→0 lim m→∞ inf γ∈B δ (µ) I (3) (m) (γ), (12)
where

B δ (µ) = {γ : d λ,M (µ, γ) ≤ δ}.
It remains for us to identify [START_REF] Faugeras | Asymptotic description of stochastic neural networks[END_REF] with the rate function in (4). We claim that for each δ > 0, lim

m→∞ inf γ∈B δ (µ) I (3) (m) (γ) = inf γ∈B δ (µ) (I (3) (γ) -Γ(γ)). (13)
To see this, we have from Lemma 11 that for all m and all γ, and constants α 1 < 1,

α 3 > 1 and α 2 , α 4 ∈ |, (note that if γ /
∈ E 2 the inequalities below are immediate from the definitions)

I (3) (γ) -Γ(γ) ≥ (1 -α 1 )I (3) (γ) -α 2 (14)
I (3) (m) (γ) ≥ (1 -α 1 )I (3) (γ) -α 2 (15)
I (3) (γ) -Γ(γ) ≤ (1 + α 3 )I (3) (γ) -α 4 I (3) (m) (γ) ≤ (1 + α 3 )I (3) (γ) -α 4
We thus see that if I (3) (γ) = ∞ for all γ ∈ B δ (µ), then ( 13) is identically infinite on both sides. Otherwise, it may be seen from ( 14) and ( 15) that it suffices to establish [START_REF] Faugeras | Large deviations of a spatially-ergodic neural network with learning[END_REF] in the case that

B δ (µ) = B δ l (µ) := {γ : d λ,M (µ, γ) ≤ δ, I (3)
(γ) ≤ l} for some l < ∞. But it follows from (29) and (34) that for all γ ∈ B δ l (µ), there exist constants (α m 5 ) which converge to 0 as m → ∞ and such that I [START_REF] Faugeras | Large deviations of a spatially-ergodic neural network with learning[END_REF]. The expression for the rate function in (4) now follows since I (3) (γ) -Γ(γ) is lower semicontinuous, by Lemma 12.

(3) (m) (γ) -I (3) (γ) + Γ(γ) ≤ α m conclude
For the second statement in the Theorem, if

I (3) (µ • β -1 ) = 0, then I (2) (µ • β -1
) Vn ||P ⊗Vn = 0 for all n. However since the relative entropy has a unique zero, this means that (µ • β -1 ) Vn = P ⊗Vn for all n. But this means that µ • β -1 = P ) d , and therefore (using Lemma 6) µ = Q.

Fourier Series, and be such that the eigenvalues of K(θ) are strictly greater than zero for all θ. We require that K is the density of a stationary sequence, which means that we must also assume that for all θ

K(-θ) = † K(θ) = K(θ).
This means, in particular, that K(θ) is Hermitian. We write

(∆(ω)) k = j∈) d † S j ω k-j , where (16)
j∈) d S j exp (-i j, θ ) = K(θ) -1 2 . ( 17
)
Here K-1 2 is understood to be the positive Hermitian square root of K-1 . The Fourier Series of K- 1 2 is absolutely convergent as a consequence of Wiener's Theorem. In this section, we determine a general expression for I (3) (ξ • ∆ -1 ). We are generalising the result for b = 1 with )-indexing given in [START_REF] Donsker | Large deviations for stationary Gaussian processes[END_REF]. These results are necessary for the proofs in the previous section.

We similarly write that

(Υ(ω)) k = j∈) d † R j ω k-j , where (18)
j∈) d R j exp (-i j, θ ) = K(θ) 1 2 . ( 19
)
As previously, K 1 2 is understood to be the positive definite Hermitian square root of K. We note that R -j = † R j and S -j = † S j . Lemma 6. For all ξ ∈ E 2 , ξ • ∆ -1 and ξ • Υ -1 are in E 2 and

ξ • ∆ -1 • Υ -1 = ξ • Υ -1 • ∆ -1 = ξ.
Proof. We make use of the following standard Lemma from [START_REF] Shiryaev | Probability[END_REF], to which the reader is referred for the definition of an orthogonal stochastic measure. Let (U j ) ∈ R b be a zero-mean stationary sequence governed by ξ ∈ E 2 . Then there exists an orthogonal

| b -valued stochastic measure Z ξ = Z ξ (∆) (∆ ∈ B([-π, π[ d ))
, such that for every j ∈ ) d (ξ a.s.)

U j = 1 2π [-π,π] d exp(i j, θ )Z ξ (dθ). (20)
Conversely any orthogonal stochastic measure defines a zero-mean stationary sequence through (20). It may be inferred from this representation that

Z ξ•∆ -1 (dθ) = † K(θ) -1 2 Z ξ (dθ), Z ξ•Υ -1 (dθ) = † K(θ) 1 2 Z ξ (dθ).
The proof that this is well-defined makes use of the fact that K 1 2 and K-1 2 are uniformly continuous, since their Fourier Series' each converge uniformly. This gives us the Lemma. We note for future reference that, if ξ has spectral measure d ξ(θ), then the spectral density of

ξ • ∆ -1 is K-1 2 (θ)d ξ(θ) K-1 2 (θ). (21)
It remains for us to determine a specific expression for

I (3) (ξ • Υ -1 ). Definition 5. If ξ ∈ E 2 , we define Γ ∆ (ξ) = 1 2 E ξ ω 0 2 -E ξ•∆ -1 ω 0 2 - 1 2(2π) d [-π,π] d log det K(θ) dθ. (22) Otherwise we define Γ ∆ (ξ) = 0. Lemma 7. If ξ ∈ E 2 , Γ ∆ (ξ) = 1 2(2π) d [-π,π] d tr Id b -K(θ) -1 d ξ(θ) - 1 2(2π) d [-π,π] d log det K(θ) dθ. (23)
Proof. We see from (21) that ξ • ∆ -1 has spectral measure K-1 2 (θ)d ξ(θ) K-1 2 (θ). We thus find that

Γ ∆ (ξ) = 1 2(2π) d [-π,π] d tr d ξ(θ) - 1 2(2π) d [-π,π] d tr K-1 2 (θ)d ξ(θ) K-1 2 (θ) - 1 2(2π) d [-π,π] d log det K(θ) dθ. = 1 2(2π) d [-π,π] d tr d ξ(θ) - 1 2(2π) d [-π,π] d tr K(θ) -1 d ξ(θ) - 1 2(2π) d [-π,π] d log det K(θ) dθ. Lemma 8. For all ξ ∈ E 2 , I (3) ξ • ∆ -1 ≤ I (3) (ξ) -Γ ∆ (ξ).
Proof. We assume for now that there exists a q such that S j = 0 for all |j| ≥ q, denoting the corresponding map by ∆ q . Let N m q : T Vm → T Vm be the following linear operator. For j ∈ V m , let (N m q ω) j = k∈Vm S (k-j) mod Vm ω k . Let ξq,n = ξ Vn • (N n q ) -1 . It follows from this assumption that the V l marginals of ξq,n and ξ • ∆ -1 q are the same, as long as l ≤ nq. Thus

I (2) (ξ • ∆ -1 q ) V l ||P ⊗V l =I (2) ( ξq,n ) V l ||P ⊗V l ≤I (2) ξq,n ||P ⊗Vn . ( 24
)
This last inequality follows from a property of the Relative Entropy I (2) , namely that it is nondecreasing as we take a 'finer' σ-algebra (it is a direct consequence of [9, Lemma 2.3]). If ξ Vn does not have a density for some n then I (3) (ξ) is infinite and the Lemma is trivial. If otherwise, we may readily evaluate I (2) ξq,n ||P ⊗Vn using a change of variable to find that

I (2) ξq,n ||P ⊗Vn = I (2) ξ Vn ||P ⊗Vn + 1 2 E ξ Vn N n q ω 2 -ω 2 + 1 2 log det N n q .
We divide (24) by (2l + 1) d , substitute the above result, and finally take l → ∞ (while fixing n = l + q) to find that

I (3) ξ • ∆ -1 q ≤ I (3) (ξ) + 1 2 E ξ•(∆) -1 ω 0 2 -E ξ ω 0 2 + 1 2(2π) d [-π,π] d log det K(θ) dθ = I (3) (ξ) -Γ ∆ q (ξ).
Here Γ ∆ q (ξ) is equal to Γ ∆ (ξ) as defined above, subject to the above assumption that S j = 0 for |j| > q. On taking q → ∞, it may be readily seen that Γ ∆ q → Γ ∆ pointwise. Furthermore the lower semicontinuity of I (3) dictates that

I (3) ξ • ∆ -1 ≤ lim q→∞ I (3) ξ • ∆ -1 q ,
which gives us the Lemma.

Lemma 9. For all ξ ∈ E 2 , I (3) (ξ • ∆ -1 ) = I (3) (ξ) -Γ ∆ (ξ).
Proof. We find, similarly to the previous Lemma, that if γ ∈ M s (T ) d ) then

I (3) γ • Υ -1 ≤ I (3) (γ)+ 1 2 E γ•Υ -1 ω 0 2 -E γ ω 0 2 - 1 2(2π) d [-π,π] d log det K(θ)dθ. ( 25
)
We substitute γ = ξ • ∆ -1 into the above and, after noting Lemma 6, we find that

I (3) (ξ) ≤ I (3) ξ • ∆ -1 + 1 2 E ξ ω 0 2 -E (ξ•∆ -1 ) ω 0 2 - 1 2(2π) d [-π,π] d log det K(θ)dθ = I (3) ξ • ∆ -1 + Γ ∆ (ξ). ( 26
)
The result now follows from the previous Lemma and (26).

We next prove some matrix identities which are needed in the proof of Lemma 11.

Lemma 10. If A, B ∈ C l×l (for some positive integer l) are both Hermitian, then

tr(AB) ≤ b A B ≤ l A tr(B) if B is positive ltr(A) B if A is positive (27)
If A and B are positive and in addition A is invertible then

tr (B) ≤ l A -1 2 tr(ABA). ( 28
)
Proof. The first part of ( 27

< a < 1 2 , aE µ [ ω 0 2 ] + b 2 log(1 -2a) ≤ I (3) (µ). (29)
For all µ ∈ E 2 , there exist constants

α 1 < 1, α 2 ∈ |, α 3 > 1, α 4 ∈ | and m > 0 such that for all m > m, Γ (m) (µ) ≤ α 1 I (3) (µ) + α 2 , (30)
Γ(µ) ≤ α 1 I (3) (µ) + α 2 . (31)
Γ (m) (µ) ≥ -α 3 I (3) (µ) + α 4 (32)
Γ(µ) ≥ -α 3 I (3) (µ) + α 4 (33)
There exist constants α m 3 , α m 4 > 0 which converge to zero as m → ∞, and such that

Γ (m) (µ) -Γ(µ) ≤ α m 3 E µ [ ω 0 2 ] + α m 4 . (34)
Proof. It is a standard result that if µ / ∈ M s (T ) d ), then I (3) (µ) = ∞, for which the first result is evident. Let µ ∈ M s (T ) d ). For w ∈ T Vn we let f (ω) = k∈Vn ω k 2 . The function f M (ω) = f (ω)1 af (ω)
≤M is bounded, and hence from definition 1 we have

a T Vn f M dµ Vn ≤ log T Vn exp(af M ) dP ⊗Vn + I (2) (µ Vn ||P ⊗Vn ).
We obtain using an easy Gaussian computation that log

T Vn exp(af ) dP ⊗Vn = - (2n + 1) d b 2 log(1 -2a).
Upon taking M → ∞ and applying the dominated convergence theorem we obtain

aE µ [ ω 0 2 ] ≤ - b 2 log(1 -2a) + 1 (2n + 1) d I (2) (µ Vn ||P ⊗Vn ).
We have used the stationarity of µ. By taking the limit n → ∞ we obtain the first inequality (29).

It follows from the definition that (30)-( 33) are true if µ / ∈ E 2 . Thus we may assume that µ ∈ E 2 . We choose m to be such that the eigenvalues of F(m) (as defined in ( 8)) are strictly greater than zero for all m > m. It may be easily verified that

b s=1 μss [-π, π] d = (2π) d E µ ω 0 2 . (35)
We observe the following upper and lower bounds which hold for all m > m (and for Γ 1 too),

Γ 1,(m) ≤ - 1 2 inf m≥ m,θ∈[-π,π] d log det G(m) (θ) < ∞ (36) Γ 1,(m) ≥ - 1 2 sup m≥ m,θ∈[-π,π] d log det G(m) (θ) > -∞ (37)
We recall that, since Lemma 10,(28), to this, obtaining

G(m) = F 2 (m) , Γ 2,(m) (µ) = 1 2(2π) d [-π,π] d tr (dμ(θ)) - 1 2(2π) d [-π,π] d tr F-2 (m) (θ)dμ(θ) . Note that tr( F -2 (m) (θ)dμ(θ)) = tr( F -1 (m) (θ)( F -1 (m) (θ)dμ(θ) F -1 (m) (θ)) F(m) (θ)) = tr( F -1 (m) (θ)dμ(θ) F -1 (m) (θ)), apply
Γ 2,(m) (µ) ≤ 1 2 (1 -α * 1 ) E µ ω 0 2 , where α * 1 = 1 b inf θ∈[-π,π] d ,m≥ m F(m) (θ) -2 > 0.
If α * 1 ≥ 1, then (30) is clear because I (3) (µ) ≥ 0, the above inequality would imply that Γ 2,(m) is negative and (36). Otherwise, we may use ( 29) and (36) to find that

Γ (m) (µ) ≤ 1 2a (1 -α * 1 )I (3) (µ) - b(1 -α * 1 ) 4a log(1 -2a) - 1 2 inf m≥ m,θ∈[-π,π] d log det G(m) .
We may substitute a >

1 2 (1 -α * 1 ), letting α 1 = 1 2a (1 -α * 1
), into the above to obtain (30). The second inequality (31) follows by taking m → ∞ in the first.

For the third inequality, we find using (27) that

Γ 2,(m) (µ) ≥ 1 2 (1 -α * 3 ) E µ [ ω 0 2 ],
where

α * 3 = b inf θ∈[-π,π] d ,m≥ m F(m) (θ) 2 .
If α * 3 ≤ 1, then (32) is clear because of the fact that I (3) (µ) ≥ 0, the above inequality would mean that Γ 2,(m) is positive, and (37). Otherwise, we may use ( 29) and (37) to find that

Γ 2,(m) (µ) ≥ 1 -α * 3 2a I (3) (µ) + b(1 -α * 3 ) 4a log(1 -2a) - 1 2 sup m≥ m,θ∈[-π,π] d log det G(m) (θ).
This yields (32) on taking a < α * 3 -1

2 . Taking limits as m → ∞ yields (33). Now, making use of Lemma 10,(27), it may be seen that

Γ 2,(m) (µ) -Γ 2 (µ) = 1 2(2π) d [-π,π] d tr G-1 (θ) -G-1 (m) (θ) dμ(θ) . ≤ b 2(2π) d G m [-π,π] d tr(dμ(θ)).
Here

G m = sup θ∈[-π,π] d G-1 (θ) -G-1 (m) (θ)
. The convergence of Γ 1,(m) to Γ 1 is clear. We thus obtain (34).

Lemma 12. I (3) (µ) -Γ(µ) is lower-semicontinuous.

Proof. Since I (3) (µ) -Γ(µ) is infinite for all µ / ∈ E 2 , we only need to prove the Lemma for µ ∈ E 2 . We need to prove that if µ (j) → µ then lim j→∞

I (3) (µ (j) • β -1 ) ≥ I (3) (µ • β -1 ).
We may assume without loss of generality that

lim j→∞ I (3) (µ (j) • β -1 ) = lim j→∞ I (3) (µ (j) • β -1 ) ∈ R ∪ ∞. If E µ (j) [ v 0 2
] → ∞, then by (29) in Lemma 11, lim j→∞ I (3) (µ (j) • β -1 ) = ∞, satisfying the requirements of the Lemma. Thus we may assume that there exists a constant l such that E µ (j) [ ω 0 2 ] ≤ l for all j. We therefore have that, for all m, because of ( 11) and ( 4),

lim j→∞ I (3) (µ (j) • β -1 ) = lim j→∞ I (3) (m) µ (j) + Γ (m) µ (j) -Γ µ (j) .
Making use of (32), we may thus conclude that

lim j→∞ I (3) (µ (j) • β -1 ) ≥ lim j→∞ I (3) (m) (µ (j) ) -ǫ * (m) ,
for some ǫ * (M ) which goes to 0 as m → ∞. In addition, lim j→∞ I

(m) (µ (j) ) ≥ I

(m) (µ) due to the lower semi-continuity of I

(m) . On taking m → ∞, since I

(m) (µ) → I (3) (µ • β -1 (3)
), we find that lim j→∞ I (3) (µ (j) • β -1 ) ≥ I (3) (µ (j) • β -1 ).

B. A lemma on the Large Deviations of Stationary Random Variables

The following lemma is an adaptation of [2, Theorem 4.9] to ) d . We state it in a general context. Let B be a Banach Space with norm • . For j ∈ ) d , let

λ j = 3 -d d δ=1 2 -|j(δ)| . We note that k∈) d λ k = 1. Define the metric d λ on B ) d by d λ (x, y) = j∈) d λ j min x j -y j , 1 .
Let the induced Prohorov metric on M(B ) d ) be d λ,M . For ω ∈ B ) d and j ∈ ) d , we define the shift operator as S j (ω) k = ω j+k . Let B δ (A) = {x ∈ B ) d : d λ (x, y) ≤ δ for some y ∈ A} be the closed blowup, and B(δ) be the closed blowup of {0}.

Suppose that for m ∈ ) + , Y (m) , W ∈ B ) d are stationary random variables, governed by a probability law P. We suppose that μn (Y (m) ) is governed by Π n (m) and μn (W ) is governed by Π n W -these being the empirical process measures, defined analogously to (2). Suppose that for each m, (Π n (m) ) satisfies an LDP with good rate function J (m) . Suppose that W = Y (m) + Z (m) for some series of stationary random variables Z (m) on B ) d . Lemma 13. If there exists a constant κ > 0 such that for all b > 0 lim

m→∞ lim n→∞ 1 (2n + 1) d log E exp b j∈Vn min Z j (m) 2 , 1 < κ, (38)
then (Π n W ) satisfies an LDP with good rate function

J(x) = lim ǫ→0 lim m→∞ inf y∈B ǫ (x) J (m) (y) = lim ǫ→0 lim m→∞ inf y∈B ǫ (x) J (
1 (2n + 1) d log P d λ,M (μ n (W ), μn (Y (m) )) > ǫ = -∞. (39) For x ∈ B ) d , write |x| λ := d λ (x, 0). Let B ∈ B(B ) d ).
Then, noting the definition of p n just above (2),

μn (W )(B) = 1 (2n + 1) d j∈Vn 1 B S j p n (Y (m) ) + S j p n (Z (m) ) ≤ 1 (2n + 1) d j∈Vn 1 B S j p n (Y (m) ) + S j p n (Z (m) ) 1 B(ǫ) (S j p n (Z (m) )) + 1 B(ǫ) c S j p n (Z (m) ) ≤ 1 (2n + 1) d j∈Vn 1 B ǫ (S j p n (Y (m) )) + 1 (2n + 1) d #{j ∈ V n : |S j p n (Z (m) )| λ > ǫ} ≤ μn (Y (m) )(B ǫ ) + 1 (2n + 1) d # j ∈ V n : |S j p n (Z (m) )| λ > ǫ . Thus P d λ,M (μ n (W ), μn (Y (m) )) > ǫ ≤ P 1 (2n + 1) d # j ∈ V n : |S j p n (Z (m) )| λ > ǫ > ǫ ≤ P j∈Vn |S j p n (Z (m) )| 2 λ > (2n + 1) d ǫ 3 ≤ exp -b(2n + 1) d ǫ 3 E P exp b j∈Vn |S j p n (Z (m) )| 2 λ ≤ exp -b(2n + 1) d ǫ 3 E P   exp   b j∈Vn,k∈) d λ k min( p n (Z (m) ) j+k 2 ,
I (3) ((Ψ ) d ) -1 (µ)) -Γ((Ψ ) d ) -1 (µ)). (41)
Clearly (Ψ ) d ) -1 (µ) is in E 2 if and only if µ is in E 2 . Let ν = (Ψ ) d ) -1 (µ). It is well known that if (Ψ ) d ) -1 (µ) Vn is absolutely continuous relative to P ⊗Vn , then the relative entropy may be written as Otherwise the relative entropy is infinite. Thus if the relative entropy is finite, (Ψ ) d ) -1 (µ) Vn must possess a density, and this means that µ Vn possesses a density which we denote by r(x) : T Vn → T Vn . We note that the density of (Ψ(P )) ⊗Vn is ρ Vn (x) = (2π) - We divide by (2n + 1) d and take n to infinity to obtain

I (3) ((Ψ ) d ) -1 (µ)) = I (3) (µ) - † m µ c + 1 2 c 2 .
If µ Vn does not possess a density for some n, then both sides of the above equation are infinite. It may be verified that the spectral density of ν is given by

dν(θ) = dμ(θ) + (2π) d δ(θ)(c † c -m µ † c -c † m µ ).
On substituting this into the expression in Theorem 1, we find that

Γ 2 (ν) =Γ 2 (µ) + 1 2 tr Id b -G(0) -1 c † c - 1 2 tr Id b -G(0) -1 m µ † c + c † m µ =Γ 2 (µ) + 1 2 † c Id b -G(0) -1 c - † c Id b -G(0) -1 m µ .
We have used the fact that G(0) is symmetric. We thus obtain [START_REF] Chiyonobu | The large deviation principle for hypermixing processes[END_REF]. This minimum of the rate function remains unique because of the bijectivity of (Ψ ) d ).

Definition 3 .

 3 A sequence of probability laws (Γ n ) on some topological space Ω equipped with its Borelian σ-algebra is said to satisfy a strong large deviation principle (LDP) with rate function I : Ω → | if I is lower semicontinuous, for all open sets O, lim n→∞ 1 (2n + 1) d log Γ n (O) ≥inf x∈O I(x) and for all closed sets F lim n→∞ 1 (2n + 1) d log Γ n (F ) ≤inf x∈F I(x).

≤lim m→∞ lim n→∞ 1 (

 1 b > 0. Since k∈) d λ k = 1 and the exponential function is convex, by Jensen's Inequalityexp -b(2n + 1) d ǫ 3 E P   exp   b j∈Vn,k∈) d λ k min( p n (Z (m) ) j+k 2 , 1) exp -b(2n + 1) d ǫ 3 E P k∈) d λ k exp b j∈Vn min( p n (Z (m) ) j+k 2 , 1) = exp -b(2n + 1) d ǫ 3 E P exp b j∈Vn min( Z j (m) 2 , 1), by the stationarity of Z (m) and the fact that k∈) d λ k = 1. We may thus infer, using (38), that 2n + 1)d log P d λ,M (μ n (W ), μn (Y (m) )) > ǫ ≤ -bǫ 3 + κ. (40)Since b is arbitrary, we may take b → ∞ to obtain (39).C. Proof of Corollary 2We now prove Corollary 2.Proof. Let φ : T → T be φ(ω) = ω + c, and φ) d : T ) d → T ) d be φ ) d (ω) j = φ(ω j ). Let Ψ ) d : M s (T ) d ) → M s (T ) d ) be Ψ ) d (µ) = µ • (φ ) d ) -1 , and Ψ : M s (T ) → M s (T ) be Ψ(ν) = ν • φ -1 .It is easily checked that these maps are bicontinuous bijections for their respective topologies. SinceΠ n c = Π n • (Ψ ) d ) -1 ,we have by a contraction principle [6, theorem 4.2.1] that (Π n c ) satisfies a strong LDP with good rate function

I ( 2 )

 2 (((Ψ ) d ) -1 (µ)) Vn ||P ⊗Vn ) = E ((Ψ ) d ) -1 (µ)) Vn (x) log d((Ψ ) d ) -1 (µ)) Vn dP ⊗Vn (x) = E µ Vn (x) log dµ Vn d(Ψ(P )) ⊗Vn (x) = I (2) (µ Vn ||(Ψ(P )) ⊗Vn ).

I ( 2 )

 2 µ Vn ||(Ψ(P )) ⊗Vn = E µ Vn (x) log r(x) ρ Vn (x) = I (2) µ Vn ||P ⊗Vn -(2n + 1) d † m µ c + (2n + 1) d 2 c 2 .

  ) follows from the von Neumann's trace inequality. Given two matrices A and B in C l×l where the α t 's and β t 's are the singular values of A and B. In the case where A and B are Hermitian the singular values are the magnitudes of the (real) eigenvalues. By ≤ l A B . If B is positive, so are its eigenvalues and B ≤ tr(B), hence (27). If A is invertible, tr(B) = tr(A -2 ABA). If moreover and A and B are both Hermitian positive, we obtain the second identity by applying (27) to the Hermitian matrix A -2 and the Hermitian positive matrix ABA.

	l
	|tr(AB)| ≤	α t β t ,
	t=1
	Cauchy-Schwartz, l t=1 α t β t ≤ t Lemma 11. For all 0 l t=1 α 2 t × l t=1 β 2

Acknowledgements: This work was partially supported by the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 269921 (BrainScaleS), no. 318723 (Mathemacs), and by the ERC advanced grant NerVi no. 227747. A. Properties of the Entropy Let K : [-π, π] d → C b×b possess an absolutely convergent