Kevin Belkacem

Transport of angular momentum by waves in stars

Keywords: Waves -Stars, oscillations -Stars, interiors -Stars, rotation -Stars, evolution

may come

Transport of angular momentum by waves in stars

Introduction

Rotation has fundamental effects on stellar evolution. For instance, it induces meridional circulations, as well as shear and baroclinic instabilities, which contribute to the redistribution of angular momentum and to the mixing of chemical elements (e.g. [START_REF] Maeder | Physics, Formation and Evolution of Rotating Stars[END_REF]). Therefore, it is essential to properly understand and model the physical processes responsible for the transport and redistribution of angular momentum. Among them, waves1 have been identified to play a major role. Indeed, the wave/mean-flow interactions has been considered for a long time for geophysical flows and more precisely the interaction between waves, angular momentum, and meridional circulation. This had been extensively studied in the 60s and 70s in the context of middle atmosphere dynamics (see for instance [START_REF] Andrews | Middle atmosphere dynamics[END_REF][START_REF] Holton | An introduction to dynamic meteorology[END_REF] for extensive reviews). All those efforts have permitted to unveil and clarify the nature of these interactions and led to the developments of adapted formalisms among which the Transformed Eulerian Mean formalism (Andrews & McIntyre 1976, 1978a) and the Generalized Lagrangian Mean formulation by Andrews & McIntyre (1978b). For stellar interiors, it has been first addressed from a theoretical point of view with the pioneer works of Press (1981) and Ando (1983). Both authors were mainly interested in addressing the problem of the redistribution of angular momentum by waves (even if with quite different motivations) and thus investigated the question of the interaction between waves and rotation through wave momentum stresses in the mean angular momentum equation. Subsequently, two set of distinct works had been developed almost in parallel.

On one had, following Ando (1983), efforts have been made to address the problem of the redistribution of angular momentum in massive stars by unstable normal modes (e.g. Ando1986, Lee & Saio 1993, Lee 2007[START_REF] Townsend | IAU Symposium[END_REF][START_REF] Ishimatsu | ASP Conference Proceedings[END_REF][START_REF] Townsend | IAU Symposium[END_REF], Lee et al. 2014[START_REF] Townsend | [END_REF]. As an illustration, one of the main motivations was to explain the episodic mass loss in Be stars. The more studied scenario can be summarized as follows; angular momentum is deposited at the surface by normal modes. Thus, the rotation locally increases to reach the break-up limit and then allows for mass loss. However, this scenario is not commonly accepted and, for instance, was questioned by [START_REF] Ishimatsu | ASP Conference Proceedings[END_REF] and [START_REF] Shibahashi | IAU Symposium[END_REF] because observations show that the break-up velocity is not reached. Consequently, the authors proposed an alternative scenario in which g-modes transfer angular momentum in the stellar surface, rotation increases (but still below the break-up limit), the critical frequency also increases, thus implying mode leakage and transfer of angular momentum to the disk. Obviously, those scenarios crucially depend on the way angular momentum is deposited at the surface by waves. Therefore, it is mandatory to properly decipher the interaction between waves and rotation.

On the other hand, and following Press (1981), many authors have considered the interaction between the internal gravity waves and the mean flow in low-mass stars (e.g. Schatzman 1993;[START_REF] Zahn | [END_REF][START_REF] Kumar | [END_REF]Charbonnel & Talon 2005;Mathis et al. 2013). The primary motivation was to explain the quasi-uniform rotation profile in the Solar radiative interior. Internal gravity waves have been shown to be able to extract angular momentum on short time scales (compared to the evolution time-scale) and had then been considered as a serious candidate to explain the solar rotation profile as well as the cool side of the Li-dip (e.g. Talon & Charbonnel 2003,2005;Charbonnel & Talon 2005). The advent of asteroseismology further strengthened the focus on internal gravity waves as it was a potential candidate for explaining the weakly increasing rotation contrast of subgiant stars and the spin-down of red giants. It has been shown that they are inefficient for evolved red giants to explain the slow-down of red giants (e.g. Fuller et al. 2014;Pinco ¸n et al. 2017) but could explain the rotation profile of subgiants (Pinco ¸n et al. 2017).

Those two set of works, while apparently disconnected, rely essentially on the same physics which require a deep understanding of the interaction between the waves and the mean flow. This is, however, not a trivial issue as we will discuss in the following. In particular, we will show that considering the divergence of the wave flux in the momentum equation is not enough and that the energy equation (more precisely the effect of the wave heat flux) is essential to account for as both equations are intricately coupled by meridional circulation. This coupling had however been too often overlooked in stellar physics, while it has been demonstrated to play a fundamental role. Our objective is thus to highlight this issue. To that end, based on the classical azimutaly-averaged equations we will first present in Sect. 2 some arguments to show the need for considering the coupling between the momentum and energy equation. In Sect. 3, the Transformed Eulerian Mean formalism will be introduced and an application to mixed modes in red giants will be presented in Sect. 4. Finally,in Sect. 3.3, concluding remarks are provided.

Setting the stage

Wave effects on the mean field: heuristic arguments

In stars, progressive or standing waves are generally considered to have small amplitudes, as (but not only) for solar-like oscillations. This is well justified and naturally motivates to adopt the linear approximation. However, given their small amplitudes, it is worthwhile to wonder if such waves can have significant effects on the mean flow and are able to modify the rotation profile. To have some hints about this issue, it is very enlightening to use some simple heuristic arguments. To do so, we here essentially follow the discussion as presented by Buhler (2009). Let us first consider the conservation of momentum equation in a schematic form as

∂ v ∂t + L ( v) + B ( v, v) = 0 , ( 1
)
where v is the velocity field, L a linear operator, and B a bi-linear operator. The velocity field can then be expanded using the wave amplitude (a w 1) as a small parameter, thus

v = v 0 + a w v 1 + a 2 w v 2 + O(a 3 w ) , (2)
where the first r.h.s term ( v 0 ) corresponds to the mean flow, the second r.h.s. term ( v 1 ) stands for the wave velocity, and the third r.h.s term ( v 2 ) is the back reaction onto the mean-flow by the waves. Using Eq. ( 2) into Eq. ( 1) then provides a hierarchical system of equations. The lowest-order is

∂ v 0 ∂t + L ( v 0 ) + B ( v 0 , v 0 ) = 0 , (3)
which corresponds to the equation governing the mean-flow. At the next order one has

∂ v 1 ∂t + L ( v 1 ) + B ( v 0 , v 1 ) + B ( v 1 , v 0 ) = 0 , (4)
which governs the linear waves. Finally, at the second-order, we obtain

∂ v 2 ∂t + L ( v 2 ) + B ( v 0 , v 2 ) + B ( v 2 , v 0 ) = -B ( v 1 , v 1 ) , (5)
which governs the effect of linear waves onto the mean flow. Then, we immediately note that the back-reaction equation is forced by wave-related terms of the second order B ( v 1 , v 1 ) and that no first-order term appears. This is an important point because, as we will see in Sect. 3.3, second-order effects such as Stokes corrections cannot be neglected when considering the effect of waves on the mean flow.

To go further, let us assume that there is no mean field velocity ( v 0 = 0). Then, Eq. ( 5) reduces to

∂ v 2 ∂t + L ( v 2 ) = R , (6)
with R ≡ -B ( v 1 , v 1 ). Now, assuming that v 2 (0) = 0 and that R is independent of time (for steady waves and time-independent mean properties) one has

v 2 R t = v 2 a 2 w t . (7)
The latter solution can then be interpreted as follows; on a large time-scale, t = O (a -2 w ), we have v 2 = O (1). Otherwise stated, waves can significantly modify the mean flow on long time scales.

Azimuthal (or zonal) Eulerian averaging of the primitive equations

In this section, our objective is to derive the equations that describe the effect of waves on the mean flow or more precisely the interaction between the mean flow and the waves. Thus, we start by considering the continuity, momentum, and energy equations in an inertial frame. They can be written such as

∂ρ ∂t + ∇ • (ρ v) = 0 (8) ∂ (ρ v) ∂t + ∇ • (ρ v v) = -∇p -ρ ∇Φ + X (9) ∂ (ρs) ∂t + ∇ • (ρ vs) = Q , ( 10
)
where ρ is the density, v the velocity field, p the pressure, s the specific entropy, Φ the gravitational potential, X is a non-conservative mechanical forcing (e.g., turbulent dissipation), and Q represents heating or cooling terms. Obviously, Eqs. ( 8) to (10) need to be complemented by the appropriate equation of state together with boundary conditions.

To go further, we need to derive the equation representing the conservation of the specific angular momentum. This equation is obtained by adopting spherical coordinates (r, θ, φ) with the associated unit vectors ( e r , e θ , e φ ) and θ = 0 corresponding to the direction of the rotation axis. Then, using the azimuthal component of Eq. ( 9), multiplied by ≡ r sin θ, leads to

∂ (ρh) ∂t + ∇ • (ρh v) = - ∂p ∂φ -ρ ∂Φ ∂φ + X φ , (11)
where the specific angular momentum is defined by h ≡ v φ and v φ is the azimuthal component of the velocity field.

To gain insight and tractability to the problem, it is useful to shift from a 3D to a 2D problem and to separate the waves and the mean-flow. To do so, we consider azimuthal (or zonal) average of the equations. Therefore, a given field can be decomposed into a mean part and a perturbation. More precisely, for a given field A, we have

A = A + A , ( 12
)
where A is the Eulerian-mean azimuthal average defined by

A = 1 2π 2π 0 A dφ . (13)
This Eulerian average then possesses the usual properties (additivity, linearity, commutativity with partial differentiation, etc. . . ) and A stands for the non-axisymmetric pertubations. Applied to the velocity field, this decomposition permits us to identify the main ingredients of the problem. More precisely, one has

v = v + v , (14)
where v = v ⊥ + Ω e φ with v ⊥ corresponding to the meridional circulation. Finally, the nonaxisymmetric perturbations are associated with non-axisymmetric waves.

Because we are interested in estimating the effect of waves on the mean flow (and primarily on h), it is necessary to obtain the mean equations. This is done by applying Eq. ( 12) to Eqs. ( 8), (10), and (11) and finally to perform an azimuthal averaging. After some manipulations, and retaining perturbations up to the second order, it gives ∂ρ ∂t

+ ∇ ⊥ • ρ v ⊥ = D , (15)
ρ ∂h ∂t + ρ v ⊥ • ∇ ⊥ h = -∇ ⊥ • ρ v φ v ⊥ + X φ + H , (16)
ρ ∂s ∂t + ρ v ⊥ • ∇ ⊥ s = -∇ ⊥ • ρ s v ⊥ + Q + S , (17)
with

H = -ρ v ⊥ • ∇ ⊥ h - ∂ρ v φ ∂t -ρ ∂Φ ∂φ -∇ ⊥ • ρ v φ v ⊥ , (18)
S = -ρ v ⊥ • ∇ ⊥ s -∇ ⊥ • ρ s v ⊥ - ∂ρ s ∂t , (19)
D = -∇ ⊥ • ρ v ⊥ , (20)
and ∇ ⊥ and v ⊥ are the gradient and velocity vector in the meridional plane, respectively. They are defined by

∇ ⊥ = e r ∂ ∂r + e θ 1 r ∂ ∂θ , and v ⊥ = v r e r + v θ e θ . (21)
Looking at Eqs. ( 15) to ( 17), one immediately note that wave-related terms are only of the second order. This is consistent with the heuristic arguments presented in Sect. 2.1. Consequently, it also implies that the wave field can modify the mean-flow only on large time-scales (compared to the wave periods). Moreover, anticipating on the following, one can already note that wave-related secondorder effects such as the Stoke drifts cannot be neglected.

Equations ( 15) to ( 17) also show that many wave-related terms are potentially able to affect the mean flow. Among them, one can distinguish terms involving a time derivative. They are nonnegligible only if the wave amplitude varies with time or if one consider the effect of wave packets. In the latter case those terms are nevertheless expected to be small because of cancellation effects. In contrast, if we consider steady waves (i.e. wave with a constant amplitude), those terms vanish. Such a steady state is often assumed when considering progressive as well as stationary waves in stars. We will thus neglect those terms.

Finally, all terms encompassed into H, S, D exhibit density perturbation (ρ ). They are generally considered to be small (e.g. Ando 1983[START_REF] Unno | Nonradial oscillations of stars[END_REF]. This is particularly the case for lowfrequency waves (σ R N , where σ R is the wave frequency and N the buoyancy frequency) where the anelastic approximation applies (e.g. Dintrans & Rieutord 2001). Therefore, we neglect the terms ρ v ⊥ , ρ s , and ρ v φ in Eqs. ( 18) to (20). In addition, we use the Cowling approximation by neglecting the perturbation of the gravitational potential, so that Eqs. ( 15) to ( 17

) become ∂ρ ∂t + ∇ ⊥ • ρ v ⊥ = 0 (22) ρ ∂h ∂t + ρ v ⊥ • ∇ ⊥ h + ∇ ⊥ • ρ v φ v ⊥ = X φ (23) ρ ∂s ∂t + ρ v ⊥ • ∇ ⊥ s + ∇ ⊥ • ρ s v ⊥ = Q . ( 24
)
These equations immediately show that wave forcing of the mean flow is not only due to the wave momentum flux in Eq. ( 23) (also sometimes named to as the Reynolds stresses in analogy to the turbulent Reynolds stress) but also to the wave heat flux in the mean entropy equation, Eq. ( 24). It is important to already note that the latter cannot be neglected compared to the former, especially for low-frequency waves. This will be discussed in details in the following sections. Moreover, we note that our system of equation is not closed because the mean specific entropy s and mean specific angular momentum h are not independent. They are connected through the baroclinic equation (also known to as the thermal wind balance equation). It is obtained by taking the curl of the hydrostatic equilibrium equation so that

ρ 2 ∇ ⊥ h 2 4 × ∇ ⊥ + ∇ ⊥ ρ × ∇ ⊥ p = 0 , (25)
together with the equation of state s = s (ρ, p). This equation is fundamental because it determines the balance between the centrifugal (first term of the l.h.s. of Eq. 25) and baroclinic (second term of the l.h.s. of Eq. 25) torques. Indeed, to ensure the torques balance, the thermal wind equation states that Ω(r) is completely determined by ρ(r) along an isobar. In other words, Ω(r) and ρ(r) (therefore s(r)) cannot evolve separately because they are connected by the meridional circulation. Finally, in principle solving Eq. ( 22) to (25) (with the appropriate equation of state and boundary conditions), together with the wave equations, permits one to properly model the wave and mean flow interactions. This is however a non-trivial task. First, the time-scales of the problem are not always commensurable in the sense that we need to solve the wave equations on a short time-scale (typically a dynamical time-scale) while we need to solve mean-equations on a longer time-scale (typically the Kelvin-Helmholtz time-scale). Second, and maybe more importantly, this system of equations made very difficult to make clear what physical mechanism is responsible for the mean-flow driving. This is illustrated by the fact that the second and third terms of the l.h.s. of Eq. ( 23) and Eq. ( 24) often nearly compensate each others so that the resulting driving of the mean flow is a residual.

Coupling between wave fluxes and meridional circulation: an essential ingredient

Before going further, it is necessary to insist on a crucial point: the wave heat flux must not be neglected, as well as meridional circulation. If neglected, most of the physics of the problem is lost and this potentially leads to incorrect estimates of the effect of waves on the mean-flow and subsequently on the evolution of the rotation profile.

To figure out the problem, let us consider the conservation of angular momentum and neglect the wave heat flux, meridional circulation, and any external mechanical forcing. After meridional averaging (denoted by . ), it gives

ρ ∂h ∂t = - 1 r 2 ∂ ∂r r 2 ρ v φ v r . ( 26
)
One then immediately concludes that the wave momentum flux modifies the distribution of angular momentum and thus the rotation profile. Let us now consider adiabatic waves. If the effect of rotation is neglected then the r.h.s of Eq. ( 26) vanishes because v φ and v r are in phase quadrature, thus waves do not modify the rotation profile. However, if we include the effect of rotation onto the waves (even in the adiabatic limit), the r.h.s of Eq. ( 26) does not vanish and one must conclude that adiabatic waves affect the mean flow. This is indeed a weird conclusion because it would mean that adiabatic waves should be able to modify the angular momentum of the mean flow without energy exchanges.

To solve this apparent physical issue, one must consider both the wave heat flux and the meridional circulation

ρ ∂h ∂t + ρ v ⊥ • ∇ ⊥ h = - 1 r 2 ∂ ∂r r 2 ρ v φ v r (27) ρ ∂s ∂t + ρ v ⊥ • ∇ ⊥ s = - 1 r 2 ∂ ∂r r 2 ρ s v r . ( 28
)
Still in the limit of adiabatic waves but considering the effect of rotation on the waves, both terms in the r.h.s of Eq. ( 27) and Eq. ( 28) do not vanish. However, waves do not modify the mean flow as shown by the non-acceleration theorem (Andrews and Mclntyre 1976). This apparent paradox is solved by realizing that the wave fluxes of heat and momentum produce a meridional circulation that cancels their tendency to affect the mean flow. It is thus important to stress on several points;

• the heat and momentum wave fluxes do not act independently of each other.

• the connexion is performed by the meridional circulation, which must ensures that h and s satisfy the baroclinic equation (Eq. 25).

• overlooking the effect of the wave heat flux or the coupling between wave fluxes and meridional circulation does not permit to grasp the main picture of the problem and potentially leads to unphysical or incorrect results.

This non-trivial interplay between the mean-flow and the wave forcing terms was properly addressed decades ago in the context of geophysical flows (e.g. Andrews & McIntyre 1976,1978a,1978b). In contrast, in the stellar context, the picture is made confused because the wave heat flux is too-often overlooked (e.g., Press 1981, Ando 1983, Lee & Saio 1993[START_REF] Zahn | [END_REF][START_REF] Kumar | [END_REF], Pantillon et al 2007, Mathis 2013[START_REF] Townsend | IAU Symposium[END_REF][START_REF] Townsend | [END_REF]) even if not always (e.g. Lee 2013, Belkacem et al. 2015a,b). In the following sections, we therefore introduce two formalisms, which nicely clarify the physical picture, namely: the Transformed Eulerian Mean (TEM) and the Generalized Lagragian Mean (GLM) formalisms.

Transformed Eulerian Mean (TEM) formalism

As discussed in Sect. 2, the wave fluxes of heat and momentum produce a meridional circulation which, in the limit of steady and adiabatic waves, cancels their tendency to affect the mean flow. Therefore, the principle of the TEM is to eliminate the advective part of the divergence of the wave heat flux in the entropy equation so as to inject it into the a newly defined meridional circulation that will be named to as the residual circulation. This formalism was introduced by Andrews & McIntyre (1976,1978a).

Incorporating the divergence of the skew flux into meridional circulation

The first step consists in splitting the divergence of the wave heat flux into an advective and a diffusive part. To do so, let us consider an isentropic surface with a normal vector define as n = ∇ ⊥ s/| ∇ ⊥ s|. Then, the wave heat flux R = s v ⊥ can be split into a component along an isentropic surface (the skew flux) and a component perpendicular to it, such as

R = n × R × n + n • R n . ( 29
)
The divergence of the skew flux can then be rewritten as

∇ ⊥ • n × R × n = ∇ ⊥ • ∇ ⊥ s × R | ∇ ⊥ s| 2 × ∇ ⊥ s = ∇ ⊥ × ∇ ⊥ s × R | ∇ ⊥ s| 2 • ∇ ⊥ s = ˜ v • ∇ ⊥ s , (30)
where we used the relation ∇ • ( a × ∇α) = ∇α • ( ∇ × a), valid for any scalar α and vector a. Equation (30) shows that the skew flux behaves like an advection by the velocity ˜ v. The main motivation underlying the TEM is thus to incorporate the advective part of the wave heat flux into the mean meridional velocity field ( v). To do so, the residual meridional circulation is thus defined by

ρ v † ⊥ = ρ v ⊥ + ∇ ⊥ × ρ ψ e φ , (31)
where the stream function ψ is deduced from Eq. ( 30) as

ψ = ∇ ⊥ s × R | ∇ ⊥ s| 2 • e φ = 1 | ∇ ⊥ s| 2 ∂s ∂r s v θ - 1 r ∂s ∂θ s v r . (32)
Now, inserting Eqs. ( 31) and (32) into Eq. ( 22) to Eq. ( 24), we have

∂ρ ∂t + ∇ ⊥ • ρ v † ⊥ = 0 , (33)
ρ ∂h ∂t + ρ v † ⊥ • ∇ ⊥ h + ∇ ⊥ • ρ F = X φ , (34)
ρ ∂s ∂t + ρ v † ⊥ • ∇ ⊥ s + ∇ ⊥ • ρ G = Q , ( 35
)
where the components of the vectors F and G are given by

F r = v φ v r + ψ r ∂h ∂θ , F θ = v φ v θ -ψ ∂h ∂r , (36)
G r = s v r + ψ r ∂s ∂θ , G θ = s v θ -ψ ∂s ∂r .
(37)

Equation ( 25) is left unmodified. Note that F is often named to as the Eliassen-Palm flux because of the pioneer paper by Eliassen & Palm (1961). Equations ( 33) to ( 35) seem very similar to Eqs. ( 22) -( 24) but clarifies the exact role of the wave flux. To make clear the advantages of the TEM, it is useful to place ourselves in the limit of quasi-shellular approximation.

The TEM within the quasi-shellular approximation

Shellular rotation, an often used approximation in stellar radiative zones, assumes that an efficient horizontal transport of angular momentum is at work so that angular velocity is almost constant on isobars. Without discussing here the validity of this assumption (but see Maeder 2009 for a comprehensive discussion), it is a useful framework that permits us to exhibit the advantages of the TEM equations.

Here we first remind the basic properties of this approximation. To do so, let us first introduce the following decomposition

Ω(r, θ) = Ω 0 (r) + Ω (r, θ) , (38)
with

Ω 0 = π 0 sin 3 θ Ω(r, θ) dθ π 0 sin 3 θ dθ = 3 4 π 0 sin 3 θ Ω(r, θ) dθ, (39)
where, within the shellular approximation, Ω 0 Ω. The scalar quantities are developed as

X(r, θ) = X (r) + X (r, θ) , (40)
with

X = π 0 sin θ X(r, θ) dθ π 0 sin θ dθ = 1 2 π 0 sin θ X(r, θ) dθ . (41)
For slow to moderate rotation, we assume that X X . Following these assumptions, and to have more insight about the advantage of using the TEM formulation, it is useful to emphasize that in stars, we are generally in the situation where the isentropic and isopycnal surfaces are nearly parallel. Therefore, we are in the situation for which we have 1 r ∂s ∂θ ∂s ∂r ,

(42) so that the stream function reduces to

ψ ∂s ∂r -1 s v θ , (43)
which is the same as what is found within the quasi-geostrophic approximation. Then, we immediately show that G θ = 0 and that

G r s v r . ( 44
)
To clarify the advantage of the TEM, let us consider the perturbed equation of the entropy by using Eq. ( 10) and Eq. ( 24). It gives

∂s ∂t + v • ∇ s + v • ∇ s = Q + O a 2 w , (45)
where Q ≡ Q/ρ. Multiplying by s and performing azimuthal average then gives

1 2 ∂s 2 ∂t + 1 2 v ⊥ • ∇ ⊥ s 2 + s v • ∇ ⊥ s = Q + O a 3 w . ( 46
)
The first term corresponds to the growth rate of the wave, which can be assumed to be small except in the driving regions and the second term is of the order of O (a 4 w ) because meridional circulation is O (a 2 w ) as it is predominantly wave-driven. Finally, assuming Q = 0 one gets

s v • ∇ ⊥ s = 0 + O a 3 w , (47)
which means that the entropy wave flux is nearly perpendicular to the mean entropy gradient. In other words, the flux is almost horizontal in the quasi-shellular approximation because using Eq. (42) one gets

s v θ s v r 1 . ( 48
)
Coming back to the TEM formalism and Eqs. ( 33) to ( 35), means that most of the wave heat flux has been removed from the entropy equation and the remaining flux (Eq. 44) is a residual. Consequently, we end up with

ρ ∂h ∂t + ρ v † ⊥ • ∇ ⊥ h + ∇ ⊥ • ρ F = X φ (49) ρ ∂s ∂t + ρ v † ⊥ • ∇ ⊥ s = Q , (50)
Therefore, the TEM formalism allows us to gather in a single equation and in the single term F both the wave momentum and wave heat fluxes. It explicitly shows that wave momentum fluxes and wave heat fluxes do not influence the mean flow separately, but only in the combination given by F .

Figure 1: Schematic sketch of a fluid particle trajectory with an initial position X 0 (t = t 0 ). v L and v ξ are the Lagrangian mean as defined by Eq. ( 52) and the Lagagian velocity as defined by Eq. ( 51).

Finally, when F is specified, the residual velocity v † ⊥ becomes part of the solution of Eqs. ( 49) and ( 50) together with Eq. ( 25). They are strictly equivalent to Eqs. ( 22) -( 24), but the TEM equations enable distinguishing the advective and diffusive parts of the wave heat flux and incorporating the advective component in the mean velocity field. More importantly, Andrews & McIntyre (1978a) showed that ∇ • ρ F only depends on wave dissipation and non-steady terms. It is the non-acceleration theorem. Therefore, the TEM makes the adiabatic and non-adiabatic contributions of waves more explicit, the latter being only able to modify the mean flow.

The TEM as seen by the Generalized Lagragian Mean (GLM) formalism

In this section, we very briefly discuss the Generalized Lagrangian Mean (GLM) formalism, which has been developed by Andrews & McIntyre (1978b,c). More details on this theory can be found in [START_REF] Craik | Wave Interactions and Fluid Flows[END_REF], Grimshaw (1984), andBuhler (2009). Also in stellar physics, this formalism was first used by Lee (2013) and in later papers such as in Lee et al. (2016) for massive stars.

The GLM formalism can be understood as a generalization of the TEM formalism for finite amplitude waves but, even for small wave amplitudes, it provides an enlightening theoretical framework for understanding the wave/mean-flow interactions and the TEM. Let us first define the Lagrangian velocity following the wave displacement

v ξ ≡ v x + ξ ( x, t) , t , (51)
where ξ is the wave displacement. The main idea underlying the GLM is to adopt a mixed Eulerian-Lagragian approach in which mean Lagrangian velocity v ξ is defined by

v L ≡ v x + ξ ( x, t) , t . (52)
Here, the Eulerian average can be either a spacial, temporal, or ensemble average and as such the GLM is a very general theory. A schematic sketch of the two above-defined velocities is shown in Fig. 1 for a fluid particle trajectory.

To illustrate one of the advantage of the GLM, let us consider the conservation of specific entropy. It can be written such as

Ds Dt = Q . ( 53
)
Equation ( 53) can then be averaged using the classical Eulerian average so that it gives

Ds Dt -Q = -v • ∇s , (54)
which immediately shows that the mean entropy is not conserved (even if Q = 0) due to a secondorder wave related term. In sharp contrast, by adopting the GLM, one can easily show that

D L s L = Q L , ( 55
)
where the Lagrangian derivative is defined by

D L ≡ ∂ ∂t + v L • ∇ . ( 56
)
Equation ( 55) demonstrates that using the GLM permits the mean equation to conserve its conservative form. Then, for Q L = 0, the mean entropy (here s L ) is conserved. Note that, as for the TEM but here exactly, the only wave contributions to the mean equation appear in the momentum equation. In the limiting case X = 0 and Q = 0 this only contribution appears through the pseudo-momentum.

For adiabatic waves, this is a conserved quantity as shown by Andrews & McIntyre (1978b,c) thus demonstrating the non-acceleration theorem.

When compared to the TEM, this is an exact result and there is no residual term compared to Eq. ( 35). To understand this, it must be realized that the difference between the Lagrangian and Eulerian quantities are the Stokes corrections. For a scalar quantity, the Stokes correction is defined by

φ S ≡ φ L -φ . (57)
Therefore, the wave-related term in Eq. ( 56) can be written

v • ∇s = v S • ∇s + ∂ ∂t + v S • ∇ s S -Q S . (58)
The latter equation permits us to understand the difference between Eq. ( 53) and Eq. ( 55) as Stokes corrections are second-order quantities in term of wave amplitude. Furthermore, it also gives us an interpretation for the advective part of the divergence of the wave heat flux appearing in Eq. ( 30) within the TEM. This velocity can be understood as due to the Stokes correction. We also stress that in the case Q = 0 and if the wave is adiabatic one has

v † ⊥ = v L , (59)
which means that GLM and residual circulations are the same in the limiting case of Q = 0 and adiabatic waves.

4 Application to mixed modes in the low-mass evolved stars

In this section, based on the work of Belkacem et al. (2015a,b), we illustrate how the TEM can be used to estimate the transport of angular momentum by waves. This work focused on mixed modes but can also be easily generalized for progressive waves.

The problem of angular mometum reditribution in low-mass evolved stars

The stars and concluded that the core of subgiant stars spins up, while their envelope decelerates. Due to the contraction of their core and the expansion of their envelope, while such a result could appear as not being surprising, the problem is that local conservation of angular momentum would have produced a much more pronounced differential rotation. Therefore, a physical mechanism (which is to be identified) is needed for extracting angular momentum and thus explaining the observations. For more evolved stars, red giants, Mosser et al. (2012) and Gehan et al. (2018) analysed a sample of hundreds of red-giant stars observed by Kepler and found that, surprisingly, the mean core rotation rate decreases significantly during the red-giant phase. The actual observational picture is summarized by Fig. 2. All those observational constraints emphasize the need of angular momentum redistribution between the core and the envelope. However, current models of red-giant stars including angular momentum redistribution processes are unable to explain such low core rotation rates in subgiant and red giant stars. They are also unable to explain the deceleration of the core during the ascent of the red-giant branch. Consequently, the quest for a physical mechanism able to extract angular momentum from the core of evolved stars is currently under way. Meridional circulation and shear instabilities have been shown to be inefficient (e.g., Eggenberger et al. 2012, Marques et al. 2013). The same conclusion was reached by [START_REF] Cantiello | [END_REF] for a magnetic field generated through the Tayler-Spruit dynamo (Spruit 1999[START_REF] Spruit | [END_REF]. However, Fuller et al. (2019) investigate the magnetic Tayler instability and argue they were able to reproduce the nearly rigid rotation of main sequence stars and the core rotation rates of low-mass red giants. For internal gravity waves, it has been found that they are unable to explain the slow-down of red-giants (see Fuller et al. 2014 andPinc ¸on et al. 2017). In contrast, Pinc ¸on et al. ( 2017) had demonstrated that it is able to explain the differential rotation as observed in subgiants.

In this context, Belkacem et al. (2015a,b) investigated the ability of mixed modes to extract angular momentum within the TEM formalism. The aim was to estimate the influence of mixed modes on the angular momentum evolution compared to the structural effect of contraction of the core. To do so, they started from Eq. ( 49) and after meridional average, the equation governing angular momentum evolution can be written as

ρ d (r 2 Ω 0 ) dt = - 1 r 2 ∂ ∂r r 2 F waves ≡ J , (60)
where the symbol denotes the horizontal average (i.e. azimuthal and meridional). The wave flux is defined by

F waves = ρ v φ v r + 2 cos θ Ω 0 v θ s d s dr -1 , (61)
where = r sin θ, the prime denotes perturbations associated with the non-radial oscillations, so that v φ , v r , v θ are the azimuthal, radial, and meridional component of the wave velocity field, and s the wave Eulerian perturbation of entropy. We have also introduced the Lagrangian derivative d/dt = ∂/∂t + ṙ ∂/∂r.

Computing the wave field

To go further and to quantify the transport of angular momentum by mixed modes in the radiative region of evolved low-mass stars (i.e. subgiants and red giants), one has first to compute the wave field and thus to estimate the wave flux as given by Eq. ( 61). As these modes have a dual nature; they behave as acoustic modes in the upper layers and as gravity modes in the inner layers, several assumptions can then be adopted by using shellular rotation and focusing on the inner radiative regions. These restrictions allow the use of several approximations to describe the wave field:

1. The quasi-adiabatic approach: It consists in neglecting the difference between adiabatic and non-adiabatic eigenfunctions in the full wave equations. This approximation is valid when the local thermal time-scale is much longer than the modal period. This is the case in the radiative region of evolved low-mass stars.

2. The low-rotation limit: the modal period is assumed to be much shorter than the rotation period. This is justified by inferences of the rotation in the core of subgiants (Deheuvels et al. 2012(Deheuvels et al. , 2014) ) and red giants Mosser et al. (2012) using seismic constraints from Kepler.

3. The asymptotic limit: an asymptotic description for gravity modes (e.g., Dziembowski et al. 2001;Godart et al. 2009) is valid for mixed modes in the inner radiative region of subgiants and red giants (e.g., Goupil et al. 2013).

The next requirement is the determination of mode amplitudes. This is a fundamental point since it determines the amount of angular momentum transported by mixed modes. To do so, there are mainly two approaches. The first is based on a full non-adiabatic computation including a time-dependent treatment of convection. This procedure is time-consuming, however, and still suffers from theoretical uncertainties (see Dupret et al. 2009, Grosjean et al. 2014, for details). The second is based on recent CoRoT and Kepler observations that made it possible to establish scaling relations that provide mode amplitudes versus global stellar parameters (e.g., Mosser et al. 2012, Samadi et al. 2012). Belkacem et al. (2015b) followed the latter approach as it provides more reliable results.

Efficiency of angular momentum extraction

Using the simplifying assumptions as described in Sect. 4.2 and after quite involving calculations, Belkacem et al. (2015a,b) computed the divergence of the wave flux appearing in Eq. ( 61). It led to the conclusion that prograde modes (m < 0) extract angular momentum in the regions near the maximum of the buoyancy frequency and the maximum of the rotation rate, therefore slowing down the core. Conversely, the retrograde modes (m > 0) tend to spin up the core.

However, the net angular momentum flux is related to the asymmetry between prograde and retrograde modes. Therefore, it is worthwhile computing the net contribution of prograde and retrograde modes for a given angular degree. A simplified, but still accurate, expression is given by (see Belkacem et al. 2015a for details)

J( , -|m|) + J( , |m|) ≈ 2|m| 2 ρk 2 r a 2 ,|m| |ξ ,|m| r | 2 Ω 0 σ R N 2 σ R α , ( 62
)
where m is the azimuthal order, k r the radial wave number, a ,|m| the mode amplitude, ξ

,|m| r the radial eigen-function, σ R the pulsational frequency, and

α = - L 4πr 2 ρT ∇ ad ∇ -1 ds dr -1 (63)
with L the luminosity, T the temperature, ∇ and ∇ ad the temperature gradient and its adiabatic counterpart, and s the entropy.

Because the term α (see Eq. 63) is negative, Eq. ( 62) is also negative so that the sum of the contributions of prograde and retrograde mixed modes implies an extraction of angular momentum. Therefore, the collective effect of mixed modes is to decrease the mean angular momentum and thus to slow down the core rotation of the star. To go further and to assess the efficiency of the transport of angular momentum, four benchmark models have been chosen as depicted in Fig. 3 (right panel). For those models, it is useful to define two characteristic time-scales, namely;

• the timescale associated to the efficiency of the transport of angular momentum by mixed modes, defined as

T -1 m = J ρr 2 Ω 0 . (64)
• the timescale associated to the contraction of the star, defined as

T -1 c = - 1 ρr 4 Ω 0 ∂ ∂r ρr 4 Ω 0 ṙ ≈ - ṙ r , (65)
where ṙ = dr/dt.

Both timescales are shown in Fig. 3 (left panel). In all the radiative regions of models M0 and M1, the extraction of angular momentum remains negligible compared to the star contraction. In contrast, for the model M2, the timescale of angular momentum extraction is of the same order of magnitude as the contraction timescale and even lower in the hydrogen shell burning region for which the contraction is maximal. Therefore, one can conclude that the extraction of angular momentum by mixed modes is negligible in subgiants and early red giants while it becomes important in the hydrogen burning shell in stars higher on the red-giant branch. In such cases, mixed modes are able to counterbalance the spin-up due to the star contraction and can thus enforce a spin-down in those layers. The amount of angular momentum extracted by mixed modes increases with the evolutionary stage of the star. This effect is the result of several factors. As shown in Fig. 3, mode amplitudes increase from models M0 to M2 and thus more energy is available to transport angular momentum. Moreover, the buoyancy frequency also increases and so does the radial wave number in the gravity-mode cavity. As the energy exchanges between modes and the background is proportional to the radial wave number, the amount of angular momentum extracted also increases. Finally, the number of mixed modes between two p-dominated modes significantly increases between M0 and M2.

Concluding remarks

In this article, we tried to show that investigating the transport of angular momentum by waves in stars demands to address the more wide issue of the interaction between the waves and the meanflow. Indeed, the classical zonal or azimuthal average introduces wave-related terms in both the mean momentum and mean energy equations. Both are to be considered as they do not act independently because they are connected through meridional circulation. We therefore described the Transformed Eulerian formalism which ensures a proper modelling of the wave/mean-flow interactions. This formalism, adapted for small wave-amplitudes, allows to properly address the problem. This has been illustrated in the case of angular momentum transport by mixed modes in low-mass red giant stars. A more general theory, the Generalized Lagrangian Mean formalism has also been briefly introduced as it provides a nice framework for understanding the problem. Finally, we stress that considering only the mean momentum equation using a classical Eulerian average is a too naive approach which does not permit to grasp the physics of the problem. It can potentially leads to incorrect results because the coupling with the wave-related terms in the energy equation as ensured by meridional circulation is overlooked. Therefore, it would be desirable in the future to reanalyze previous works that only considered the mean momentum equation in the TEM and GLM frameworks.

  Figure 2: Rotation rate of the core as a function of the stellar radius. The open dots and open squares correspond to the RGB and clump stars, respectively, investigated by Mosser et al. (2012). The filled symbols indicate the subgiant stars studied by Deheuvels et al. (2012,2014). Figure from Deheuvels et al. (2015)

Figure 3 :

 3 Figure 3: Right panel: Evolutionary track on the HR diagram of a 1.3 M model showing the location of the three considered models in Sect. 4.3. Figure from Belkacem et al. (2015b). Left panel: Timescales versus normalised radius (i.e., normalised by the radius of the base of the convective envelope) for models M0, M1, and M2. The solid lines correspond to the timescale associated with the transport of angular momentum by mixed modes (see Eq. 64)) and the dotted lines correspond to the timescale associated with the contraction of the star (see Eq. 65). Figure from Belkacem et al.(2015b)

Note that by waves we denote both progressive waves and normal modes."How Much do we Trust Stellar Models?", held in Liège (Belgium), 10-12 September

Acknowledgements

I thank Marc-Antoine Dupret for his invitation to give a talk to the Liege conference, dedicated to Arlette Noels 75th birthday, and for his careful reading of the manuscript. I was honoured to spend several years in her research group and it was a very nice period when I have learnt a lot thanks to Arlette. I thus sincerely thanks Arlette Noels for our enlightening discussions and her constant support. I finally acknowledge financial support from the Programme National de Physique Stellaire (PNPS) of CNRS/INSU.


Ali Ak

Patrick Le-Callet

INVESTIGATING EPIPOLAR PLANE IMAGE REPRESENTATIONS FOR OBJECTIVE QUALITY EVALUATION OF LIGHT FIELD IMAGES

Keywords: light field, epipolar plane image, objective quality evaluation

With the ongoing advances in Light Field(LF) technology, research in LF acquisition, compression, processing has gained momentum. This increased the need for objective quality evaluation of LF content. Many processing algorithms are still optimized against peak signal to noise ratio(PSNR). Lately, several attempts have been made to improve objective quality evaluation such as extending 2D metrics to 4D LF domain. However, there is still a great room for improvement. In this paper, we experiment with existing 2D image quality metrics on the Epipolar Plane Image representations of LF content to reveal characteristics of LF related distortions. We discuss the challenges and suggest possible directions towards a LF image quality evaluation on EPI representations.

INTRODUCTION

The goal to produce ultra-realistic displays brought a lot of attention to Light Fields (LF) in the research community. Light Fields were introduced by Michael Faraday in 1846 on his lecture titled as "Thoughts on Ray Vibrations" [START_REF] Faraday | Liv. thoughts on ray-vibrations[END_REF]. In the last century, 7D plenoptic function is introduced. It describes the radiance received from every point to every point in 3D space at any given time at every wavelength. Later, it was simplified to 4D plenoptic function which allowed to easier computations. It is defined by 2 planes in space, parallel to each other. 4D function defines emitted light from any given point in the first plane to any point in the second plane. A comprehensive overview of the current literature for LF imaging has been done by Wu et al. [START_REF] Wu | Light field image processing: An overview[END_REF].

Light Field content can be represented in many different ways. In Fig. 1 sub-aperture views and a corresponding Epipolar Plane Image representation is presented. Subaperture view representation is the most common one where each view is an image from a slightly distant point of view depending on the baseline. The baseline is defined as the distance between each view. In this representation, we can observe the spatial information over each view. EPI representations, on the other hand, reveals information about the angular domain. It describes the shift of pixel information From acquisition to the display of the LF content, there are different processing steps. Most of these processing tools add a different type of distortion to the processed LF image. Captured LF data needs to be coded, transmitted, and decoded at the display device. This process also results in unwanted distortions which lower perceived quality. Quality evaluation metrics are necessary to evaluate the performance of the developed tools for each stage.

Objective image quality assessment for 2D images has attracted considerable attention over the last several decades [START_REF] De Angelis | Image quality assessment: an overview and some metrological considerations[END_REF]. Quality evaluation of videos and sparse multi-view content also attracted attention [START_REF] Chikkerur | Objective video quality assessment methods: A classification, review, and performance comparison[END_REF]. There have been attempts to create a metric for LF image quality assessment [START_REF] Roopak | Supermultiview content with high angular resolution: 3D quality assessment on horizontal-parallax lightfield display[END_REF]. Authors combined estimated spatial quality with angular quality to provide a score.

In this paper, we experiment with some of the existing 2D quality evaluation metrics on EPI representations. Rest of the paper is organized as follows: In Section 3, we investigate existing objective image quality metrics and select a set of metrics to experiment on EPI representations. In Section 2, we investigate existing LF image quality datasets and LF related distortions. In Section 4, we share our experiments and discuss the acquired results. Finally in Section 5 we conclude with our findings.

LIGHT FIELD RELATED DISTORTIONS

Different stages in LF visual chain may introduce a different type of distortion [START_REF] Wu | Light field image processing: An overview[END_REF]. Even though there is numerous type of distortions in common with 2D images or videos, there exist some distortions specific to LF content as well. Recently, very few LF image quality datasets have been published. Some of them consider common distortions [START_REF] Roopak | Supermultiview content with high angular resolution: 3D quality assessment on horizontal-parallax lightfield display[END_REF] such as gaussian blur, additive noise, etc. Some publicly available datasets consider LF specific distortions [6][7][8] such as compression, view synthesis related distortions. There is still a need for datasets who offers acquisition and display related distortions with the corresponding subjective scores.

In Free Viewpoint Video dataset(FVV) [START_REF] Sandi-Stankovi | Free viewpoint video quality assessment based on morphological multiscale metrics[END_REF], seven different compression algorithms have been used to code depth maps. Coded depth maps are used to synthesize intermediate views to generate horizontal parallax videos with 50 frames. Generated videos are mirrored to 100 frames and then evaluated subjectively. Even though FVV dataset consists only horizontal parallax, it has high angular resolution and LF related distortions. In Fig. 2, an EPI slice from the "balloons" image is presented. 7 mid-level distortions are also presented with different annotations. Extracted edge maps from each distorted image are also displayed next to them. Each coding algorithm introduces a different type of distortion to the EPI image. Some algorithms produce more subtle distortions than other algorithms which result in higher perceived quality in the subjective experiment. In the dataset by Adhikarla et al. [START_REF] Vamsi Kiran Adhikarla | Towards a quality metric for dense light fields[END_REF], LF resolution is chosen as 960 × 720 × 101. Several types of distortions are considered, such as transmission, reconstruction, and display related distortions. As in FVV dataset, there is only horizontal parallax but with a high angular resolution. In Fig. 3, an EPI slice with 4 different reconstruction distortions along its nondistorted version is presented. All of the presented distortions are highly visible on EPIs

OBJECTIVE IMAGE QUALITY EVALUATION

Numerous metrics have been proposed to assess the quality of images and videos over the last several decades [START_REF] De Angelis | Image quality assessment: an overview and some metrological considerations[END_REF]. On the other hand, PSNR is still used by averaging PSNR on subaperture images to evaluate the performance of many LF processing tools.

This paper focuses on full-reference image quality assessment (FR-IQA), where the reference image exists along with the evaluated distorted image. There are mainly two strategies for FR-IQA model design. One is using a bottom-up approach where different stages in the Human Visual System(HVS) is modeled separately to evaluate the image quality. These stages include visual masking, contrast sensitivities, and just noticeable differences. This results in computationally complex models. The second approach is called top-down, which tries to estimate the whole HVS behaviour with broad mathematical assumptions. This provides a computationally simple evaluation procedure while sacrificing accuracy. In our work, we focused on existing top-down approaches, and those utilize the structural information on the images.

Non-shift edge-based ratio(NSER) [START_REF] Zhang | Non-shift edge based ratio (nser): An image quality assessment metric based on early vision features[END_REF] evaluates the im-age quality by the edge differences between the reference and distorted images. Natural image counter evaluation(NICE) [START_REF] Rouse | Natural image utility assessment using image contours[END_REF] operates similarly, where the comparison is made between morphologically dilated edge maps. Gradient magnitude similarity deviation (GMSD) works with standard deviation of the image gradients to evaluate the image quality. Morphological Wavelet Peak Signal-to-Noise Ratio (MW-PSNR) [START_REF] Sandi-Stankovi | Dibr-synthesized image quality assessment based on morphological multi-scale approach[END_REF] calculates the image quality based on morphological wavelet decomposition. PSNR is calculated between each band of decomposed reference and distorted image pairs. GMSD works on directional image gradients to produce the final quality score. Prewitt filters along horizontal and vertical directions are used to convolve the image. They are defined as:

h x =   1/3 0 -1/3 1/3 0 -1/3 1/3 0 -1/3   , h y =   1/3 1/3 1/3 0 0 0 -1/3 -1/3 -1/3  
Both reference(r) and distorted(d) images are convolved with filters above. Gradient magnitudes at each pixel location for both images are calculated as follows:

m r (i) = (r h x ) 2 (i) + (r h y ) 2 (i) (1)
m d (i) = (d h x ) 2 (i) + (d h y ) 2 (i) (2)
These gradient magnitudes for each map then used to calculate gradient magnitude similarity(GMS) map as follows:

GM S(i) = 2m r (i)m d (i) (m r ) 2 (i) + (m d ) 2 (i) (3)
Standard deviation pooling is used on GMS to calculate the final quality score. Convolving the input images with Prewitt filters reveals the structural information on the image. This explains our motivation to choose GMSD to experiment.

NICE is slightly simpler, and it relies on edge detection algorithms. Both reference and distorted images are inputted to an edge-detection algorithm. Extracted edge maps are then dilated with a plus-sign kernel. Non-zero elements in the XOR maps between the dilated image pair is then used to calculate the final quality score.

EXPERIMENTS

Non-Uniform Distortion Distribution on EPI Slices

In Fig. 4, quality evaluation of 3 different distortions on the same LF image is plotted. NICE [START_REF] Rouse | Natural image utility assessment using image contours[END_REF] metric is used as a quality score. A lower number represents better quality. We have observed that the same EPI slices contribute to lower quality despite the fact that different source of distortions is introduced. EPI slices which belong to occluded regions contribute more to lower the perceived quality. This suggests that averaging quality scores across all EPI slices may not be the optimal solution for pooling. Additionally, lower quality EPI slices might be used as a selective process to identify distorted regions on the spatial domain.

Edge Detection on EPIs

PCC SROCC Canny 0.7145 0.6552 Sobel 0.6457 0.5951 Prewitt 0.6450 0.6038 Table 1. NICE [START_REF] Rouse | Natural image utility assessment using image contours[END_REF] metric performance on the dataset by Adhikarla et al [START_REF] Vamsi Kiran Adhikarla | Towards a quality metric for dense light fields[END_REF]. Metric is incorporated with different edge detector algorithms [START_REF] Canny | A computational approach to edge detection[END_REF] [START_REF] Kanopoulos | Design of an image edge detection filter using the sobel operator[END_REF].

Few edge detectors are commonly used in literature such as Canny [START_REF] Canny | A computational approach to edge detection[END_REF], Sobel [START_REF] Kanopoulos | Design of an image edge detection filter using the sobel operator[END_REF]. Even though most of the methods have been reviewed extensively in the literature [START_REF] Shrivakshan | A comparison of various edge detection techniques used in image processing[END_REF], their performance on EPI for quality evaluation has not been investigated yet. We have experimented with three different edge detector algorithms and evaluated their performance on revealing distortions on EPIs. We incorporated different edge detectors on NICE [START_REF] Rouse | Natural image utility assessment using image contours[END_REF] metric. In table 1, the effect of different edge detectors on the performance of NICE metric is reported. For evaluation, 2 different source images have been chosen from the dataset by Adhikarla et al [START_REF] Vamsi Kiran Adhikarla | Towards a quality metric for dense light fields[END_REF], which resulted in 800 comparisons. Even though it does not reflect the full performance of the metric on the dataset, it is enough comparisons to evaluate the edge-detector performances for NICE algorithm. In table 2, Pearson correlation coefficients with subjective scores are presented for MW-PSNR [START_REF] Sandi-Stankovi | Dibr-synthesized image quality assessment based on morphological multi-scale approach[END_REF], GMSD [START_REF] Xue | Gradient magnitude similarity deviation: A highly efficient perceptual image quality index[END_REF] and NICE [START_REF] Rouse | Natural image utility assessment using image contours[END_REF] metrics on FVV dataset [START_REF] Sandi-Stankovi | Free viewpoint video quality assessment based on morphological multiscale metrics[END_REF]. All metrics have been tested on the same dataset with both LF representations. We observed that GMSD and NICE provide higher correlation when it is used with EPI images as input. On the other hand, MW-PSNR provides a lower correlation with subjective opinions when used with EPI inputs. Low spatial resolution of EPIs provides only 2 scales for morphological wavelet decomposition, which is utilized by MW-PSNR metric. This results in low accuracy in terms of prediction for the MW-PSNR metric. Gradient provides structural information regarding the image, and in EPIs this information reveals the distortions in the angular domain. This results in higher accuracy for GMSD and NICE, which relies on filtering the image with Prewitt and Canny kernels. Another experimentation has been made on individual distortions available on the FVV dataset. Some of the available distortions are more subtle on the EPIs. Both metrics, MW-PSNR and GMSD, has been used to evaluate the quality of individual distortions with both EPI slices and GMSD as an input. Pearson correlation coefficients of the evaluations are reported in table 3. Low visibility of 3D-HEVC coding distortions on the EPI slices results in lower accuracy of the quality evaluation for both metrics when used with EPI slices. Multi-View video coding distortions is very easy for both of the metrics to predict, and there are no statistically significant differences between different input types. We observed that HEVC test model coding and JPEG2000 coding results in more visible distortions on many of the EPI slices. This provides a better correlation with subjective scores for both metrics on EPIs.

Individual Distortions

MW

CONCLUSION

We have run several experiments to evaluate existing metric performances on EPI representations. In our experiments, we have adopted two different metrics (NICE, GMSD) which rely on structural information such as gradient and extracted edge information. Although both metrics are not designed for EPI domain, we have observed that their performances increased compared to sub-aperture view performances.

On the other hand, metrics which utilize multi-scale decomposition cannot fully benefit from the low spatial resolution of EPIs. Multi-scale decomposition algorithms are designed to separate the image into several components based on spatial frequency and orientation ranges. This separation does not provide any information in the angular domain because an EPI slice is a line across an angular direction and does not contain any information regarding its neighborhood pixels in the opposite direction. In order to benefit from EPI representations, structural feature extraction, which is not based on natural image statistics is required. Currently, simple edge extraction algorithms do not provide a satisfactory correlation with subjective evaluations. Experimental results suggest that with a carefully designed feature extraction scheme, EPI representations might be useful for identifying numerous type of angular domain distortions.

Additionally, when both horizontal and vertical parallax is taken into consideration, the same pixel values need to be evaluated twice. This doubles the computational complexity. A selective mechanism on the spatial domain to locate the occluded areas can be beneficial. This may allow identifying EPI slices which contribute most to distortions on the angular domain.

Fig. 1 .

 1 Fig. 1. 2 different Light Field image representation. On the left, each sub-aperture views. On the right, an EPI slice

Fig. 2 .

 2 Fig. 2. 7 different distortions on a single EPI slice from the FVV dataset[6] and their corresponding edge maps. a)Reference, b)3D-HEVC, c)Multi-View Video, d)HEVC Test Model, e)JPEG200, f)Lossless edge depth coding, g) Color channel corrolation based coding, h)Z-LAR-RP

Fig. 3 .

 3 Fig. 3. 4 of 7 different distortions on a single EPI slice from the dataset by Adhikarla et al.[9] and their corresponding edge maps a)Reference, b)Depth map interpolation, c)Optical flow interpolation, d)Linear interpolation, e)Nearest neighborhood interpolation

Fig. 4 .

 4 Fig. 4. Each line in the plot represents a different type of distortions on the image("balloons" from FVV dataset). Every point in the X axis represents a single EPI slice. Quality scores are calculated based on NICE.

Table 2 .

 2 Pearson Correlations of MW-PSNR, GMSD, NICE metrics on EPIs and Sub-Aperture views. FVV dataset has been used for evaluation.

	4.3. Sub-aperture view vs EPI
		EPI	View
	MW-PSNR 0.7698 0.7921
	GMSD	0.7410 0.6715
	NICE	0.5122 0.4310

Table 3 .

 3 Pearson Correlations of MW-PSNR and GMSD metrics on EPIs and Sub-Aperture views for individual distortions available in FVV dataset.

		-PSNR	GMSD
		EPI	View	EPI	View
	3D-HEVC	0.4665 0.7137 0.3644 0.5110
	Multi-View Video	0.9549 0.9484 0.8649 0.8746
	HEVC Test Model	0.9584 0.8662 0.7974 0.7059
	JPEG2000	0.8268 0.7687 0.7109 0.6643
	Lossless edge depth 0.8016 0.8635 0.8288 0.8186
	Color corr. based	0.7940 0.9099 0.9058 0.7315
	Z-LAR-RP	0.7195 0.8147 0.7811 0.7317

Acknowledgment

This project has received funding from the European Unions Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Grant Agreement No. 765911 (Re-alVision)


Ali Ak
email: ali.ak@univ-nantes.fr

Patrick Le-Callet
email: patrick.lecallet@univ-nantes.fr

Towards Perceptually Plausible Training of Image Restoration Neural Networks

Keywords: visual perception, neural networks, just noticeable difference

Learning-based black-box approaches have proven to be successful at several tasks in image and video processing domain. Many of these approaches depend on gradient-descent and back-propagation algorithms which requires to calculate the gradient of the loss function. However, many of the visual metrics are not differentiable, and despite their superior accuracy, they cannot be used to train neural networks for imaging tasks. Most of the image restoration neural networks rely on mean squared error to train. In this paper, we investigate visual system based metrics in order to provide perceptual loss functions that can replace mean squared error for gradient descentbased algorithms. We also share our preliminary results on the proposed approach.

I. INTRODUCTION

Visual metrics has numerous use cases in visual processing domain. They play an important role in the development, evaluation, and optimization of many visual processing algorithms. There are various approaches to develop visual metrics. While some of the metrics focus on signal driven calculations [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF] [8] [START_REF] Aydin | Extending quality metrics to full luminance range images[END_REF], some focus on modeling the visual system [START_REF] Mantiuk | HDR-VDP-2: A calibrated visual metric for visibility and quality predictions in all luminance conditions[END_REF] [START_REF] Daly | The visible differences predictor: an algorithm for the assessment of image fidelity[END_REF]. Metrics which relies on signal driven calculations model the quality perception as a continuous function. On the other hand, Visual Model based metrics, such as VDP [START_REF] Daly | The visible differences predictor: an algorithm for the assessment of image fidelity[END_REF] and HDR-VDP [START_REF] Mantiuk | HDR-VDP-2: A calibrated visual metric for visibility and quality predictions in all luminance conditions[END_REF], can predict the perceptual quality of the images more accurately and tuned on the Just Noticeable Differences around near threshold values. Although they are more accurate, they have a high computational complexity since they are derived from different components of Human Visual System (HVS) where the data is collected from a set of psychophysical measurements. Additionally, this complexity results in nondifferentiable models which prevents them to be used in many visual processing applications.

Visual model based metrics can estimate the probability of the detection of differences between a pair of images for an observer. They are generally tuned towards the near-threshold Just Noticeable Differences. This provides an accurate judgment of perception by providing information regarding statistically significant differences between a pair of stimuli.

HVS is not sensitive to the low amount of differences on a pair of stimuli. Several factors affect the perception of a difference. Under certain conditions, the minimum amount of difference required to be perceived is called the Just Noticeable Difference(JND). However, most of the existing visual metrics provide a continuous quality score for a varying amount of distortions. Visual model based metrics such as HDR-VDP provides more perceptually plausible predictions by focusing on just noticeable differences.

Neural networks which are trained for image restoration tasks such as denoising, compression, super-resolution, view synthesis uses loss functions during the training. Most of the current approaches rely on MSE. This results in blurry highfrequency details, and noisy flat regions depending on the task.

There are several attempts to tackle the lack of perceptual training of image restoration tasks. Zhao et al. provide an indepth analysis of loss functions in image restoration tasks [START_REF] Zhao | Loss Functions for Image Restoration With Neural Networks[END_REF]. In their paper, authors also provide a differentiable version of SSIM as an alternative loss function which also fails to provide perceptually plausible results.

Lately, Perceptual loss has been introduced by Johnson et al [START_REF] Johnson | Perceptual losses for real-time style transfer and super-resolution[END_REF]. They used layer activations of deep neural networks that have been trained for image classification tasks on millions of images such as VGG [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF] and ResNet [START_REF] Wu | A New Approach to Compute CNNs for Extremely Large Images[END_REF]. Activations of different layers have been used as a loss function in image restoration tasks which provides perceptually more plausible results compared to MSE.

Goodfellow et al. proposed a unique network scheme called Generative Adversarial Networks(GANs) [START_REF] Goodfellow | Generative adversarial nets[END_REF]. GANs are deep neural network architectures contains two individual networks, computing with each other. Generative network produces an image to fool the discriminator that it is a real image while discriminator tries to predict the inputted image is from the real distribution or a generated one. This training scheme allows for capturing any distribution of data. For image restoration tasks, GANs combined with traditional loss functions provide perceptually plausible results compared to former alternatives. To keep output image similar to input, GANs need to incorporate with a loss function such as perceptual loss or MSE.

Despite the success of learning-based black-box approaches at image restoration tasks, there is still a need for perceptually plausible training schemes. Complexity and non-differentiable 978-1-7281-3975-3/19/$31.00 ©2019 European Union nature of the HVS based visual metrics prevent them from being used as a loss function for the training of neural networks. In this paper, we will investigate different modules from HVS based visibility metric, HDR-VDP [START_REF] Mantiuk | HDR-VDP-2: A calibrated visual metric for visibility and quality predictions in all luminance conditions[END_REF], and how we can replicate them with neural networks. Later, we will share our preliminary results and conclude with our findings.

II. HDR-VDP OVERVIEW

While an arbitrary amount of distortion is perceivable under specific viewing conditions, it might not be noticeable under another. There are many factors involved in the perception of distortions [START_REF] Foley | Human luminance pattern-vision mechanisms: masking experiments require a new model[END_REF], such as viewing density, display luminance, background luminance, and individual observer differences. Many of these factors have been well explored and modeled to some extent in the literature.

We will investigate some of these factors and how they are modeled in HDR-VDP [START_REF] Mantiuk | HDR-VDP-2: A calibrated visual metric for visibility and quality predictions in all luminance conditions[END_REF]. HDR-VDP is a comprehensive model of the human visual system. It is built with several modules, such as Intra-ocular Light Scattering, Photoreceptor Spectral Sensitivity, Luminance Masking, Multi-scale Decomposition, Neural Contrast Sensitivity, Masking, Multiscale Reconstruction. It is designed to estimate perceivable differences for near-threshold distortions.

Working with gamma-corrected images as input and assuming that the pixel values are scaled perceptually uniform may mislead the perceptual evaluation of produced images. The following formula is used to convert the pixel intensity values into absolute luminance values.

Luma = 0.2126 × R + 0.7152 × G + 0.0722 × B L = (L P eak -L Black ) × Luma (g) + L Black
Where R, G, B values are first converted into grayscale (Luma) and according to display device's Peak Luminance (L P eak ), Black level Luminance (L Black ) and gamma (g), final Luminance values can be calculated for each pixel in terms of cd/m 2 .

Optical and Retinal Pathway of the HVS modeled as a combination of separate modules in HDR-VDP [START_REF] Mantiuk | HDR-VDP-2: A calibrated visual metric for visibility and quality predictions in all luminance conditions[END_REF]. The first module is intra-ocular light scatter in which the scattered light is passing through the cornea, lens, and eye chamber is modeled [START_REF] Ritschel | Temporal Glare: Real-Time Dynamic Simulation of the Scattering in the Human Eye[END_REF]. This process results in a loss of contrast of the perceived image(scene).

After the intra-ocular light scatter module, photo-receptor spectral sensitivities are modeled as non-linear functions. It describes the probability of L, M, S cones, and rods on the retina sensing each particular wavelength.

Photo-receptor spectral sensitivity module is followed by Luminance masking in HDR-VDP [START_REF] Mantiuk | HDR-VDP-2: A calibrated visual metric for visibility and quality predictions in all luminance conditions[END_REF]. Luminance masking accounts for the local luminance adaptation of photo-receptors. Most of the visual metrics assume a global adaptation of global adaptation to luminance values, while photo-receptors has a highly non-linear response to different luminance levels [START_REF] Aydin | Extending quality metrics to full luminance range images[END_REF]. This adaptational sensitivity to luminance is captured in the contrast sensitivity function, and it is derived from the measurements by Mantiuk et al [START_REF] Mantiuk | Predicting visible differences in high dynamic range images: model and its calibration[END_REF]. Afterward, achromatic responses are calculated as the sum of L, M cones, and rod responses acquired from Luminance adaptation module.

Psychophysical studies suggest that HVS responses to each spatial frequency and orientation with a different sensitivity [START_REF] Foley | Human luminance pattern-vision mechanisms: masking experiments require a new model[END_REF]. Many visual metrics employ one of several multi-scale approaches. In Daly's VDP Cortex Transform [START_REF] Watson | The cortex transform: Rapid computation of simulated neural images[END_REF] is used while Steerable Pyramid [START_REF] Simoncelli | The steerable pyramid: a flexible architecture for multi-scale derivative computation[END_REF] is used in HDR-VDP [START_REF] Mantiuk | HDR-VDP-2: A calibrated visual metric for visibility and quality predictions in all luminance conditions[END_REF]. After multi-scale decomposition, neural CSF function is applied to each spatial frequency and orientation bands individually. Neural CSF is signal independent and can be derived from CSF measurements which are acquired from psychophysical experiments. Before applying, optical and luminance dependent components are excluded (by division) from fitted CSF function since they are already modeled in earlier stages. Contrast masking is also applied in three different parts as self masking and masking across orientations and masking due to two neighboring frequencies. Afterward, the psychometric function is applied, and probability summation is done as summation across all frequencies.

This visual model provides a probability of detection map as well as taking the max value over an image, produces a probability of detection score. Due to the lack of benefit in our task, we skipped the quality score branch of the HDR-VDP. Authors also provided an extension for the quality score.

III. OPTICAL AND RETINAL PATHWAY NEURAL NETWORK IMPLEMENTATION

Optical and Retinal Pathway modules are jointly modeled and trained with a combined loss function of each three individual modules. General scheme of the model is summarized in Fig. 1. Each network receives its inputs from the previous module as it has been implemented in HDR-VDP. An equally weighted loss function has been used during the training where the ground truth data has been acquired from HDR-VDP. In each of the Optical and Retinal Pathway neural networks, reference images and test images are inputted separately as it is implemented in HDR-VDP. For the following modules, we have used an auto-encoder neural network architecture since we want each module to represent their corresponding white-box outputs accurately. Using auto-encoder for each module in the Optical and Retinal Pathway allowed us to keep the input and output size the same. Additionally, it allowed us to train all three networks with a combined loss function, and by restricting the networks on individual tasks, it allowed us to subdivide the problem. In each of the following subsections, we will provide our incentives for the neural network structures and provide details about the individual neural network structures.

A. Intra-Ocular Light Scatter

In HDR-VDP, this module is modeled as a modulated transfer function. Intra-ocular light scatter module works on spectral values, and in an ideal model, there should be a small wavelength dependency, however, most studies show that there is little to none wavelength dependency in intra-ocular light scatter [START_REF] Whitaker | Light scatter in the normal young, elderly, and cataractous eye demonstrates little wavelength dependency[END_REF]. This allows using the same MTF for each input radiance map. In HDR-VDP, authors used the psychophysical measurements to fit the function. Then, this function is used to create MTF kernels for the required size.

MTF kernel is multiplied by the input spectral values in the frequency domain. As it can be observed in Fig. 2, we have used six convolution layers in this network to capture the same receptive field with the MTF kernel in HDR-VDP. We have used an auto-encoder structure and used bilinear interpolation instead of deconvolution in the decoder network to prevent checkerboard artifacts [START_REF] Odena | Deconvolution and Checkerboard Artifacts[END_REF]. ReLu activation functions have been used for non-linearity. This helped us to achieve a computationally more efficient and differentiable replicate of the module with similar accuracy.

B. Photo-receptor Spectral Sensitivities

In HDR-VDP, Photo-receptor Spectral sensitivities are modeled as non-linear functions for each photo-receptor type separately over the wavelength range. Photo-receptor sensitivities are acquired via multiplying the fitted functions with each channel.

As it can be observed from Fig. 3, in our neural network implementation, the output of the Intra-Ocular Light Scatter Network is inputted directly to this network as it is implemented in HDR-VDP. We have modeled this module as a combination of convolution layers and ReLu activation functions (for non-linearity) in the network. Combination of ReLus provided us a fair estimation of the non-linearity of spectral sensitivity curves.

C. Luminance Masking

To learn the non-linear adaptation responses, we have used a similar auto-encoder scheme again. Instead of ReLu, we have used SeLu activation functions as non-linearity after each convolution layers. We have adapted SeLu slope by adjusting α value in the SeLu function according to the data acquired from the HDR-VDP [START_REF] Mantiuk | HDR-VDP-2: A calibrated visual metric for visibility and quality predictions in all luminance conditions[END_REF]. We have combined Luminance Masking and Achromatic Responses modules together in our implementation since Achromatic responses are implemented as an equally weighted summation of luminance masking output channels in HDR-VDP [START_REF] Mantiuk | HDR-VDP-2: A calibrated visual metric for visibility and quality predictions in all luminance conditions[END_REF].

IV. PROBABILITY OF DETECTION NETWORK

Instead of predicting a probability of detection map, in our first implementation, we chose to acquire just a probability of detection score. This allowed us to try and train image restoration networks by just replacing common loss functions.

Currently, we are working on implementing Steerable Pyramid spatial frequency scales and orientation filters. Instead of training the network freely, we aim to use the existing knowledge from the white-box approach, which we already used as ground truth.

To train our probability of detection neural network, we have used MCL-JCI dataset [START_REF] Jin | Statistical Study on Perceived JPEG Image Quality via MCL-JCI Dataset Construction and Analysis[END_REF]. MCL-JCI dataset contains 50 source images compressed where each of them compressed by JPEG algorithm with 100 different Q levels. Subjective tests are conducted on the dataset where subjects are asked if they notice any distortion a reference image and a distorted version of it. JND staircase functions are acquired via subjective experiment. After experimenting with the dataset, to focus on near-threshold distortions, we have used only the first JND step information. We have extracted 50 × 100 image pairs to compare with 1920 × 1080px size. We have assigned each image pair with a score between [0,1] depending on the agreement of the population about the visibility of the distortion.

We have used a Siamese neural network architecture in our implementation. Siamese networks have several advantages. They don't require as much as data compared to other neural networks architectures. They learn similarity functions, so instead of a classification task, they are trained on discriminating which directly related to our task.

Achromatic responses of both reference and distorted images, which acquired from the previous neural networks, are inputted to the probability of detection network. The network is trained with data generated from MCL-JCI dataset. Due to large size of input images (1920 × 1080px), we used a narrow architecture. Pairwise distance between the output vector of input images is trained against the visibility score extracted from MCL-JCI dataset.

V. TRAINING DETAILS

Optical and Retinal Pathway implementation contains three similar auto-encoder neural networks. Each auto-encoder pro- duces an intermediate result. Their ground truth data is acquired from HDR-VDP [START_REF] Mantiuk | HDR-VDP-2: A calibrated visual metric for visibility and quality predictions in all luminance conditions[END_REF]. We have created a dataset by using HDR-VDP. We have created a dataset of 300000 images with 128×128px size, and varying background luminance, pixel per degree. We have used 0.0005 learning rate during the training and adaptive momentum optimizer(Adam) for the training. We have used 16 as batch size. We have trained the network over 500 epochs.

After training of the Optical and Retinal Pathway, we have used the network in evaluation mode to process the input reference and distorted images for the training of probability of detection network. We have used the same learning rate of 0.0005 and Adam optimizer for the training. Due to the high resolution of input images, we have used 8 batch size during training. We have used Pytorch library for the implementation.

VI. PRELIMINARY RESULTS

We have processed 100 images (1920 × 1080px sized) with both HDR-VDP 2 original implementation and our proposed method. On average we have acquired a 120 times speed improvement with the same computer specifications.

In Fig. 5, we can see the visibility score predictions of HDR-VDP and proposed approach. The horizontal axis represents the source image with 100 Q levels of JPEG distortions. From left to right, the Q parameter of the JPEG compression algorithm increases. On the vertical axis, quality score predictions are displayed. Blue line is subjective scores from MCL-JCI dataset, the red line is our prediction, and the yellow line is HDR-VDP 2 predictions. As can be observed, our predictions are better aligned with subjective opinions.

After confirming our perceived quality predictions, we have used our model to train image restoration neural network algorithms. In Fig. 6, We have shared our preliminary results for a denoising network on MNIST dataset. Input images with noise, ground truth images, the output of a neural network trained with our model as an objective function and the same network trained to minimize MSE has been provided. The neural network has an auto-encoder structure, and while trained with our proposed metric and MSE, all the hyperparameters are kept same. The network trained with proposed metric outputs smoother results. Solid white regions have

VII. CONCLUSION

Developing a perceptually plausible visual metric is a challenging task. Due to the complex nature of HVS, proposed metrics have mathematically complex and non-differentiable models. This prevents to use existing visual metrics to train neural networks with gradient-descent based optimization methods. In this paper, we examined the shortcomings of existing loss functions and discussed some of the proposed solutions to overcome these shortcomings. Then, we investigated several critical elements of existing visual system based models and presented our approach to implement them in a neural network scheme. We shared our preliminary results on training a denoising auto-encoder with the current state of the model as well as JND based dataset MCL-JCI. As future work, we plan to integrate a visual difference probability map into our visual model network and improve our training scheme of image restoration networks. Ultimately, we plan to achieve a perceptually plausible GAN-like training scheme. We also work on using the proposed metric as a loss function in different image restoration tasks, such as super-resolution, compression, etc.

Fig. 1 .

 1 Fig. 1. HDR-VDP [2] box diagram with our proposed scheme. Red dashed rectangles represent individual networks in our approach. Each colored box represents different modules from HDR-VDP [2]. Optical and Retinal Pathway of HDR-VDP is color coded with darker orange.

Fig. 2 .

 2 Fig. 2. Intra-Ocular Light Scatter Neural Network box diagram.

Fig. 3 .

 3 Fig. 3. Photo-receptor Spectral Sensitivity Neural Network box diagram.

Fig. 4 .

 4 Fig. 4. Probability of detection neural network architecture. Both reference and distorted images go trough the network separately and network is trained against the pairwise distance between the extracted vectors.

Fig. 5 .

 5 Fig. 5. Predicted Visibility Scores from proposed metric vs HDR-VDP on a source image from MCL-JCI dataset.

Fig. 6 .

 6 Fig.6. Denoising neural network trained on MNIST dataset with two different objective functions. On the 3rd columns, outputs of the neural network trained with proposed metric is displayed, and on the 4th column, outputs of the same neural network trained with MSE is displayed.

ACKNOWLEDGMENT

This project has received funding from the European Unions Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Grant Agreement No. 765911 (Re-alVision)