S Clain

D Rochette

R Touzani

A multislope MUSCL method on unstructured meshes applied to compressible Euler equations for swirling ows

Keywords: Axisymmetric compressible Euler equations, Swirling ow, Finite volumes, MUSCL method, Unstructured mesh

A nite volume method for the numerical solution of axisymmetric inviscid swirling ows is presented. The governing equations of the ow are the axisymmetric compressible Euler equations including swirl (or tangential) velocity. A rst-order scheme is introduced. In this one, convective uxes at cell interfaces are evaluated by the Rusanov or the HLLC numerical ux and geometric source terms are discretized by the explicit Euler method. Extension to the second-order space approximation using a multislope MUSCL method is derived. A stationary solution of the uid ow following the radial direction has been established with a zero and non-zero tangential velocity. Numerical and exact solutions are compared for the Riemann problem. Effectiveness of the multislope MUSCL scheme is demonstrated for strongly shocked axially symmetric ows as the forward-facing step and the spherical bubble compression problems.

1 Introduction Axisymmetric Euler system using cylindrical coordinates is used in numerous applications such as axisymmetric ows in a nozzle [START_REF] Jacob | Transient, hypervelocity ow in an axisymmetric nozzle[END_REF], supersonic jets [START_REF] Leschziner | Computation of strongly swirling axisymmetric free jets[END_REF], turbomachine modelling [START_REF] Wu | A general theory of three-dimensional ow in subsonic and supersonic turbomachine in radial, axial and mixed ow types[END_REF][START_REF] Ch | A nite element method for the axisymmetric ow computation in a turbomachine[END_REF]. More recently, inductive plasma ows are modelled with the axisymmetric Euler formulation taking into account the tangential velocity, the so-called swirling ow [START_REF] Vanden Abeele | An ecient computational model for inductive plasma ows[END_REF][START_REF] Xue | Extended-led electromagnetic model for the inductively coupled plasma[END_REF].

From a numerical point of view, the nite volume method [79] is a popular technique to compute numerical approximations of the Euler system solution for axisymmetric geometries. A particular issue concerns the choice of the variables to conserve. In a rst approach, the mean value approximation of any generic function v on the cell C i is performed by using the classical average

v i ≈ C i v drdz C i drdz
where we employ the measure (metric) dr dz [START_REF] Cumming | Supersonic, turbulent ow computation and drag optimization for axisymmetric afterbodies[END_REF]. A second approach consists in computing the mean value using the measure r dr dz [START_REF] Guardone | Finite element/volume solution to axisymmetric conservation laws[END_REF]:

v i ≈ C i v r drdz C i r drdz
which leads to a better formulation. Indeed, v i corresponds to the mean value of v in the three-dimensional context, i.e. the mean value on the axisymmetric torus with section C i . On the other hand, the formulation simplies the boundary condition at r = 0 since a null ux value naturally derives from the ux integration and no additional constraint is required on the symmetry axis [START_REF] Guardone | Finite element/volume solution to axisymmetric conservation laws[END_REF].

To provide an approximation of the solution of the axisymmetric Euler system, we use a fractional step technique where we split the formulation into a conservative homogeneous equation and the source term. To solve the conservative part, the traditional technique consists in using a rst-order solver (Roe, HLLC or Riemann solver) combined with a second-order reconstruction such as the MUSCL method [START_REF] Van Leer | Towards the ultimate conservative dierence scheme. V. A secondorder sequel to Godunov's method[END_REF][START_REF] Hubbard | Multidimensional Slope Limiters for MUSCL-Type Finite Volume Schemes on Unstructured Grids[END_REF] to improve accuracy. The classical MUSCL technique uses a piecewise linear reconstruction with a slope limiting procedure to ensure the L ∞ stability. Then two new approximations are computed on both sides of each edge and are employed in the numerical ux evaluation. We propose here to use a new reconstruction technique: the multislope MUSCL method [1315] where the reconstructed values are obtained using an approximation of specic directional derivatives instead of the full gradient. The main advantage is that the reconstruction can be rewritten as a one-dimensional MUSCL method at each interface leading to a simple and ecient scheme.

In the axisymmetric context, there are few numerical tests to validate the scheme for compressible Euler equations. For example, we are not able to compute the exact solution of the Riemann problem excepting in very particular situations. We propose a new numerical test for the swirling ow based on the steady-state situation. We manage to reduce the Euler system to an ordinary dierential equation and a simple approximation based the Euler forward algorithm is proposed to provide an accurate numerical solution.

The organization of the paper is as follows. In section 2, we present the governing Euler equations in cylindrical coordinates assuming rotational symmetry.

In section 3, we present the numerical scheme and its second-order extension using a multislope MUSCL method. In section 4, we establish a stationary solution assuming that the ow depends only on the radial direction. Finally, we present numerical experiments to test the obtained scheme.

Axisymmetric Euler equations for swirling ows

We rst present compressible Euler equations using the cylindrical coordinates and simplify them under the axisymmetric invariance assumption. For any point X = (x, y, z) ∈ R 3 we denote by (r, θ, z) the associated cylindrical coordinates. We start by giving the compressible Euler equations in the domain Ω using cartesian coordinates:

Let P = R + × R = {r ≥ 0, z ∈ R}
∂ρ ∂t + ∇ • (ρu) = 0, (1)
∂(ρu) ∂t + ∇ • (ρu ⊗ u + P I) = 0, (2)
∂E ∂t + ∇ • ((E + P )u) = 0, (3)
where ρ is the uid mass density, P is the pressure, u is the velocity vector and E is the total energy per unit volume. The tensors u ⊗ u and I stand for the tensor product of u by u and the unit tensor respectively.

To close the system, we add a state equation which in general form reads

P = P (ρ, e),
where e stands for the specic internal energy related to the total energy by:

E = ρ e + 1 2 ρ|u| 2 .
In the sequel, we restrict ourselves to an ideal gas, that is,

P = (γ -1)ρ e, ( 4
)
where γ is the ratio of specic heats at constant pressure and volume.

We introduce the mapping dening cylindrical coordinates:

(r, θ, z) → (r cos θ, r sin θ, z).

To any scalar function f : Ω → R, we associate the function f (r, θ, z) = f (r cos θ, r sin θ, z) for (r, z) ∈ Ω and θ ∈ [0, 2π). Due to the axial symmetry we look for solutions (density, velocity, pressure, . . . ) that are independent of θ, i.e. all the involved functions will depend on space variables r and z and on the time variable t.

In the sequel, for the sake of simplicity, we shall omit the tilde symbol on the functions f = f (r, z). Applying the divergence operator in cylindrical coordinates and eliminating the θderivatives, we deduce from (1)(3) the following system:

∂ρ ∂t + 1 r ∂ ∂r (rρu r ) + ∂ ∂z (ρu z ) = 0, (5)
∂ ∂t (ρu r ) + ∂ ∂r (ρu 2 r + P ) + ∂ ∂z (ρu r u z ) + 1 r ρ (u 2 r -u 2 θ ) = 0, (6)
∂ ∂t (ρu z ) + ∂ ∂r (ρu z u r ) + ∂ ∂z (ρu 2 z + P ) + 1 r ρ u z u r = 0, (7)
∂ ∂t (ρu θ ) + ∂ ∂r (ρu θ u r ) + ∂ ∂z (ρu θ u z ) + 2 r ρ u θ u r = 0, (8)
∂E ∂t + 1 r ∂ ∂r (ru r (E + P )) + ∂ ∂z (u z (E + P )) = 0. (9)
Note that assuming rotational symmetry does not reduce the number of unknowns since the velocity eld has three nonzero components that are respectively the radial u r (r, z), the tangential u θ (r, z) and the axial component u z (r, z).

To derive a conservative form of the system, we multiply equations (5)(9) by the radial coordinate r like in [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics, A practical Introduction[END_REF][START_REF] Glaister | Flux dierence splitting for the Euler equations with axial symmetry[END_REF]. We then obtain the nal form of the Euler equations:

∂ ∂t (rρ) + ∂ ∂r (rρu r ) + ∂ ∂z (rρu z ) = 0, (10)
∂ ∂t (rρu r ) + ∂ ∂r (rρu 2 r + rP ) + ∂ ∂z (rρu r u z ) = ρu 2 θ + P, (11)
∂ ∂t (rρu z ) + ∂ ∂r (rρu z u r ) + ∂ ∂z (rρu 2 z + rP ) = 0, (12)
∂ ∂t (rρu θ ) + ∂ ∂r (rρu θ u r ) + ∂ ∂z (rρu θ u z ) = -ρu θ u r , (13)
∂ ∂t (rE) + ∂ ∂r (ru r (E + P )) + ∂ ∂z (ru z (E + P )) = 0. (14)
Equivalently, we can write this system under the conservative form:

∂(rU ) ∂t + ∂(rF r (U )) ∂r + ∂(rF z (U )) ∂z = G(U ), ( 15
)
where U is the conservative variable vector dened by

U =              
3 A nite volume scheme with multislope MUSCL reconstruction 3.1 A rst-order nite volume scheme

To address numerical approximation, we describe the nite volume method we use for the conservative part of the Euler system. We consider a conforming unstructured mesh T h of the two-dimensional domain Ω, made of triangles (or cells) C i ⊂ Ω, i = 1, . . . , I. We denote by ν(i) the index set of the neighbouring triangles C j which share a common edge S ij with the cell C i and by n ij = (n ij,r , n ij,z ) the outward unit normal vector to C i .

We rst integrate the system (15) over the cell C i and use the Green formula to get

d dt C i U (r, z, t) r dr dz + ∂C i (F r (U )n ij,r + F z (U )n ij,z ) r dσ = C i G(U ) dr dz,
where ∂C i denotes the boundary of the cell C i . Let (t n = n ∆t) n∈N be a uniform subdivision of the time interval [0, +∞). Integrating on the interval [t n , t n+1 ] we obtain the equation:

C i U (r, z, t n+1 ) r dr dz = C i U (r, z, t n ) r dr dz - t n+1 t n ∂C i (F r (U ) n ij,r + F z (U ) n ij,z ) r dσ dt + t n+1 t n C i G(U ) dr dz dt.
We aim at computing a constant piecewise approximation of the solution U on each cell C i . To this end, U n i represents an approximation of the average of U on C i at time t n . More precisely, for a given cell C i , we introduce the following measures

|C i | = C i dr dz, |C i | r = C i r dr dz,
and we dene the approximation

U n i ≈ 1 |C i | r C i U (r, z, t n ) r dr dz. (16)
In the same way, we dene the measures of an edge

|S ij | = S ij dσ, |S ij | r = S ij r dσ,
and we dene the approximation of the ux across the interface S ij during the interval [t n , t n+1 ] by

F n ij ≈ 1 ∆t |S ij | r t n+1 t n S ij (F r (U ) n ij,r + F z (U ) n ij,z ) r dσ dt.
We eventually dene an approximation of the right-hand side contribution by

G n i ≈ 1 ∆t |C i | t n+1 t n C i G(U ) dr dz dt.
It results from [START_REF] Buard | Monoslope and multislope MUSCL methods for unstructured meshes[END_REF] that the integrals involve two kinds of measures: r dr dz and dr dz. To compute the mean value on the cell C i employing the r dr dz measure, we have to divide by |C i | r since the approximation ( 16) becomes an equality for constant functions. On the other hand, we divide by |C i | to provide a mean value of G n i since we use the dr dz measure in the integral (see [START_REF] Guardone | Finite element/volume solution to axisymmetric conservation laws[END_REF], p. 495).

To dene an explicit scheme, we use a numerical ux such that

F n ij = F(U n i , U n j , n ij ), G n i = G(U n i ).
The scheme reads then

|C i | r U n+1 i = |C i | r U n i -∆t j∈ν(i) |S ij | r F(U n i , U n j , n ij ) + ∆t |C i |G(U n i ). (17)
In the present study, we present two numerical uxes. We rst use the Rusanov ux

F(U i , U j , n ij ) = F r (U i ) + F r (U j ) 2 n ij,r + F z (U i ) + F z (U j ) 2 n ij,z -λ ij (U j -U i )
with λ ij = λ(U i , U j ) large enough to ensure stability (see [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics, A practical Introduction[END_REF]). The main advantage of such a ux is its ability to handle real gases. The drawbacks are on one hand a signicant diusion amount and, on the other hand, contact discontinuities are not well approximated when steady-state problems are considered (see [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics, A practical Introduction[END_REF] p. 293, [START_REF] Batten | on the choice of wavespeeds for the HLLC Riemann solver[END_REF], and [START_REF] Godlewski | Hyperbolic systems of conservations laws[END_REF] p. 204). We then resort to a second less diusive numerical ux: the HLLC ux (see [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics, A practical Introduction[END_REF] for a technical description).

For this one, we obtain a better resolution of the contact discontinuity for stationary solutions while ux computations remain reasonable for complex applications (real gas, multispecies and multiphase problems for instance).

Let us give a new expression of (17) that uses the |C i | and |S ij | measures only. For a given triangle C i , let (r i,1 , z i,1 ), (r i,2 , z i,2 ), (r i,3 , z i,3 ) stand for the coordinates of its three vertices. A similar notation is used for a given edge S ij . Using the identities

|C i | r = |C i | (r i,1 + r i,2 + r i,3 ) 3 , |S ij | r = |S ij | (r ij,1 + r ij,2 ) 2 ,
we obtain from [START_REF] Batten | on the choice of wavespeeds for the HLLC Riemann solver[END_REF],

U n+1 i = U n i -∆t j∈ν(i) |S ij | |C i | 3(r ij,1 + r ij,2 ) 2(r i,1 + r i,2 + r i,3 ) F(U n i , U n j , n ij ) + 3∆t r i,1 + r i,2 + r i,3 G(U n i ).
This new formulation enables proving that the numerical scheme preserves the situation when the gas is at rest. Indeed, consider the initial condition U 0 where the pressure P 0 and the density ρ 0 are constant and the velocity is null, the resulting numerical ux is given by

F(U 0 i , U 0 j , n ij ) =               0 P 0 n ij,r P 0 n ij,z 0 0              
.

We easily deduce that we have after the rst step

ρ 1 i = ρ 0 , u 0 θ,i = u 1 θ,i = 0 and E 1
i = E 0 since the rst and the two last components of vector G are null. It remains to study radial and axial components of the velocity u 1 r,i and u 1 z,i . From the momentum equations, the ux contribution is

C ux = j∈ν(i) |S ij | |C i | 3(r ij,1 + r ij,2 ) 2(r i,1 + r i,2 + r i,3 ) P 0 n ij = P 0 |C i | r j∈ν(i) |S ij | r n ij = P 0 |C i | r j∈ν(i) S ij rn ij dσ.
Using the Green formula we obtain

C ux = P 0 |C i | r C i ∇ r,z (r) dr dz = P 0 |C i | |C i | r    1 0    .
The right-hand side term for the radial and axial equations reduces to

C rhs = 3 r i,1 + r i,2 + r i,3 P 0    1 0    = |C i | |C i | r P 0    1 0    .
Since C ux = C rhs , then u n r,i = u n z,i = 0 and the gas stays at rest.

3.2 A second-order scheme using the multislope MUSCL method

In the early 70's, Van Leer [START_REF] Van Leer | Towards the ultimate conservative dierence scheme. V. A secondorder sequel to Godunov's method[END_REF] introduced the MUSCL technique (Monotonic Upwind Schemes for Conservation Laws) to get a more accurate approximation with less diusion eect while maintaining stability. Extensions to multidimensional situations for unstructured meshes have been proposed (see [START_REF] Hubbard | Multidimensional Slope Limiters for MUSCL-Type Finite Volume Schemes on Unstructured Grids[END_REF][START_REF] Godlewski | Hyperbolic systems of conservations laws[END_REF]). We present here a new extension of the MUSCL technique on triangles where we use approximations of the directional derivative of U as proposed in [1315] instead of an approximation of ∇U .

The goal is to compute a more accurate ux F(U n ij , U n ji , n ij ) where U n ij and U n ji are better approximations of U on both sides of the edge S ij . In the following, we detail the method to construct the U ij values where we skip the time index n for the sake of simplicity.

The fundamental decomposition

For a given volume C i , we denote by B i the centroid and Q ij the intersection of the segment [B i , B j ] with the common edge S ij for all j ∈ ν(i) (see gure. 1). We introduce the barycentric coordinates (ρ ij ) j∈ν(i) , dened by

j∈ν(i) ρ ij B j = B i , j∈ν(i) ρ ij = 1. ( 18
)
We assume that the point B i is strictly inside the triangle formed by the three other points. Therefore ρ ij > 0. Dening the normalized direction

t ij = B i B j |B i B j |
and using relation [START_REF] Godlewski | Hyperbolic systems of conservations laws[END_REF], we obtain a decomposition (referred to as the fundamental decomposition) of t ij in function of the two other directions:

t ij = k∈ν(i) k =j β ijk t ik , (19)
with the explicit expression of the β coecients

β ijk = - ρ ik ρ ij |B i B j | |B i B k | .

The slope limiter

We now focus on the reconstruction of U ij . To this end, we consider the physical variable vector

V =               ρ u r u z u θ P              
, and denote by v a generic component of the vector V , i.e. v = ρ, u r , u z , u θ or P . We construct a rst set of slopes, that will be referred to as downstream slopes, and dened by

p + ij = v j -v i |B i B j | for all C i ∈ T h , j ∈ ν(i).
The slopes represent an approximation of the directional derivative following t ij . Obviously, we have to construct ve slopes p + ij , one for each variable. We dene the upstream slopes by

p - ij = k∈ν(i) k =j β ijk p + ik for all C i ∈ T h , j ∈ ν(i).
We nally compute the slopes p ij using a limiting procedure, for example

p ij = minmod (p + ij , p - ij ).
We eventually construct the new value v ij by setting

v ij = v i + p ij |B i Q ij |. (20)
With V ij and V ji in hand, we compute the conservative variable vector U ij and U ji on both sides of the interface S ij to provide the ux evaluation

F(U ij , U ji , n ij ).
Remark 1 Thanks to the fundamental decomposition [START_REF] Vasil'evskii | Numerical simulation of equilibrium induction plasma ows in a cylindrical plasmatron channel[END_REF] we can prove that this reconstruction is exact for linear functions, i.e. v(

Q ij ) = v ij , if v is piece- wise linear.
Remark 2 The positivity of the barycentric coordinates implies that β ijk < 0.

An important consequence is that if v i is a local extremum then p + ij p - ij ≤ 0. This implies p ij = 0, which prevents the extrema from increasing (see [START_REF] Clain | The multislope MUSCL method[END_REF]). Remark 3 The main issue is that we consider only onedimensional reconstruction in each direction. This enables using any classical 1-D limiter to compute the slope p ij (see [START_REF] Clauzon | Analyse de Schémas d'ordre élevé pour les écoulements compressibles[END_REF] for more details). In the case where C i has a common edge with the boundary, we set p ij = 0 which results in a rstorder scheme.

Remark 4 It is important to perform the reconstruction with the physical variables instead of the conservative variables. Indeed, in some situations a reconstruction based on the conservative variables can yield a negative internal energy since it is obtained by subtracting the kinetic energy from the total energy.

Steady-state radial solutions

Let us investigate the particular case of steady-state radial solutions. This case is helpful in the sense that it provides a reference solution to test the chosen numerical scheme. In practice, a stationary solution is obtained by supplying adequate boundary conditions to the time dependent problem and letting the solution converge to a stationary one. In addition, as far as some compressible ows for industrial applications like in induction plasma problems (see [START_REF] Vasil'evskii | Numerical simulation of equilibrium induction plasma ows in a cylindrical plasmatron channel[END_REF] for instance) are concerned, stationary solutions are to be considered.

It turns out that it is rather dicult to compute an exact solution for stationary ows. To provide a simple case, we resort then to restrict this section to radial solutions, i.e. that do not depend on z and such that u z = 0. Under these assumptions, equations [START_REF] Cumming | Supersonic, turbulent ow computation and drag optimization for axisymmetric afterbodies[END_REF] [START_REF] Clauzon | Analyse de Schémas d'ordre élevé pour les écoulements compressibles[END_REF] reduce to:

d dr (rρu r ) = 0, (21)
d dr (r(ρu 2 r + P )) = ρu 2 θ + P, (22)
d dr (ρu θ u r ) = -ρu θ u r , (23)
d dr (ru r (E + P )) = 0. ( 24
)
Using the state equation of an ideal gas (4), the integration of equations ( 21), ( 23) and (24) gives:

ρu r = α 1 r , (25)
ρu θ u r = α 2 r 2 , ( 26
) u r γ γ -1 P + 1 2 ρ (u 2 r + u 2 θ ) = α 3 r , (27)
where α 1 , α 2 and α 3 are constants.

Thanks to relations (25) and ( 26), the expression of the tangential velocity u θ becomes

u θ = α 2 α 1 r . ( 28
)
Let us introduce two new variables: a = ρu 2 r and b = ρu 2 θ . The system becomes then:

                 ab = α 2 2 r 4 , γ γ -1 P + 1 2 (a + b) = α 1 α 3 ar 2 , d (a + P ) dr = b -a r . (29)
We now show how to reduce system (29) to an ordinary dierential equation that we solve numerically. We rst deal with the simple case where u θ = 0, i.e. b = 0. Then, we consider the more general case of a swirling ow.

4.1

First case: u θ = 0

Let us rst consider the case where the uid ow has the property u θ = 0.

Then system (29) is reduced to:

       γ γ -1 P + 1 2 a = α 1 α 3 ar 2 , d (a + P ) dr = - a r . (30)
From this we deduce

d dr γ -1 γ α 1 α 3 ar 2 + γ + 1 2γ a = - a r .
Therefore

da dr = - a r 2A -1 A -γ+1 2γ , (31)
where

A = γ -1 γ α 1 α 3 a 2 r 2 = γ -1 γ α 3 α 3 1 r 2 ρ 2 .
By dierentiating the identity α 2 1 = r 2 ρa, we deduce

dρ dr = -α 2 1 2a + ra ′ r 3 a 2 .
From (31), we get

dρ dr = α 2 1 γar 3   1 A -γ+1 2γ   .
This eventually gives the ordinary dierential equation for the density:

dρ dr = ρ α 3 α 3 1 ρ 2 r 2 - γ + 1 2(γ -1) (γ -1)r . ( 32
)
The numerical solution of equation ( 32) is described in section 4.3.

Remark 5 An analytical solution of equation ( 32) can be obtained if we assume that

α 3 α 3 1 ρ 2 r 2 ≫ γ + 1 2(γ -1)
.

This assumption is valid in particular for subsonic ows at ambient temperature for instance. Neglecting the higher order term, the dierential equation (32) reduces to

2ρ dρ dr = α 3 1 (γ -1)α 3 2 r 3 .
Using the initial condition ρ(r = R 0 ) = ρ 0 , we obtain the solution

ρ(r) = ρ 2 0 + α 3 1 α 3 (γ -1) 1 R 2 0 - 1 r 2 .

The general case (swirling ow)

We now deal with the general case where u θ = 0. Dierentiating the rst equation of system (29) yields

b da dr + a db dr = - 4ab r . (33)
Thus

db dr = - b a da dr - 4b r . ( 34
)
Let us substitute P in the system (29) and use ( 33) and (34). We get

d dr γ -1 γ α 1 α 3 ar 2 - γ -1 2γ b + γ + 1 2γ a = b -a r .
We compute the derivative of the previous expression and eliminate the derivative of b by using (34),

γ -1 γ a b 2 - α 1 α 3 ar 2 da dr = a r b a - γ -1 γ a 4b 2 - 2α 1 α 3 ar 2 -1 .
Since α 2 1 = r 2 ρa, we obtain

dρ dr = -C 2 1     2ra + r 2 da dr r 4 a 2     .
After calculation, the two previous equations give:

dρ dr = - α 2 1 ar 3       a + b (γ -1) b 2 - α 1 α 3 ar 2 + (γ + 1) 2 a      
.

Using the denitions of ρ and P , we nally obtain the ordinary dierential equation for the uid density:

dρ dr = ρ r a + b γP -a . (35)
Remark 6 In equation ( 35), functions a, b, P depend only on ρ, r and the constants. Indeed, for a given ρ, we can compute b = ρu 2 θ using (28). We then deduce a with the rst relation of system (29). We nally compute P with the help of the second relation of system (29).

A numerical method for stationary radial solutions

We consider a uniform subdivision of the domain [R 0 , R 1 ] with K elements where r k = R 0 + k∆r, k = 0, ..., K with ∆r = (R 1 -R 0 )/K. To avoid the singularity at r = 0, we have taken R 0 > 0. To obtain an approximation of the density ρ for the equations (32) and (35), we use the explicit Euler method:

case u θ = 0 : ρ k+1 = ρ k + ∆r ρ k α 3 α 3 1 ρ 2 k r 2 k - γ + 1 2(γ -1) (γ -1)r k , (36)
case u θ = 0 :

ρ k+1 = ρ k + ∆r ρ k r k a k + b k γP k -a k . ( 37
)
The other parameters, namely the radial velocity u r , the tangential velocity u θ and the pressure P , are deduced using:

u r,k+1 = α 1 ρ k+1 r k+1 , u θ,k+1 = α 2 α 1 r k+1 , P k+1 = γ -1 γ α 1 α 3 ρ k+1 u 2 r,k+1 r 2 k+1 - 1 2 ρ k+1 (u 2 r,k+1 + u 2 θ,k+1 ) .
The a and b values are then updated by relations:

a k+1 = ρ k+1 u 2 r,k+1 , b k+1 = ρ k+1 u 2 θ,k+1 ,

Numerical tests

We present several numerical tests to show convergence and accuracy properties of the MUSCL scheme presented throughout this paper. The rst test shows that the algorithm preserves the situation when the gas is at rest which proves that the scheme is well-balanced in this sense.

A second series of tests concerns radial steady-state solutions. We observe convergence for large time to the stationary solution computed with the numerical method proposed in subsection 4.3.

A third series of numerical experiments is dedicated to the Riemann problem.

We consider a domain Ω = [0, 1] × [0, 1] and impose an initial condition made of two constant states for z < 1/2 and z > 1/2. We obtain a classical onedimensional Sod tube problem in the Oz direction. We then compare the numerical solutions to the analytical ones.

Finally, we describe two more signicant test problems. The rst one is an adaption of the well-known front-facing step problem of [START_REF] Woodward | The Numerical Simulation of Two-Dimensional Fluid Flow with Strong Shocks[END_REF] for an axisymmetric conguration where a stationary Mach 3 ow hits a step. The second test produces a converging spherical shock simulating a gas bubble compression. This test problem allows to check the capacity of the numerical scheme to preserve the spherical symmetry of the problem.

Computations have been carried out using the nite element / nite volume library OFELI described in [START_REF] Touzani | An Object Finite Element Library, Copyright c 19982008 Rachid Touzani[END_REF]. We detail all these numerical experiments in the sequel.

Conservation of the gas at rest

In this test, we check the capacity of the algorithm to preserve a situation where the gas is at rest. We consider a domain Ω = [0, 1] × [0, 1] and prescribe the following initial conditions:

u r = u θ = u z = 0, P = 1 × 10 5 Pa, ρ = 1.2 kg•m -3
We also consider a second test where the initial null axial velocity is given by u z = 100 m•s -1 .

On the boundaries r = 0 and r = 1 we impose a reection condition using the ghost-cell technique [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics, A practical Introduction[END_REF] while we enforce a Dirichlet condition for the boundaries z = 0 and z = 1 equal to the initial condition. In both situations, we observe the preservation of the gas at rest. A few numerical artifacts of order 10 -14 which correspond to the precision of the computer are present.

Stationary solution

The goal of the study is to consider steady-state radial solution approximations and show that any non-stationary solution of the system (15) using judicious boundary conditions U inlet and U outlet converges asymptotically to the stationary solution given in section 4.

To compute the stationary solution, we use a subdivision following the onedimensional radial direction (R 0 = 0.1, R 1 = 0.4) of 10000 cells. Since the problem reduces to an ordinary dierential equation, we prescribe the initial condition at point R 0 (named the inlet condition U inlet ). Using the algorithm given in subsection 4.3, we obtain a numerical approximation on domain [R 0 , R 1 ] and we denote by U outlet the solution at point R 1 .

To compute the non-stationary solution, we use the rectangular domain Ω = [0.1, 0.4] × [0, 0.1] discretized in 6800 triangle cells (approximately 100 cells along the radial direction). We prescribe the reection condition at the boundaries z = 0 and z = 0.1 while we impose the inlet and outlet conditions at the boundaries r = R 0 and r = R 1 respectively. We initialize the solution with the inlet condition. Two numerical experiments have been performed, the rst one concerns the simple case where u θ = 0. The second test deals with the general case which corresponds to the swirling ow.

5.2.1

Case u θ = 0

We use the numerical scheme (36) and impose the following inlet conditions at r = R 0 :

C inlet =               ρ inlet = 2 u r,inlet = 200 u θ,inlet = 0 u z,inlet = 0 P inlet = 2 × 10 5              
.

We compute all the variables of vector U at each point of the subdivision, in particular we obtain the outlet condition at r = R 1

C outlet =               ρ outlet = 2.283222 u r,outlet = 43.79776 u θ,outlet = 0 u z,outlet = 0 P outlet = 2.407435 × 10 5              
.

We now introduce the two boundary conditions in the non-stationary scheme and run until a stationary situation is achieved.

In gure 3, we show the density, the pressure and the radial velocity in function of r at z = 0.05. We also propose a zoom of the density distribution to highlight the dierence between the Rusanov and the HLLC solver but also between the rst and the second order scheme. These two solutions clearly coincide and as we can expect, the second order scheme with the HLLC ux provides the best approximation.

5.2.2

Case u θ = 0

We now deal with the swirling ow where we use the algorithm (37) and the new inlet conditions:

C inlet =               ρ inlet = 2 tu r,inlet = 200 u θ,inlet = 10 u z,inlet = 0 P inlet = 2 × 10 5               .
We also provide the outlet condition:

C outlet =               ρ outlet = 2.283956 u r,outlet = 43.78367 u θ,outlet = 2.5 u z,outlet = 0 P outlet = 2.408519 × 10 5               .
With the two boundary conditions in hand, we compute the non-stationary solution until a stationary solution is obtained. Figure 4 shows density and velocity components u r and u θ . We also include a zoom of the density to highlight the second order scheme eciency. As in the previous case, we succeed in reaching the steady-state solution and the second order scheme with the HLLC ux provides the best solution.

A shock tube test

We decompose the domain

Ω = [0, 1] × [0, 1] into two subdomains Ω L = [0, 1] × [0, 1 2 ] and Ω R = [0, 1] × [ 1 2
, 1] and denote by D the line D = {(r, 1

2 ); r ∈ [0, 1]}. We then consider the Sod tube problem prescribing the initial condition

U (t = 0) = U L in Ω L and U (t = 0) = U R in Ω R .
We also impose the reection condition at the boundary of Ω. It results that the solution is invariant with respect to r and we obtain a classical one-dimensional Sod tube problem in the Oz direction with an initial discontinuity at z = 1/2.

We compare exact solutions of the Riemann problem to the numerical approximations obtained using the rst and second-order schemes and the Rusanov and HLLC schemes. All the computations have been performed with an unstructured mesh of 20748 elements.

First conguration: Rarefaction and shock

We consider a rst conguration providing a rarefaction for the u -c characteristic eld and a shock for the u+c characteristic eld. We denote by U * L , the state just after the rarefaction while U * R represents the state just before the shock. Of course, a contact discontinuity occurs between the two intermediate states. We sum up in table 1 the used values in the test. We solve the Riemann problem with the algorithm for axisymmetric congu- Internal energy (J.kg -1 ) Fig. 5. Comparison of density and internal energy obtained by rst-order Rusanov scheme (dashed) and rst-order HLLC scheme (dotted) with the exact solution (continuous) at time t = 0.8 ms for the rst conguration. energy between the exact solution and the approximations using the Rusanov and HLLC ux with the rst order scheme. In gure 6, a similar comparison is shown but we use the second order MUSCL technique. The solutions are ) Fig. 6. Comparison of density and internal energy obtained by second-order Rusanov scheme (dashed) and second-order HLLC scheme (dotted) with the exact solution (continuous) at time t = 0.8 ms for the rst conguration. in good agreement with the exact solution and the second order technique provides the best approximations.

Second conguration: double shock

Here; we are concerned with the double shock conguration. Table 2 lists the states we employ in this test. We compute the solution approximations with the algorithm dedicated to axisymmetric geometries and visualize the density and the internal energy using the rst order scheme (gure 7) or the second order scheme (gure 8).

As expected, the MUSCL technique reduces the diusion eect close to the discontinuities. We also remark that the Rusanov and HLLC schemes provide the same approximation quality in the second order context whereas the HLLC ux is less diusive if one uses a rst-order scheme. Internal energy (J.kg -1 ) Fig. 8. Comparison of density and internal energy obtained by a second-order Rusanov scheme (dashed) and a second-order HLLC scheme (dotted) with the exact solution (continuous) at time t = 1.984 ms for the second conguration.

Third conguration: double rarefaction

The third conguration is composed of two symmetric rarefaction waves where the solution involves a state near vacuum. We list in table 3 the states we use for the calculation. Approximations have been computed using the rst (gure 9) and the second order technique (gure 10). The main issue is the internal energy approximation. First order methods give a non physical peak at z = 0.5 while second order methods reduce this numerical artifact. We note that the HLLC ux with the MUSCL technique reduces signicantly the undesired peak. Internal energy (J.kg -1 ) Fig. 9. Comparison of density and internal energy obtained by rst-order Rusanov scheme (dashed) and rst-order HLLC scheme (dotted) with the exact solution (continuous) at time t = 0.6 ms for the third conguration. ) Fig. 10. Comparison of density and internal energy obtained by second-order Rusanov scheme (dashed) and second-order HLLC scheme (dotted) with the exact solution (continuous) at time t = 0.6 ms for the third conguration. This numerical test deals with a Mach 3 ow in a wind tunnel with a forwardfacing step. This test is analogous to the well-known test problem of [START_REF] Woodward | The Numerical Simulation of Two-Dimensional Fluid Flow with Strong Shocks[END_REF]. We have considered here a discretization of the domain into 21186 nite volume cells. A stable shock wave pattern develops after a few time units. We impose the inow boundary conditions with ρ = 1.4, u z = 3, u r = u θ = 0 and P = 1. Figure 11 shows the snapshot of ow eld at t = 15 s.

(kg • m -3 ) u z (m • s -1 ) P (

  denote an axial cut of the three dimensional space (the set of parameters) and let Ω be an open set of P. The open set Ω ⊂ R 3 will denote the threedimensional volume obtained by rotation around the axial direction 0z, i.e. Ω := {(r cos θ, r sin θ, z); (r, z) ∈ Ω, 0 ≤ θ < 2π}.

Fig. 1 .

 1 Fig. 1. Geometrical ingredients and notations.

Fig. 2 .

 2 Fig.2. Radial, axial and tangential velocity distribution following the radial direction at z = 0.5 for the rst test (Left) and radial velocity following the radial direction at z = 0.5 for the second test (Right).

Fig. 3 .

 3 Fig. 3. Comparison of density, pressure and the radial velocity distributions versus the radial direction between the stationary solution and numerical results.

Fig. 4 .

 4 Fig. 4. Comparison of density, radial and tangential velocity distributions versus the radial direction between the stationary solution and numerical results.

ρ

  (kg • m -3 ) u z (m • s -1 ) P (Pa) 1.08767 × 10 4 rations using the Rusanov and HLLC ux and using the rst and second order techniques. We present in gure 5 a comparison of the density and internal

ρ

  (kg • m -3 ) u z (m • s -1 ) P (Pa)

Fig. 7 .

 7 Fig.7. Comparison of density and internal energy obtained by rst-order Rusanov scheme (dashed) and a rst-order HLLC scheme (dotted) with the exact solution (continuous) at time t = 1.984 ms for the second conguration.

5. 4 A

 4 Mach 3 wind tunnel with a forward-facing step

Fig. 11 .

 11 Fig.11. Forward-facing step using a second-order HLLC scheme. 100 contours of the density (left) from 0.762 to 5, pressure (centre) from 0.462 to 11.8 and Mach number (right) from 0.077 to 3.02.5.5Converging spherical shock test

Fig. 12 .

 12 Fig. 12. Density at t = 0.2, 31 contours from 1 to 4. Rusanov rst-order (top left) and second-order (top right), HLLC rst-order (bottom left) and second-order (bottom right).

Fig. 13 .

 13 Fig.[START_REF] Clain | The multislope MUSCL method[END_REF]. Comparison of density at t = 0.2 for r = 0 and z = 0 using rst-order (Left) and second-order (Right) schemes.

Table 1

 1 Data for the rst Riemann problem

Table 2

 2 Data for the second Riemann problem

Table 3

 3 Data for the third Riemann problem

	ρ

Conclusion

We have proposed a new second order cell-centered nite volume formulation for compressible Euler equations using cylindrical coordinates to compute axisymmetric solutions. The nite volume scheme is based on a mean value approximation using the r dr dz metric instead of the traditional dr dz metric to simplify the boundary condition on the Or axis. An original MUSCL method (multislope MUSCL method) is introduced where an estimate of the directional derivatives is performed to construct more accurate approximations of the solution on both side of the edges. New numerical tests based on the steady-state situation has been proposed to validate the numerical scheme and to compare the dierent schemes performances.

and Oz axis.

To simulate the spherical bubble compression using the axisymmetric model, we consider a quarter of disk D = {(r, z); √ r 2 + z 2 ≤ 1/2} in the unit square Ω = [0, 1] × [0, 1] and we state the following initial conditions:

(1, 0, 0, 0, 1) in D, (4, 0, 0, 0, 4) in Ω \ D.

We also prescribe reecting boundary conditions on the whole boundary. We have performed numerical tests with a discretization of the domain into 20640 nite volume cells.

In gure 12 we reproduce the density distribution at time t = 0.2 using the Rusanov and the HLLC uxes. A small asymmetry appears when the Rusanov ux is employed but it disappears when we use the MUSCL procedure. This is conrmed by the visualization of the density value along the Or and Oz axis in gure 13.


C Opitom

A Guilbert-Lepoutre

S Besse

B Yang

C Snodgrass

MUSE observations of comet 67P/Churyumov-Gerasimenko A reference for future comet observations with MUSE

Keywords: Comets, individual, 67P/Churyumov-Gerasimenko,Techniques, imaging spectroscopy

Aims. Observations of comet 67P/Churyumov-Gerasimenko were performed with MUSE at large heliocentric distances postperihelion, between March 3 and 7, 2016. Those observations were part of a simultaneous ground-based campaign aimed at providing large-scale information about comet 67P that complement the ESA/Rosetta mission. Methods. We obtained a total of 38 datacubes over 5 nights. We take advantage of the integral field unit (IFU) nature of the instrument to study simultaneously the spectrum of 67P's dust and its spatial distribution in the coma. We also look for evidence of gas emission in the coma. Results. We produce a high quality spectrum of the dust coma over the optical range that could be used as a reference for future comet observations with the instrument. The slope of the dust reflectivity is of 10%/100 nm over the 480-900 nm interval, with a shallower slope towards redder wavelengths. We use the Afρ to quantify the dust production and measure values of 65±4 cm, 75±4 cm, and 82±4 cm in the V, R, and I bands respectively. We detect several jets in the coma, as well as the dust trail. Finally, using a novel method combining spectral and spatial information, we detect the forbidden oxygen emission line at 630 nm. Using this line we derive a water production rate of 1.5 ± 0.6 × 10 26 molec./s, assuming all oxygen atoms come from the photo-dissociation of water.

Introduction

The ESA/Rosetta mission was the first mission to approach and follow a comet for a large part of its orbit, as well as to land on its nucleus. The Rosetta orbiter and the Philae lander were equipped with a series of instruments to study in-situ the gas and dust coma surrounding comet 67P/Churyumov-Gerasimenko (hereafter 67P) and the properties of its nucleus. The mission provided unprecedented insights into the structure, composition, and activity evolution of 67P, and comets in general (Taylor et al. 2017).

In order to complement the mission, an ambitious groundbased campaign was put together, with dozens of telescopes across the globe and in space performing observations of the comet at the same time as the Rosetta mission (Snodgrass et al. 2017b). Those observation probed the large scale coma of 67P and provided context to the measurements performed by the Rosetta mission. More importantly, they provide a link with ground-based observations of a large number of comets that have been performed over the last century. As part of this campaign, observations were performed with the Multi Unit Spectroscopic Explorer (MUSE) instrument at the Very Large Telescopes (VLT). The observations were performed in March 2016, after the comet's perihelion passage, when the comet was moving away from the Sun. In this paper, we present the results of the MUSE observations of comet 67P.

Observations and data reduction

MUSE is an integral field unit spectrograph mounted on the UT4 telescope of the VLT, in Chile [START_REF] Bacon | Proc. SPIE[END_REF], nominally covering the 480-930 nm range. In Wide Field mode, which was used for our observations, MUSE has a field of view (FoV) of 1 ×1 covered without gaps. The instrument has a platescale of 0.2 /pix and a spectral resolving power between 1770 (at 480 nm) and 3590 (at 930 nm). Observations of comet 67P were performed between 2016 March 3 and March 7. In total, 38 datacubes were obtained over 5 nights. For 4 of those cubes, the comet was not visible, or very badly centred, so that we discarded the data. The sky during those 5 nights was either clear or photometric, but the seeing was variable, between 0.6 and 2.3 . For all observations, we used an exposure time of 600s, and the position angle of the instrument was set to 0 • (we did not apply any rotation between the exposures). At the time of the observations, even though the comet was still active, and thus extended, it did not fill the entire MUSE FoV. Because of that, no dedicated sky observations were performed. The observing circumstances are presented in Table 1.

The data reduction was performed using the ESO pipeline [START_REF] Weilbacher | MUSE-DRP: MUSE Data Reduction Pipeline[END_REF], with the sky estimated from regions near the edge of the cubes that are free from the comet contamination. In addition to the cube reconstruction, bias subtraction, flatfield correction, wavelength calibration, and sky subtraction, the ESO MUSE pipeline also corrects for the telluric absorption and flux-calibrates the science data, using a standard star observed the same night as the science observations. Even though the sky was estimated directly on the science cubes and subtracted by the pipeline, while examining the reduced cubes we noticed relatively strong sky residuals. In order to reduce those, we then used the ZAP (Zurich Atmosphere Purge) software (Soto et al. 2016). ZAP is a Principal Component Analysis-based software designed to perform sky subtraction on IFU data. The use of ZAP allows us to remove most of the sky residuals left after the pipeline reduction. As will be discussed later, we also performed a full data reduction without applying any sky subtraction, in order to search for forbidden oxygen emission lines in the coma of 67P.

As mentioned before, the ESO pipeline performs a correction of the telluric absorption using a standard star observed on the same night as the science data. However, the standard stars are primarily used for flux calibration and are not optimal for telluric correction. Also, they are not always observed just after the science observations nor at the same airmass. Because of that, in the reduced cubes, there are residuals of the strong O 2 telluric band around 760 nm if the telluric correction is performed with the pipeline. For all the analysis focused on the 2D structure of the dust coma, and the detection of forbidden oxygen lines, this does not impact the quality of our measurements, as those residuals are relatively constant over the field and we extracted the cubes over wavelength ranges that mostly avoid that region of the spectrum. However, in Section 3.1, we focus on the spectrum of the dust in the coma of 67P and residual from the telluric correction could impact the quality of the spectrum presented. For the spectra presented in that section, we thus used the Molecfit software (Smette et al. 2015;Kausch et al. 2015) to perform a better telluric correction. Molecfit has proven to be very useful to provide an accurate telluric correction even when dedicated observations of a standard star were not available. We applied Molecfit directly on the extracted 1D spectra. The fit of the telluric features was performed for each extracted spectrum individually.

Analysis

A reference dust spectrum

We first study the spectrum of comet 67P over the optical range. To do so, we extract all the data cubes over a 5 pixels radius (1 ) aperture around the comet optocenter. We chose such a small aperture to focus on the part of the coma where the signal is the strongest, to avoid having strong sky residuals in the presented spectrum. After the extraction, the spectra are corrected for the telluric features using the Molecfit software, as explained above. All 34 spectra are then median-combined to produce a high quality spectrum of the dust in the coma of 67P. This spectrum is shown in the upper part of Fig. 1. We do not see any emission band in the spectrum, because of the large heliocentric distance and low activity level of the comet at the time of the observations. In the same part of the Figure, we show a reference solar spectrum obtained using the SOLar SPECtrometer (SOL-SPEC) instrument of the SOLAR payload on board the International Space Station (ISS) (Meftah et al. 2018), re-sampled to match the sampling of the 67P spectrum. Given that the aperture we chose is very small, the spectrum shown in Fig. 1 probably contains a non-negligible nucleus contribution. For comparison purposes, we also extract the dust spectrum over a 10 aperture, in which the dust coma dominates (blue dashed curve in the top part of Fig. 1). The two spectra match well, except for a small difference around 880 nm. This region is severely affected by telluric absorption and has strong sky emissions, that might explain this difference. In general, we can thus say that the spectrum extracted over the 1 aperture is representative of the coma dust and we will discuss this spectrum only in the following text.

In the bottom part of Fig. 1, we divide the comet spectrum by the Solar reference spectrum to compute the relative reflectance of the dust in the coma of 67P. The reflectance spectrum is normalised at 600 nm. Overlaid to the reflectance spectrum measured with MUSE, is the reflectance measured using the X-Shooter spectrograph in November 2014 over the same wavelength range (Snodgrass et al. 2016). The reflectance spectra measured in November 2014 and March 2016 are consistent with each other, indicating that the dust reflectance as measured from the ground is similar at large heliocentric distance pre-and post-perihelion. We do not see any sign of absorption bands in the optical spectrum of 67P. We measure a spectral slope for the reflectance, or dust reddening, of 10%/100 nm in the 480-900 nm interval. We also notice that the slope becomes shallower at longer wavelengths. We measure a slope of 13 %/100 nm in the 500-700 nm interval but only of 5%/100 nm in the 700-900 nm interval. This is fully consistent with what is reported from X-Shooter observations performed in 2014, from which values between 10%/100 nm and 20%/100 nm were reported in the 550 to 1000 nm interval, the shallower values corresponding to the red end of the wavelength range.

Our measurements are also consistent with in-situ measurements of the dust from the ESA/Rosetta mission. Bertini et al. (2017) report a reddening in the interval 376-744 nm ranging between 11 and 14%/100 nm from measurements with the OSIRIS cameras. Similarly, La [START_REF] Forgia | EPSC-DPS Joint Meeting 2019[END_REF] report slopes measured in the inner coma ranging from 12 to 16%/100 nm between 480 and 649 nm. Those slopes are similar to those measured for the nucleus of 67P, both in-situ and from groundbased observations. Fornasier et al. (2015) report average slopes of 11-16%/100 nm over 250-1000 nm with shallower slopes toward longer wavelengths, from observation of 67P's nucleus with the Rosetta/OSIRIS camera. From ground-based observations, Tubiana et al. (2011) measure a slope for the nucleus of 12 ± 1%/100 nm over 430-850 nm, slightly shallower in the 500-850 nm range. Finally, the dust reddening measured in the coma of 67P is consistent with what is measured usually in the coma of active comets, typically between 0 and 20%/100 nm, and with a shallower slope towards the near-IR (Solontoi et al. 2012;Jewitt & Meech 1986).

Dust coma morphology and activity

For each cube in our dataset, we extract maps over the bandpasses of the V, R, and I Johnson-Cousins filters. All maps in the same bandpass are re-centred and then co-added. The resulting co-added maps are displayed in the top part of Fig. 2. The dust coma morphology is the same for all three bandpasses. It is asymmetrical, probably due to the presence of dust jets. To investigate further the presence of jets, we divide the dust maps by an azimuthal median profile [START_REF] Samarasinha | Cometary Coma Image Enhancement Facility[END_REF]). Enhanced maps are displayed on the bottom part of Fig. 2. On those we can clearly see two jets, which are most likely the cause of the apparent asymmetry of the dust coma. The first jet is located close to the anti-sunward direction, towards the South-West. The second jet is located about 90 • away, towards the South-East. Finally, we see a faint feature towards the sunward direction. Those jets are consistent with what was observed on previous passages and what is reported by Knight et al. (2017) and Snodgrass et al. (2017b) from observations at the same epoch. They are also consistent with the modelling of the pole orientation and active region location done by Vincent et al. (2013) prior to the Rosetta mission. In addition to the two jets mentioned above, we see an enhancement towards the North-West. This corresponds to the dust trail, that was reported to be at least two degrees long at that epoch (Snodgrass et al. 2017b;Boehnhardt et al. 2016;Knight et al. 2017). We do not see changes of the coma morphology over the 5 nights during which we have observations, nor over a single night.

To constrain the comet activity at that time, we compute the Afρ value, which is a proxy for dust production, as defined by A' Hearn et al. (1984). We compute the Afρ for the V, R, and I bands over a 2500 km physical aperture. We use a 2500 km aperture instead of the more commonly used 10000 km aperture because the comet signal drops significantly at 10000 km. We obtain values of 65±4 cm, 75±4 cm, and 82±4 cm in the V, R, and I bands respectively. Those values have not been corrected for the phase angle effect and are consistent with those reported by Boehnhardt et al. (2016) at the same epoch, and comparable (even though slightly lower) to those reported by Knight et al. (2017) once the phase angle effect is taken into account. Using the Afρ values in the different bands, we can also compute the dust reflectivity gradient defined as (Jewitt & Meech 1986):

S = Afρ 2 -Afρ 1 Afρ 2 + Afρ 1 20000 λ 2 -λ 1
where Afρ x and the corresponding wavelength λ x are expressed in cm and nm, respectively. For the I-V combination, this gives a reflectivity gradient of 7±1 %/100 nm, while it is of 11±2 %/100 nm for the I-R and 7±1 %/100 nm for the R-V combination. This confirms the trend outlined in section 3.1 of higher reflectivity gradient at lower wavelengths. To check for trends with the aperture size, we computed the Afρ in apertures of 5000 km, 7500 km, and 10000 km. Within the error bars, we do not see significant changes in the Afρ or reflectivity gradi- ent values with the aperture size in our data, as expected for a steady-state coma.

Gas detection

Among the features usually observed in the optical spectrum of comets, the oxygen forbidden line at 630 nm is one of the brightest. There are three forbidden oxygen lines that can be detected in the coma of comets at optical wavelengths: the green line at 557.7 nm, and the red doublet at 630 and 636.4 nm. Those lines are emitted by the decay of oxygen atoms in a metastable ( 1 D) or ( 1 S) state. Exited atomic oxygen is mainly produced through the photo-dissociation of H 2 O, CO, CO 2 , or even O 2 (see, e.g. Cessateur et al. (2016)). Observing forbidden oxygen lines then represents an opportunity to constrain the production of those potential parent species from optical observations. However, the cometary forbidden oxygen lines are often blended with the equivalent atmospheric lines, unless observed at high spectral resolution.

The spectral resolution of MUSE does not allow us to resolve the telluric and cometary lines for the geocentric velocity of comet 67P at the time of our observations. In theory, the sky subtraction should subtract the atmospheric contribution, allowing us to recover the cometary signal. However, the noise introduced by the subtraction of the very strong atmospheric features prevents us from detecting any cometary forbidden oxygen emission lines. At the time of our observations the comet was at 2.5 au post-perihelion and was only weakly active, as confirmed by the fact that we do not detect any emission lines in the spectrum presented in Fig. 1. If the oxygen lines are present, they are thus very faint and masked by the noise introduced by the sky subtraction.

We therefore attempt to detect the forbidden oxygen lines using another method. We reduce the full dataset without performing any sky subtraction. For each spaxel of each cube, we then subtract the continuum sky contribution and the dust contribution underlying the oxygen forbidden lines by defining continuum region on both sides of the red doublet and the green oxygen lines, fitting a line through those regions and subtracting it. Finally, we extract the datacubes over a very narrow wavelength range centred on the wavelength of the three oxygen lines. In low/medium resolution spectra, forbidden oxygen lines can be contaminated by NH 2 or C 2 lines, but this is highly unlikely in the case of those observations since the comet is weakly active and no emission lines are detected in Fig. 1. The result is maps of the flux contained in the forbidden lines at 557.7, 630, and 636.4 nm over the whole MUSE FoV. Since the atmospheric and cometary lines cannot be resolved by MUSE, the maps contain the sum of the atmospheric and potential cometary contribution. All maps are re-centred so that the optocenter of the comet is placed at the same position. We average all 34 maps for each line, performing a 3-σ clipping. The result in shown in Fig. 3 In those maps, we expect the atmospheric contribution to be relatively uniform over the 1 ×1 FoV of MUSE. The comet contribution, on the other side, would be concentrated around the optocenter of the comet, with a limited spatial extension since the states responsible for the emission lines are metastable states. In Fig 3, we can see that for the 557.7 and 636.4 nm lines, the maps are mostly uniform over the whole FoV. Small variations are observed but they most likely come from non-optimal corrections of the detector-to-detector effects. We thus conclude that the signal in those maps comes from the atmospheric oxygen lines. In the 630 nm map, however, in addition to inhomogeneities similar to those seen in the other two maps, we have a clear over-density located around the position of the comet optocenter. This indi-cates that, in addition to the atmospheric line, we detect signal from the cometary forbidden oxygen line at 630 nm. It is not surprising that we only detect the 630 nm line, since it is the strongest of the three oxygen lines. The signal we measure in the comet aperture for the 630 nm map is barely 3-σ above the sky background variation (measured from four adjacent sky apertures, see later). Since the ratio between the 630 nm and 636.4 nm maps is expected to be 3 (as it is the case for the background in the two maps, see Fig. 3), it is then consistent with the fact that we do not detect significant cometary signal in the 636.4 nm map. We note that the [OI] cometary signal is enhanced towards the South, similarly to what is observed for the dust. However, we do not believe the signal is due to residuals from the dust continuum. Indeed, the technique used for the dust subtraction is the same of all 3 lines and only the 630 nm map (with the strongest oxygen feature) shows an over-density around the optocenter. Control maps built by using a wavelength range adjacent to the 630 nm line do not show a similar signal.

The cometary signal we detect in the 630 nm map is faint, but we use this signal to estimate the water production rate of the comet at that time. The 630 nm line comes from the decay of atomic oxygen in ( 1 D) state. As mentioned before, oxygen in that state is mainly produced by the photo-dissociation of H 2 O, CO, CO 2 , or O 2 . For comets at 1 au from the Sun, H 2 O photodissociation is the dominant source for producing metastable oxygen, but at distance above 2.5 au other molecules such as CO and CO 2 start to contribute (Decock et al. 2013;[START_REF] Mckay | [END_REF]. Long-term measurements of the production rate of all 4 species of interest derived from the ROSINA instrument onboard Rosetta are presented by [START_REF] Laeuter | [END_REF] and Combi et al. (2020). At the time of our observations, the CO and O 2 production rates are more than a factor 10 and 100 lower than that of H 2 O, respectively. We can thus assume that those species will contribute little to the production of metastable oxygen. The CO 2 production rate is only about a factor 2 to 4 lower than the H 2 O production rate. The emission rate for O( 1 D) production from CO 2 is 1.5 times higher than that of H 2 O (Bhardwaj & Raghuram 2012). Taking this into account, a significant part (up to half) of O( 1 D) atoms could be produced by the photo-dissociation of CO 2 . For the measurements presented below, we will assume that water is the main source for the production of metastable oxygen, so that the water production rate we derive might be overestimated.

In order to derive water production rates, we follow the procedure described in [START_REF] Schultz | [END_REF]; Morgenthaler et al. (2001). We consider a photo-chemical model including the following three reactions,

H 2 O + hν -→ H 2 + O( 1 D) (1)
H 2 O + hν -→ H + OH (2) OH + hν -→ H + O( 1 D) (3)
the water production rate is given by:

Q(H 2 O) = Q([OI]) BR1 + (BR2)(BR3)
where BR1 and BR2 are the branching ratios for reactions 1 and 2 (Huebner et al. 1992) (equal to 0.05 and 0.855 for the quiet Sun, respectively) and BR3 is the branching for reaction 3 from Morgenthaler et al. (2001) (equal to 0.094). The production rate of oxygen atoms in ( 1 D) state, Q([OI]), can be obtained from:

Q([OI]) = 4 3 4π∆ 2 I 630 AC
where ∆ is the comet geocentric distance (in cm), I 630 is the intensity of the [OI] emission (in photons/s/cm 2 ), and AC is the aperture correction factor correcting for the [OI] emission not encompassed in the aperture. We measure the comet flux in the 630 nm map using a 25 pixels (5 ) radius circular aperture centred on the comet optocentre. In order to measure and subtract the sky flux, we measure the flux in four 5 apertures located at four different positions with respect to the comet optocenter ((+50 pix,+50 pix), (+50 pix,-50 pix), (-50 pix,+50 pix), (-50 pix,-50 pix)), computed the average, and subtracted it from the flux measured in the comet aperture. Given the faintness of the cometary emission, the main uncertainty in the determination of the water production rate comes from the measurement of the sky. To estimate this uncertainty, we used the upper and lower values measured in the individual sky apertures. The resulting intensity of the comet 630 nm line is 1.8 ± 1.0 × 10 -3 photons/cm 2 /s. The aperture we used encompasses all the visible emission from the comet, so that we set the aperture correction factor to 1. This results in a water production rate of 1.5 ± 0.7 × 10 26 molec./s. This represents the actual water production rate only if all the oxygen atoms in ( 1 D) state are produced by the photo-dissociation of water.

No other measurement of the comet water production rate was reported from ground-based observations at a similar epoch because the comet was too faint and weakly active. In-situ measurements made with the ROSINA DFMS instruments determined a water production rate around 2.3 × 10 26 molec./s at the same heliocentric distance (Hansen et al. 2016). Similar values are reported by Combi et al. (2020); [START_REF] Laeuter | [END_REF] using the same instrument. Biver et al. (2019) report a water production rate of 8.5 ± 2.5 × 10 25 molec./s from measurements with the MIRO instrument. Given that CO 2 could potentially contribute up to half of the 630 nm [OI] line brightness, the actual water production rate we measure could be as low as 0.7 × 10 26 molec./s, which is in very good agreement with the Biver et al. (2019) measurement. Given the difference in technique, scale of the observations and models used, our measurement is close to those reported using Rosetta instruments in particular the MIRO instrument.

Discussion and Conclusions

Comet 67P was observed with the MUSE IFU over 5 nights, between March 3 and 7, 2016 when the comet was at 2.5 au from the Sun and 1.5 au from the Earth. The spectrum of the dust coma of 67P presented here is of high quality. It has a good signal to noise ratio given the faintness and the distance of the comet at the time of our observations. Due to the use of the Molecfit software it is little affected by telluric features at near-IR wavelengths, which is not the case of most comet spectra in the 800-900 nm range. It matches very well other ground-based observations of 67P, as well as in-situ measurements from the Rosetta spacecraft. It is also representative of what is usually observed for active comets in general. Finally, it does not contain gas emission features above the noise level (some [OI] signal is present, merged with the sky line; see below). For all those reasons, this spectrum could be used in the future as a 'template' dust spectrum to help perform the subtraction of the dust-reflected continuum for observations of other comets with the MUSE IFU. This spectrum could replace the observation of solar analogs since so far very few good solar analogs have been observed with MUSE. The dust reddening can vary from one comet to another, so that this spectrum would need to be corrected from any slope difference between the target and 67P. Nonetheless, this dust-only spectrum of 67P represents a good tool to help for the analysis of future comet observations with MUSE.

The morphology of the coma, the dust activity, and the dust reflectivity gradient as measured with MUSE are consistent with other measurements performed with ground-based telescopes at the same epoch as well as with in-situ measurements from the Rosetta mission. About two weeks before our observations, on February 2 2016, an outburst was detected by several instruments on-board the Rosetta orbiter (Grün et al. 2016). A similar but less intense event was also detected in NAVCAM images on 1 March 2016. Since the outburst(s) happened shortly before our observations, we could have expected to detect some indications of such event, either in terms of activity level or coma morphology, in our observations. However, we do not detect any sign of the outburst in the MUSE data presented here. In Fig. 4, we show images obtained with the OSIRIS WAC on-board Rosetta. In the left part, gas jets are easily detected during the February 19 outburst. The right side shows the comet at the same epoch as our observations. No gas jets are visible. This is consistent with the comet being back to a quieter state and the fact that we do not detect the outburst in the MUSE observations. Through careful reprocessing and separating comet and sky signals spatially (rather than by resolving them spectroscopically), we detect the 630 nm forbidden oxygen line and derive a water production rate of 1.5 ± 0.7 × 10 26 molec./s in the coma of 67P, if all [OI] atoms in ( 1 D) state are produced by water. This value is consistent with the Rosetta measurements (Hansen et al. 2016;Combi et al. 2020;[START_REF] Laeuter | [END_REF]Biver et al. 2019). It is the only measurement of the water production rate of 67P at such large heliocentric distance from remote observations. Ground-based detection of water (or water products) in the coma of 67P only happened close to perihelion (Snodgrass et al. 2017b), so that this point is precious to compare measurements and models from the Rosetta spacecraft to ground-based measurements over a larger portion of the comet orbit. In general, detecting water at such low production rate for comets at 2.5 au from the Sun is extremely difficult (see e.g. discussion on water detection in Snodgrass et al. (2017a)). MUSE had been shown in the past to have a huge potential to study species parentage in the coma of brighter comets (Opitom et al. 2019). Our present work demonstrate the efficiency of MUSE to detect low level of water production around distant solar system bodies. This opens up future opportunities, at a time where an intensive search for evidence of water ice in the main asteroid belt and elsewhere in the Solar System is being performed.

Fig. 1 .

 1 Fig. 1. Top: Spectrum of 67P dust coma (black), together with the solar spectrum used to compute the reflectance spectrum (red). The blue dotted line is the spectrum extracted over a 10 aperture. The solar spectrum has been shifted arbitrarily over the y-axis for better visibility. Bottom: Relative reflectance spectrum of 67P (black), compared to the one measured with X-Shooter in November 2014 (red; Snodgrass et al. 2016).

Fig. 2 .

 2 Fig. 2. Top row: combined maps of 67P in V, R, and I bands. The maps are centred on the comet and the FoV is 1 ×1 . Bottom row: V, R, and I maps enhanced by dividing by an azimuthal median profile (Samarasinha et al. 2013).

Fig. 3 .

 3 Fig. 3. Combined maps of the sky+comet flux at the wavelengths of the 3 forbidden oxygen lines at 557.7 nm (left), 630 nm (middle), and 636.4 nm (right). The FoV is 52 ×52 and is oriented North up and East left and the black arrow in the centre points to the position of the comet optocenter. The unit of the color bars is 10 -20 erg/s/cm 2 .

Table 1 .

 1 Observing circumstances of the 67P MUSE campaign

	Article number, page 1 of 7

Acknowledgements. Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme 096.C-0160(A). Datasets of the ESA/Rosetta OSIRIS instrument have been downloaded from the ESA Planetary Science Archive. The authors acknowledge the OSIRIS Principal Investigator H. Sierks (MPS, Goettingen, Germany) and the OSIRIS Team for providing images and related datasets and the ESA Rosetta Project for enabling the science of the mission. -Article number, page 6 of 7 C. Opitom et al.: 67P MUSE