Gaël Servignat

Jérôme Novak

Isabel Cordero-Carrión

A new formulation of general-relativistic hydrodynamic equations using primitive variables

Keywords: hydrodynamics, primitive variables, general relativity, 3+1 formalism, pseudospectral methods, neutron stars

We present the derivation of hydrodynamical equations for a perfect fluid in General Relativity, within the 3+1 decomposition of spacetime framework, using only primitive variables. Primitive variables are opposed to conserved variables, as defined in the widely used Valencia formulation of the same hydrodynamical equations. The equations are derived in a covariant way, so that they can be used to describe any configuration of the perfect fluid. Once derived, the equations are tested numerically. We implement them in an evolution code for spherically symmetric self-gravitating compact objects. The code uses pseudospectral methods for both the metric and the hydrodynamics. First, convergence tests are performed, then the frequencies of radial modes of polytropes are recovered with and without the Cowling approximation, and finally the performance of our code in the black hole collapse and migration tests are described. The results of the tests and the comparison with a reference corecollapse and neutron star oscillations code suggests that not only our code can handle very strong gravitational fields, but also that this new formulation helps gaining a significant amount of computational time in hydrodynamical simulations of smooth flows in General Relativity.

Introduction

General-relativistic hydrodynamics is a key ingredient in many astrophysical models, most often combined with a numerical approach to obtain a solution (see [START_REF] José | Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity[END_REF] and [START_REF] Rezzolla | Relativistic Hydrodynamics[END_REF] for reviews on the subject). As examples, one can think of several phenomena in the field of high-energy astrophysics, such as binary neutron star (NS) mergers [START_REF] Baiotti | Binary neutron star mergers: a review of Einstein's richest laboratory[END_REF], NS oscillations [START_REF] Krüger | Fast Rotating Relativistic Stars: Spectra and Stability without Approximation[END_REF], core-collapse supernovae [START_REF] Janka | Physics of Core-Collapse Supernovae in Three Dimensions: A Sneak Preview[END_REF] or accretion flows around black holes [START_REF] Marek | Foundations of Black Hole Accretion Disk Theory[END_REF]. These simulations have undergone tremendous progress in the last two decades, particularly since the first three-dimensional fully-relativistic code has simulated the merger of a binary NS system [START_REF] Shibata | Simulation of merging binary neutron stars in full general relativity: Γ=2 case[END_REF]. Later major improvements have been obtained in the numerical modeling thanks to the development of so-called upwind high-resolution shock-capturing schemes written in conservation form and based on the characteristic fields of the system of hydrodynamics equations [START_REF] Banyuls | Numerical {3 + 1} General Relativistic Hydrodynamics: A Local Characteristic Approach[END_REF]. This type of methods possesses the great advantage of being able to sharply and stably resolve shocks and is thus able to model ultrarelativistic flows. However, the fact of using conservation schemes means that time evolution is done with conserved variables, which are different from primitive variables (velocity and thermodynamic quantities). Therefore, it is necessary, at each time-step, to pass from conserved to primitive variables with the solution of a set of equations and the call to the equation of state (EoS). Although some optimization is possible [START_REF] Noble | Type II critical phenomena of neutron star collapse[END_REF], this step requires a non-negligible amount of computing time and may be the source of code failure, in particular if the EoS is not given by an analytic expression. Note that Lagrangian methods often rely on this recovery in their algorithm too [START_REF] Rosswog | SPHINCS BSSN: a general relativistic smooth particle hydrodynamics code for dynamical spacetimes[END_REF].

In contrast to the ultrarelativistic cases mentioned above, numerical studies of subsonic flows may not require such sophisticated approaches and, in order to speed-up the numerical process, an interesting alternative can be to evolve directly the primitive variables, not considering the hydrodynamic equations in a conserved form. By doing so, the recovery step is skipped, and the potential problems associated with it are avoided. The aim of this article is to derive such explicit hydrodynamics equations for primitive variables, within the so-called 3+1 formalism of General Relativity, and to test them in some simplified setting (one-dimensional fluid with analytic EoS) in the simulation of the evolution of a self-gravitating compact star. The paper is organized as follows: in Section 2, we introduce the fluid variable notations and thermodynamic properties, and the general four-dimensional covariant equations; Section 3 is devoted to the computation of hydrodynamic equations in their 3+1 form, using explicitly primitive variables; these equations are then implemented in a numerical code, presented in Section 4 with which several tests are performed to compare the new formulation and code with previous ones, and are presented in Section 5. Finally, Section 6 summarizes the results and gives some concluding remarks. In the following, we use units such that c = 1, a 4-metric g µν with a signature (-, +, +, +), and a covariant derivative (connection) ∇ µ associated to that metric. Greek indices are running from 0 to 3, whereas Latin lowercase ones are running from 1 to 3. Latin uppercase indices are used to denote various thermodynamic species present in the fluid.

Fluid variables and covariant conservation laws

Notations and definitions

We consider a perfect fluid, for which the energy-momentum tensor is:

T µν fluid = (e + p)u µ u ν + pg µν , (1)
where u µ is the timelike unitary 4-velocity that carries all species and that satisfies u µ u µ = -1. e is the energy density in the fluid frame, and p the pressure. The thermodynamics of the fluid is described with the following parameters: the number densities n X where X spans all the types of species in the fluid, with the particular case of the subscript B that refers to the baryons, and the entropy density S which is conjugate to the temperature θ:

θ = θ(S, n X ) = ∂e(S, n X ) ∂S n X . (2)
We shall denote by Ŝ = S/n B the entropy per baryon and Y A = n A /n B the abundancies, where A is a subscript that refers to any species except the baryons. Finally let

µ Z = µ Z (S, n X ) = ∂e(S, n X ) ∂n Z S,n Y =Z (3)
be the chemical potential associated with the species Z.

Linking the thermodynamic variables

The differential of the e(S, n X ) is:

de = θ dS + Z µ Z dn Z . (4)
The relation for the pressure p as a function of S and n X reads:

p = p(S, n X ) = θS + Z µ Z n Z -e. (5)
The differential of p(S, n X ) is:

dp = ∂p ∂S n X dS + Z ∂p ∂n Z S,n Y =Z dn Z . ( 6
)
The definition of the squared sound velocity is given e.g. by [START_REF] Typel | CompOSE CompStar online supernova equations of state harmonising the concert of nuclear physics and astrophysics compose[END_REF] :

c 2 s = dp de ds=dN B =dN A =0 , (7)
where s = S V , N X = n X V are the entropy and number of X particles, with V the volume. Eq. ( 7) can be rewritten as:

c 2 s = 1 e + p ∂p ∂S n X S + Z ∂p ∂n Z S,n Y =Z n Z . (8)

Relativistic conservation laws

The conservation of energy and momentum is written with the following Ansatz for the stress-energy tensor:

∇ µ (T µν fluid + T µν emission ) = 0, (9)
where T µν fluid is defined by Eq. ( 1) and T µν emission encodes all possible emissions such as electromagnetic radiation or neutrino emission via weak reactions. The latter can be incorporated as bulk viscosity [START_REF] Camelio | Simulating bulk viscosity in neutron stars I: formalism[END_REF]. The equation is then rewritten as:

∇ µ T µν fluid = Q ν , (10)
denoting

Q ν = -∇ µ T µν emission .
The conservation of a species X flux is, with a possible source term σ X that account for possible processes such as chemical or nuclear reactions:

∇ µ (n X u µ ) = σ X . (11)
Projecting Eq. ( 9) along u ν gives, also denoting

Q = u ν Q ν : θ ∇ µ (S u µ ) = -Q - Z µ Z σ Z . (12)
Then, in the case θ = 0:

∇ µ (S u µ ) = - 1 θ Q + Z µ Z σ Z , (13)
and in the case θ = 0 the entropy density in the fluid is also 0 as stated by the third law of thermodynamics, and the evolution equation for S is irrelevant.

Substitution of pressure gradients with the log-enthalpy

As we shall see in Secs. 3.6 and 3.7, Eqs. ( 32), [START_REF] Lin | A new spectral apparent horizon finder for 3D numerical relativity[END_REF], (34) and (44) involve the ratio between the gradient of p and e + p, which is numerically ill-defined for applications to self-gravitating objects, as both terms can vanish near the border of such an object. Let H be the log-enthalpy:

H = ln e + p m B n B , (14)
with m B a baryon mass. Then:

dp e + p = dH -e -H θd Ŝ + A =B µ A dY A , (15)
where B denotes the baryons. We thus see that the problematic terms can be replaced with a well-behaved expression that does not involve any ratio of vanishing quantities.

Hydrodynamic equations in 3+1 form

3+1 decomposition

In this subsection, we introduce the notations that we use within the 3+1 formalism framework. For a full introduction to 3+1 formalism and 3+1 hydrodynamics, the reader is referred to [START_REF] Gourgoulhon | 3+1 Formalism in General Relativity[END_REF] and references therein.

We consider an asymptotically flat spacetime. The spacetime metric is given in the standard 3+1 form:

g µν dx µ dx ν := -N 2 dt 2 + γ ij (dx i + β i dt)(dx j + β j dt), ( 16
)
where N is the lapse function, β i the shift 3-vector and γ ij the induced 3-metric on spacelike 3-hypersurfaces, associated with its covariant derivative D i . The extrinsic curvature tensor is defined by

K ij := - 1 2 L n γ ij , (17)
where L n is the Lie derivative along the vector n µ , a unitary, future-oriented, timelike 4-vector normal to a hypersurface Σ t of constant time t:

n µ := N ∇ µ t. (18)
K := K i i denotes the trace of the extrinsic curvature tensor. Finally, we consider the following hydrodynamical quantities: the Lorentz factor with respect to the Eulerian observer will be denoted Γ, the Eulerian velocity U µ and the coordinate velocity v i . These quantities are related according to:

Γ := -n µ u µ = (1 -U i U i ) -1/2 , (19)
and

u α = Γ(n α + U α ) = Γ N (1, v i ), U i = 1 N (v i + β i ). (20)

3+1 formulation of conservation laws

In this subsection we write the conservation laws [START_REF] Typel | CompOSE CompStar online supernova equations of state harmonising the concert of nuclear physics and astrophysics compose[END_REF] and [START_REF] Gourgoulhon | 3+1 Formalism in General Relativity[END_REF] in the framework of the 3+1 formalism, and write the 3+1 relativistic equivalent of Euler's equation. Rewriting Eq. ( 13) by reformulating the divergence with γ = det(γ ij ), the determinant of the 3-metric:

∇ µ (S u µ ) = 1 N √ γ ∂ µ (N √ γ S u µ ) = - 1 θ Q + Z µ Z σ Z . (21)
Using Γ = N u 0 and then the 3-dimensional divergence formula, the previous equation is written as:

∂ t (S Γ √ γ) + √ γD i (S Γ v i ) = - N √ γ θ Q + Z µ Z σ Z . (22)
Expanding the time derivative we get:

∂ t S = -S ∂ t Γ Γ -S ∂ t γ 2γ - D i (S Γ v i ) Γ - N Γθ Q + Z µ Z σ Z . (23)
Then we can deduce analogously that:

∂ t n X = -n X ∂ t Γ Γ -n X ∂ t γ 2γ - D i (n X Γv i ) Γ + N Γ σ X . (24)
The 3+1 Euler's equation is [START_REF] Gourgoulhon | 3+1 Formalism in General Relativity[END_REF]:

∂ t U i + v j D j U i = - 1 Γ 2 (e + p) N D i p + U i (∂ t p -β j D j p) + U j D i β j -D i N + U i U j (D j N -N K jk U k ). ( 25
)
As we can see, the evolution of the parameters is governed by the evolution of the Lorentz factor, the determinant of the 3-metric and the pressure. The following subsections will be dedicated to the derivation of the evolution equations for those three quantities. At the end, the equations will be decoupled to have evolution equations for S, n X and U i that only rely on spatial derivatives of the thermodynamical and hydrodynamical primitive quantities.

Evolution equation for γ

The evolution equation for the 3-metric combined with the derivative of a determinant yields [START_REF] Gourgoulhon | 3+1 Formalism in General Relativity[END_REF]:

∂ t γ 2γ = -N K + D i β i . (26)

Evolution equation for p

Starting from Eq. ( 6) and using Eqs. ( 23), ( 24) and ( 8):

∂ t p = ∂p ∂S n X ∂ t S + Z ∂p ∂n Z S,n Y =Z ∂ t n Z = -(e + p)c 2 s ∂ t Γ Γ + ∂ t γ 2γ + D i (Γv i ) Γ -v i D i p + N Γ - ∂p ∂S n X 1 θ Q + Z µ Z σ Z + Z ∂p ∂n Z S,n Y =Z σ Z . (27)
The last term in the equation is computed from the EoS and the microphysical reaction rates. In the following, it will be denoted with Ξ:

Ξ = - ∂p ∂S n X 1 θ Q + Z µ Z σ Z + Z ∂p ∂n Z S,n Y =Z σ Z . (28)
Note that in the special case of no emission (Q = 0) and no chemical or nuclear reactions (σ X = 0, ∀X), Ξ = 0.

Evolution equation for Γ

The evolution equation for the Lorentz factor is deduced from its relation to U i :

∂ t Γ = 1 2 Γ 3 ∂ t (γ ij U i U j ) = 1 2 Γ 3 (U i U j ∂ t γ ij + 2U i ∂ t U i ). ( 29
)
Injecting Eqs. ( 25) and ( 27), as well as the evolution equation for the inverse metric [START_REF] Gourgoulhon | 3+1 Formalism in General Relativity[END_REF] in the previous equation, gives:

∂ t Γ = Γ(N U k U l K kl -U j D j N ) -Γ 3 U i v j D j U i - ΓU i e + p N D i p + U i -v j D j p + N Γ Ξ -β j D j p -(e + p)c 2 s ∂ t Γ Γ + ∂ t γ 2γ + D i (Γv i ) Γ . ( 30
)
Bringing all the ∂ t Γ terms together on the left hand side finally yields:

∂ t Γ Γ = Γ 2 Γ 2 -c 2 s (Γ 2 -1) N U i U j K ij -U i D i N -Γ 2 U i v j D j U i - Γ 2 -1 Γ 2 (e + p) N Γ Ξ - 1 Γ 2 (e + p) N U j D j p + c 2 s Γ 2 -1 Γ 2 ∂ t γ 2γ + D i (Γv i ) Γ . (31)
3.6. Evolution equations for S, n X and U i Replacing Eqs. ( 26) and [START_REF] Thierfelder | Numerical relativity simulations of binary neutron stars[END_REF] in Eqs. ( 23) and [START_REF] Lin | Rotating star initial data for a constrained scheme in numerical relativity[END_REF] gives the evolution equations for S and n X :

∂ t S = -v i D i S - SN Γ 2 Γ 2 -c 2 s (Γ 2 -1) U i U j K ij -K - U j D j p Γ 2 (e + p) + D j U j - Γ 2 -1 Γ 2 (e + p) Ξ Γ - N Γ 1 θ Q + Z µ Z σ Z , (32)
∂ t n X = -v i D i n X - n X N Γ 2 Γ 2 -c 2 s (Γ 2 -1) U i U j K ij -K - U j D j p Γ 2 (e + p) + D j U j - Γ 2 -1 Γ 2 (e + p) Ξ Γ + N Γ σ X . (33)
For the evolution equation for U i , we use Eqs. ( 27), ( 26) and (30):

∂ t U i = -v j D j U i + U j D i β j -D i N + U i U j D j N + c 2 s N U i Γ 2 -c 2 s (Γ 2 -1) (D j U j -K) + U i Γ 2 (c 2 s -1) Γ 2 -c 2 s (Γ 2 -1) N U l U j K lj - N Γ 2 D i p e + p - Γ 2 (1 -c 2 s ) Γ 2 -c 2 s (Γ 2 -1) U i U j D j p e + p - U i Γ 2 -c 2 s (Γ 2 -1) N Ξ Γ(e + p) . ( 34
)
Note that the causality condition is c 2 s ≤ 1, and Γ 2 ≥ 1 by definition, thus Γ 2c 2 s (Γ 2 -1) ≥ 1, and therefore the denominator Γ 2c 2 s (Γ 2 -1) never vanishes.

These equations can be analytically compared with some particular cases: we checked that in the Newtonian limit, in a non reactive system (i.e. with Q = 0 and σ X = 0, ∀X), Eqs. ( 33) and (34) yield their Newtonian counterparts, and that using the polar slicing and radial gauge, these equations exactly correspond to that of [START_REF] Gourgoulhon | Simple equations for general relativistic hydrodynamics in spherical symmetry applied to neutron star collapse[END_REF].

Set of evolution equations for numerical simulations

Depending on the different equilibrium conditions, the EoS depends on a certain number of thermodynamical variables, that we call parameters. For practical applications, the number of parameters in the EoS is usually between one and three. The choice of three parameters typically corresponds to the description of a hot fluid where weak nuclear β reactions take place. The three parameters are usually (θ, n B , Y e ), or any thermodynamically equivalent quantity. Y e corresponds to the electron fraction in the fluid. However, the density of baryons may be discontinuous across a perfect fluid, for example in the presence of a phase transition, whereas the log-enthalpy is betterbehaved. For applications using pseudospectral methods, see Sec. 4.2, we use ( Ŝ, H, Y e ), where those parameters were defined in Secs. 2.1 and 2.4. The source term for the baryon density is σ B := 0. It encodes that the baryon number is strictly conserved. The conservation of the lepton number is:

∇ µ (n l u µ ) = 0. (35)
Therefore, the associated chemical potential is that of the leptons:

µ l := ∂e ∂n l S,n B . (36)
However, we assume that the fluid is transparent to neutrinos so that they escape instantly. Under this particular assumption, with n l = n e + n νe , the energy density does not contain the contribution from the neutrinos, and the lepton chemical potential expression is rewritten as:

µ l = ∂e ∂n e S,n B , (37)
and denoting σ = σ e := -∇ µ (n νe u µ ) the conservation law becomes:

∇ µ (n e u µ ) = σ. (
Then the evolution equation for Y e is simply:

∂ t Y e = -v i D i Y e + N Γ σ n B . (40)
Then, for the entropy per baryon, we relate the differentials of S and Ŝ:

d Ŝ = 1 n B dS -Ŝdn B . (41)
The evolution equation for Ŝ is then once again a simple advection equation:

∂ t Ŝ = -v i D i Ŝ - N Γθn B (Q + µ l σ) . (42)
Those two equations are the same as the ones of [START_REF] Pascal | Proto-neutron star evolution with improved charged-current neutrino-nucleon interactions[END_REF], that were also derived by [START_REF] Pons | Evolution of Proto-Neutron Stars[END_REF].

The corresponding differential of H in the particular case of the parameters ( Ŝ, H, Y e ) that we chose is:

dH = dp e + p + e -H (θd Ŝ + µ l dY e ), (43)
Therefore, combining Eqs. ( 27), ( 40) and (42), the evolution equation for the logenthalpy is:

∂ t H = -v i D i H -c 2 s Γ 2 N Γ 2 -c 2 s (Γ 2 -1) K j U i U j -K + D i U i - U i Γ 2 D i H -e -H θD i Ŝ + µ l D i Y e + e -H N Γ Q n B + e -H N Γ (Γ 2 -1)c 2 s Γ 2 -c 2 s (Γ 2 -1) ∂p ∂n e S,n B σ n B - ∂p ∂S ne,n B Q + µ l σ n B θ . (44)
The final set of variables that are evolved through partial differential equations are (U i , H, Y e , Ŝ), for which the corresponding equations are Eqs. (34), (44), ( 40) and ( 42), where all the terms of the form D i p/(e + p) have been replaced thanks to Eq. [START_REF] Pascal | Proto-neutron star evolution with improved charged-current neutrino-nucleon interactions[END_REF]. The system must be closed by providing an EoS, which links the pressure, energy density, sound speed, chemical potential and temperature to the parameters (H, Y e , Ŝ). In order to check the validity of the set of equations, we will perform numerical tests that will be described in the following sections.

Description of the code

Conformal decomposition and choice of coordinates

In order to fix the coordinates, we must choose a gauge and a foliation. Following [START_REF] Bonazzola | Constrained scheme for the Einstein equations based on the Dirac gauge and spherical coordinates[END_REF], we use the conformal decomposition of spacelike hypersurfaces. Let f ij be a flat metric, Di its associated covariant derivative, f = det(f ij ). We require that γ ij is f ij at spatial infinity. By introducing the so-called conformal factor

Ψ = γ f 1/12 , (45)
the conformal metric is defined as

γij = Ψ -4 γ ij . (46)
We choose the maximal slicing as a foliation, mathematically expressed as:

K = 0. ( 47
)
This choice has the effect of simplifying some of the equations and it numerically avoids possible black hole singularities, for example during a gravitational collapse, see Sec. 5.3. We also choose the Dirac gauge:

Di γij = 0. ( 48
)
In spherical symmetry, it corresponds to the isotropic gauge and γij = f ij . The Dirac gauge is straightforward to use in three dimensions and allows to compute quasistationary initial data [START_REF] Bonazzola | Constrained scheme for the Einstein equations based on the Dirac gauge and spherical coordinates[END_REF].

Numerical implementation in spherical symmetry

We will perform various numerical tests on a one-species, spherically symmetric, cold, catalysed (i.e. where β reactions are at equilibrium) NS, which is described by a oneparameter EoS. We choose the parameter to be the log-enthalpy H. In this case, the EoS only enters through the sound speed squared c 2 s (H). In particular, Ξ = Q = σ = 0 in Eqs. ( 34) and (44). The evolution equations are that of U r (the first component of U i in a spherical triad) (34) and H (44). We implement them along with the spherically symmetric version of the Einstein system of [START_REF] Cordero-Carrión | Improved constrained scheme for the Einstein equations: An approach to the uniqueness issue[END_REF] (Eqs. ( 30) to (33) therein), which is a fully constrained formulation of Einstein equations, formulated as a set of four coupled Poisson-like partial differential equations. The algorithm to perform the time evolution is the following:

• At a given time t, H and U r are known thanks to their evolution equations ( 44) and (34). Γ can be computed from U r , and along with the EoS, H allows to compute p and e. Ψ is temporarily known thanks to the following evolution equation, which can be obtained from Eq. ( 26):

∂ t ln Ψ = β i Di ln Ψ + 1 6 Di β i . (49)
It is then possible to compute the following quantities,

E * = Ψ 6 (Γ 2 (e + p) -p), (50)
S * = Ψ 6 (3p + (Γ 2 -1)(e + p)), (51)
p * r = Ψ 6 Γ 2 (e + p)U r , (52)
and solve the Einstein equations. The output is the value of Ψ consistent with the Hamiltonian constraint, which replaces the value of Ψ computed with Eq. ( 49) as soon as it is available, as well as the lapse N and the radial component of the shift vector β r .

• We use the knowledge of the radius of the star R, as well as H, U r , Ψ, N and β r together with the EoS to compute the sources of Eqs. (34), ( 44), ( 56) and ( 49). An explicit finite-difference time scheme allows to compute Ψ, R, H and U r at the next time step.

The implementation is based upon LORENE [START_REF] Gourgoulhon | LORENE: Spectral methods differential equations solver[END_REF], a dedicated C ++ library that comes with built-in Poisson solvers, based on pseudospectral methods [START_REF] Grandclément | A Multidomain Spectral Method for Scalar and Vectorial Poisson Equations with Noncompact Sources[END_REF][START_REF] Gottlieb | Numerical analysis of spectral methods: theory and applications[END_REF][START_REF] Boyd | Chebyshev and Fourier spectral methods[END_REF][START_REF] Grandclément | Spectral Methods for Numerical Relativity[END_REF].

We use the relaxation iterative method to solve the non-linear elliptic partial differential equations. Pseudospectral methods are also used to compute the spatial derivatives in the source of the evolution equations. The time evolution is done thanks to an explicit finite differences Adams-Bashforth scheme of order 3. The numerical grid comprises two domains: a nucleus centered on the star, with the domain boundary corresponding to the surface of the star, and a compactified external domain (CED) from the surface of the star to infinity, allowing to solve Poisson-like equations and imposing boundary conditions at spatial infinity. As the star is evolved, the boundary between the nucleus and the CED will move with it. This semi-Lagrangian approach was used in [START_REF] Gourgoulhon | Simple equations for general relativistic hydrodynamics in spherical symmetry applied to neutron star collapse[END_REF]. In order to do so, we perform a change of variables in the nucleus defined as follows:

t = t, ξ = r/R(t). ( 53
)
The radial coordinate is mapped from [0, R(t)] to the numerical coordinate ξ ∈ [0, 1]. This changes the differential operators with respect to r and t:

∂ ∂r = 1 R ∂ ∂ξ , (54)
∂ ∂t = - Ṙ R ξ ∂ ∂ξ + ∂ ∂t , ( 55
)
where Ṙ is the time derivative of the radius. The term -Ṙ ξ can be interpreted as a grid advection velocity inside the star. The radius is an additional parameter, therefore we have to introduce a new equation:

Ṙ(t ) := dR dt = v r (R(t ), t ). ( 56
)
It corresponds to the impermeable boundary condition, which requires that no matter crosses the external border of the domain. In the CED, a similar change of variables is performed:

t = t, ξ = 1 -2R(t)/r. ( 57
)
The radial coordinate is mapped from [R(t), +∞] to the numerical coordinate ξ ∈ [-1, 1]. The differential operators are then also changed accordingly:

∂ ∂r = (ξ -1) 2 2R ∂ ∂ξ , (58)
∂ ∂t = Ṙ R (ξ -1) ∂ ∂ξ + ∂ ∂t . ( 59
)
In what follows, we make no distinction between t and t as they are identical. The boundary conditions between the domains for the Poisson-like equations are imposed by the Poisson solvers of LORENE: the continuity of the metric and its radial derivative is required. At infinity, we impose that Ψ → 1, N → 1, β r → 0, as the metric should be asymptotically flat. For the evolution equations, the boundary conditions are the following:

• At the center of the star, i.e. ξ = 0 in the nucleus, we impose the regularity of the fields.

• Eqs. ( 44) and (34), within the assumption of subsonic flow, need one boundary condition at the border of the star, i.e. ξ = 1 in the nucleus. The impermeable boundary condition is naturally imposed by choosing a comoving grid. It also helps with the stability of the code to enforce that H is constant at the surface of the star.

• The communication between the nucleus and the CED for Eq. ( 49) is done in an upwind fashion. We refer to Appendix A for the details of the characteristic speed computation. No boundary condition needs to be imposed at ξ = 1 in the CED as the characteristic speed there is β r (r = +∞, t) = 0.

All initial data are computed with LORENE, using the solver described in [START_REF] Lin | Rotating star initial data for a constrained scheme in numerical relativity[END_REF] in spherical symmetry.

Numerical tests

Convergence

Here we describe the convergence tests that we performed. We use a polytropic EoS with parameters γ = 2, κ = 100 (in geometrized units where G = c = 1, supplemented with M = 1) that was used by [START_REF] José | Three-dimensional numerical general relativistic hydrodynamics. II. Long-term dynamics of single relativistic stars[END_REF] to perform tests, with a central log-enthalpy of H c = 0.2279. The convergence is measured on the conservation of the ADM and baryon masses. First, we fix the timestep ∆t = 3.34×10 -4 ms and vary N r the number of radial grid points between 5 and 65. We let the star at equilibrium (meaning that it oscillates only because of numerical noise) for a 1ms simulation. We find exponential convergence on the maximum of the relative error on the ADM mass with respect to the equilibrium value, as expected with pseudospectral methods. The error reaches a plateau around 10 -10 at 25 grid points. To check the convergence of the time scheme, we use the same star, but this time the simulations were run up to 10 ms, we set N r = 17, the Cowling approximation was used, and a Gaussian enthalpy profile was added to the equilibrium, initial enthalpy profile as a perturbation. The timestep was varied between 10 -4 and 3.34 × 10 -3 ms. We find an order 3 power-law convergence on the relative error of the baryon mass with respect to the initial equilibrium baryon mass, which is expected as the error for a finite difference scheme of order n must behave as O(∆t n ).

Frequency extraction

We seek to recover the frequencies presented elsewhere in the literature. We use the same star than in Sec. 5.1 and the γ = 5/3, κ = 7.308 km 4/3 EoS described in [START_REF] Hartle | Slowly Rotating Relativistic Stars. VIII. Frequencies of the Quasi-Radial Modes of an N = 3/2 Polytrope[END_REF]. The method to extract the frequencies is the following:

(i) A simulation is run up to some t max , set to 1s in this test.

Mass [M ⊙ ] Fundamental mode frequency [kHz]

γ = 3 γ = 2.5 γ = 2 γ = 1.8
Figure 2. Mass-frequency diagrams for four different polytropes. The κ coefficients are 6.23 × 10 6 (γ = 1.8), 4.79 × 10 8 (γ = 2), 1.97 × 10 13 (γ = 2.5) and 6.78 × 10 17 (γ = 3), in units of m 2(γ-1) (in geometrized units i.e. with G = c = 1). Table 1. Comparing polytropic frequencies (fundamental mode, first and second overtone) of the code with the ones found in the literature. The frequencies were extracted from a 1 s simulation. For the first polytrope [START_REF] José | Three-dimensional numerical general relativistic hydrodynamics. II. Long-term dynamics of single relativistic stars[END_REF], κ is expressed in units of G 3 M 2 /c 4 , and for the second [START_REF] Hartle | Slowly Rotating Relativistic Stars. VIII. Frequencies of the Quasi-Radial Modes of an N = 3/2 Polytrope[END_REF], κ is expressed in km 4/3 , i.e. in geometrized units where G = c = 1. The central energy density of the polytrope from [START_REF] Hartle | Slowly Rotating Relativistic Stars. VIII. Frequencies of the Quasi-Radial Modes of an N = 3/2 Polytrope[END_REF] is e c = 3.16 × 10 14 g • cm -3 . The corresponding spectrum for the star of [START_REF] José | Three-dimensional numerical general relativistic hydrodynamics. II. Long-term dynamics of single relativistic stars[END_REF] is given on the left panel of Fig. 1. (ii) A Fast Fourier Transform algorithm is performed on the R(t) signal to compute its spectrum R(f ).

κ γ H c [c 2 ] M [M ] R [km] Fund
(iii) A local maximum research is performed on R(f ). When a peak (i.e. a maximum) is detected, the point where the maximum is reached, along with its two neighbours, are interpolated with a quadratic Lagrange polynomial.

(iv) The maximum of the polynomial can be analytically extracted; it is an approximation of the frequency of the peak.

The precision on the frequency extraction is less than 1/t max . Note that we compute coordinate frequencies, as the time considered for oscillations is the coordinate time, not the proper time. For the frequencies computed here, we used N r = 17 and ∆t = 3.34×10 -4 ms. The results are reported in Table 1. The code used in [START_REF] José | Three-dimensional numerical general relativistic hydrodynamics. II. Long-term dynamics of single relativistic stars[END_REF] is a full GR 3D code with Eulerian high resolution shock capturing methods. The frequencies were also recovered by [START_REF] Hébert | General-relativistic neutron star evolutions with the discontinuous Galerkin method[END_REF] with discontinuous Galerkin methods and [START_REF] Rosswog | SPHINCS BSSN: a general relativistic smooth particle hydrodynamics code for dynamical spacetimes[END_REF] with Lagrangian Smooth Particles Hydrodynamics methods. We recover all the frequencies with a precision of less than 1%. We also look at the frequencies of the [START_REF] José | Three-dimensional numerical general relativistic hydrodynamics. II. Long-term dynamics of single relativistic stars[END_REF] polytrope in the Cowling approximation. The results are reported in Table 2. Here we once again recover the tabulated frequencies with an accuracy of less than 1%. Fig. 1 shows the spectra for our simulation both in full metric evolution and using the Cowling approximation. On Fig. 2 are plotted mass-frequency diagrams for four different polytropes of which the κ parameter has been adjusted to yield the same maximum mass, of 2.43M . To produce such a diagram, we extract frequencies on several stars described by the same EoS, from close to zero mass to maximum mass, by varying the central enthalpy. We see that, as expected [START_REF] Harrison | Gravitation theory and gravitational collapse[END_REF], the frequency of the fundamental mode of the stars drops down to zero close to the maximum mass. We use the occasion to compare our code with the reference core-collapse and NS oscillation code CoCoNuT [START_REF] Dimmelmeier | Combining spectral and shock-capturing methods: A new numerical approach for 3D relativistic core collapse simulations[END_REF] which uses pseudospectral methods for the metric equations but which is based upon finite volume methods and a conserved scheme for the hydrodynamics. For the same γ = 2, κ = 100 EoS we used to produce frequencies, a run of the same star, with comparable to worse precision (as measured on the conservation of the ADM mass), is typically five times longer with CoCoNuT than with our code.

Collapse to a black hole

The next test we perform is the collapse of a NS to a black hole. The test is rather demanding but often performed [START_REF] José | Three-dimensional numerical general relativistic hydrodynamics. II. Long-term dynamics of single relativistic stars[END_REF][START_REF] Cordero-Carrión | Improved constrained scheme for the Einstein equations: An approach to the uniqueness issue[END_REF][START_REF] Rosswog | SPHINCS BSSN: a general relativistic smooth particle hydrodynamics code for dynamical spacetimes[END_REF][START_REF] Thierfelder | Numerical relativity simulations of binary neutron stars[END_REF][START_REF] Bernuzzi | Constraint violation in free evolution schemes: Comparing the BSSNOK formulation with a conformal decomposition of the Z4 formulation[END_REF][START_REF] Baiotti | Three-dimensional relativistic simulations of rotating neutron-star collapse to a Kerr black hole[END_REF]. We use the γ = 2, κ = 100 EoS that we used in Secs. 5.1 and 5.2, and set H c = 0.61, on the unstable branch, and perturb it with an enthalpy profile so that it exceeds the maximum mass. The star then collapses into a black hole. We use N r = 129 and ∆t = 3.34 × 10 -6 ms. We see on Fig. 3 that we indeed find a black hole: the value of the central lapse goes as low as a few 10 -6 . We find that the hydrodynamical quantities such as the central pressure, baryon density, logenthalpy, sound speed squared and the circumferential radius seem to freeze; this is due to the maximal slicing condition, which avoids black hole singularities singularities [START_REF] Gourgoulhon | 3+1 Formalism in General Relativity[END_REF].

The last evidence is the existence of an apparent horizon, which we found using the Apparent Horizon Finder described in [START_REF] Lin | A new spectral apparent horizon finder for 3D numerical relativity[END_REF]. It appears around t = 0.18 ms at the corresponding Schwarzschild radius of the star: M ADM = 2M , R AH = 6 km. This demonstrates the ability of the code to handle very strong gravitational fields.

Migration test

The migration test is the final standard test for numerical isolated NS evolution codes that we perform. We try to reproduce the test originally presented by [START_REF] José | Three-dimensional numerical general relativistic hydrodynamics. II. Long-term dynamics of single relativistic stars[END_REF] which has also been done by [START_REF] Cordero-Carrión | Improved constrained scheme for the Einstein equations: An approach to the uniqueness issue[END_REF][START_REF] Rosswog | SPHINCS BSSN: a general relativistic smooth particle hydrodynamics code for dynamical spacetimes[END_REF][START_REF] Thierfelder | Numerical relativity simulations of binary neutron stars[END_REF][START_REF] Bernuzzi | Constraint violation in free evolution schemes: Comparing the BSSNOK formulation with a conformal decomposition of the Z4 formulation[END_REF][START_REF] Baiotti | Three-dimensional relativistic simulations of rotating neutron-star collapse to a Kerr black hole[END_REF]. Unfortunately, as can be seen in the simulation run with CoCoNuT, a shock is formed and we are able to compute the migration only up to the shock formation, as one would expect from using pseudospectral methods. This is illustrated in Fig. 4 where we compare our simulation with a migration test performed with CoCoNuT [START_REF] Dimmelmeier | Combining spectral and shock-capturing methods: A new numerical approach for 3D relativistic core collapse simulations[END_REF]. The simulation is run with N r = 17 and ∆t = 3.34 × 10 -4 . The migration being triggered by numerical truncation error, it has no reason to start at the same time, so we have to shift the time of the one performed with CoCoNuT to align the curves. We see that because our simulation develops Gibbs phenomenon, the star is not able to oscillate. It is worth noting that, when run on a laptop, our codes reaches the point where it crashes in approximately ten seconds, while CoCoNuT takes six minutes to reach that same point.

Summary and conclusions

We have presented a new formulation of the hydrodynamical equations in General Relativity, using the 3+1 formalism, aimed for numerical simulations. Those equations only rely on the primitive variables, meaning that, as opposed to the widely used conservative form of the Valencia hydrodynamics formalism [START_REF] Banyuls | Numerical {3 + 1} General Relativistic Hydrodynamics: A Local Characteristic Approach[END_REF], no recovery procedure is needed to compute the sources for the metric equations. The derivation of the equation relies on a procedure that shows the evolution of the variables n X , S, U i , Γ, γ and p is actually degenerate and decouples them to end up with equations for n X , S and U i only. This formulation allows for the simulation of any relativistic perfect fluid, not only selfgravitating objects, and is expressed in a 3+1 covariant form. Then, we have summarized the equations for a fluid with an EoS that uses at most three parameters. The extension to more parameters, for example accounting for the presence of muons, is straightforward as the new equations needed are already encoded in the general derivation. We then applied them to the simulation of NS oscillations. The code used for numerical tests implements a fully constrained formalism for the computation of the metric. It uses an Adams-Bashforth finite-differences scheme for the time evolution, and pseudospectral methods for the spatial coordinates. Pseudospectral methods are well-suited for smooth flows and help to save even more computational resources. A notable feature is that, like in the Lagrangian approach of [START_REF] Rosswog | SPHINCS BSSN: a general relativistic smooth particle hydrodynamics code for dynamical spacetimes[END_REF], the surface of the star does not need any specific treatment; we use the fact that it is well-defined at equilibrium, and then evolve the grid according to the star evolution. The code is still at an early stage of development, hence only spherically symmetric configurations have been implemented. Nonetheless, the frequencies of polytropic NS are recovered with excellent precision, and the code is also able to perform a black hole collapse. Only the migration test was not fully run because of a shock formation, which also highlights the limitations of our numerical approach. Future work includes the extension of the code to three spatial dimensions, as well as two or three parameters in the EoS.

Figure 1 .

 1 Figure1. Fourier spectrum of the radius over a 1 s simulation of a γ = 2 polytropic NS with full spacetime evolution (left) and Cowling approximation (right). The frequencies of the three first modes are reported in Tables 1 (full spacetime evolution) and 2 (Cowling approximation). The spectrum is given in arbitrary units.

Figure 3 .

 3 Figure 3. BH collapse: lapse function at the center of the star. The collapse of the lapse is characteristic of the formation of a black hole.

Figure 4 .

 4 Figure 4. Comparing our code with CoCoNuT for a migration test. Please see main text for details.

  38)Then, the knowledge of ∂ t S, ∂ t n B and ∂ t n e with Eqs. (32),[START_REF] Lin | A new spectral apparent horizon finder for 3D numerical relativity[END_REF], allows to compute ∂ t Ŝ, ∂ t H and ∂ t Y e from the differential relation between those variables. Let us start with the electron fraction:

	dY e =	1 n B	dn e -	Y e n B	dn B .

Table 2 .

 2 Comparing polytropic frequencies (fundamental mode, first and second overtone) of the code with the ones found in the literature using the Cowling approximation. The frequencies were extracted from a 1 s simulation. κ is expressed in units of G 3 M 2 /c 4 . The corresponding spectrum is given on the right panel of Fig.1.

		κ	γ H c [c 2 ] M [M ] R [km] Fund. [kHz] 1st ov. [kHz] 2nd ov. [kHz]
	Font et al. [25]	100 2 0.2279 1.4	14.15	2.696	4.534	6.346
	This work	100 2 0.2279 1.401	14.16	2.685	4.548	6.339
	Relative difference				0.4%	0.3%	0.1%

Acknowledgments

The authors gratefully acknowledge the Italian Instituto Nazionale de Fisica Nucleare (INFN), the French Centre National de la Recherche Scientifique (CNRS) and the Netherlands Organization for Scientific Research for the construction and operation of the Virgo detector and the creation and support of the EGO consortium. This work has been supported by the Spanish Agencia Estatal de Investigación / Ministerio de Ciencia, Innovación y Universidades through the Grants No. PGC2018-095984-B-

Appendix A. Characteristic speed for the conformal factor evolution Eq. ( 49) is an advection equation, therefore taking the moving grid into account, the characteristic speed is

in the nucleus and

in the CED. At the exterior boundary of the nucleus, the characteristic speed reads

and at the interior boundary of the CED, the characteristic speed is

Therefore, the lapse function being positive, when U r (R(t), t) < 0, the value of ln Ψ at the last grid point of the nucleus is copied to the first grid point of the CED, and vice-versa for U r (R(t), t) > 0.


Tristan Cragnolini

Samuela Pasquali

Philippe Derreumaux

Marc Baaden

UnityMol: Interactive Scientific Visualization for Integrative Biology S ébastien Doutreligne *

Keywords: I, 6, 6 [Simulation and Modeling]-Simulation Output Analysis; J, 3 [Computer Applications]: Life and Medical

A broad challenge facing scientists today is the availability of huge amounts of data from various sources. Computers are required to store, analyze, explore and represent these data in order to extract useful information. With UnityMol, we pursue the ambitious goal to create an interactive virtual laboratory enabling researchers in biology to visualize biomolecular systems, run simulations and interact with physical models and data. Visual effects can enrich the dynamic and immersive aspects. Ultimately, we want to combine an appealing visual feedback already in place with a set of analysis features to extract information about properties of the fascinating biomolecular systems under study.

INTRODUCTION

Nowadays, an increasing number of structures of molecular assemblies is available in a variety of databases. With the advances in experimental structure determinations, the size of these assemblies keeps expanding. Descriptive data from various measurements is abundant. Studying and exploring the high-dimensional set of porperties of such huge structures is part of today's challenges to modern biology.

Visualization of molecular structures is commonly achieved by specialized software. UnityMol [START_REF] Lv | Game on, science -how video game technology may help biologists tackle visualization challenges[END_REF] in particular, uses a raycasting technique called HyperBalls to render evolving atoms and bonds [START_REF] Chavent | GPU-accelerated atom and dynamic bond visualization using hyperballs: a unified algorithm for balls, sticks, and hyperboloids[END_REF]. This innovative method is well suited for visualizing molecules moving according to accurate molecular simulations. In the last months, many features were added to turn UnityMol into an interactive virtual biology laboratory tuned for this purpose. Among them, we can cite the support of haptic interaction with the molecule, rendering on large high-resolution display walls and implementation of custom representations for RNA and sugar molecules.

Here, we present UnityMol's current features and, to some extent, our future developments for the interactive visualization and analysis of molecular and biological data.

MOLECULAR VISUALIZATION IN UNITYMOL

includes both original and more traditional implementations of molecular graphics algorithms based on polygonal static and ray-casted dynamic rendering techniques. The former provide common tools for a static view of the system, the latter elegantly render dynamic biomolecular events such as the formation of bonds * e-mail:doutreligne@ibpc.fr † e-mail:tristan.cragnolini@ibpc.fr ‡ e-mail:samuela.pasquali@ibpc.fr § e-mail:philippe.derreumaux@ibpc.fr ¶ e-mail:baaden@smplinux.de and interactions. Complex and large-scale data may be rendered on advanced displays such as the 33 megapixel 4,3x2,4 m display wall at Maison de la Simulation in Saclay. See Figure 1.

Secondary Structures, SugarRibbons and other Representations

UnityMol makes use of polygons to render molecular features, including among others secondary structures based on splines, molecular isosurfaces determined by marching cubes and a novel polysaccharide representation for complex sugar molecules. These representations are shown in Figure 2.

HyperBalls Representation

UnityMol comes with an HyperBalls GPU implementation. The heart of this ray-casting technique is the use of quadratic surface equations, perfectly compatible with a traditional rendering pipeline. Here, we successfully extend the representation to coarsegrained (CG) HiRE-RNA [START_REF] Pasquali | HiRE-RNA: a high resolution coarsegrained energy model for RNA[END_REF] models. Based on this coarse-grained model developed in the laboratory [START_REF] Cragnolini | Coarse-grained simulations of RNA and DNA duplexes[END_REF], nucleo-bases can be rendered as ellipsoids as shown in Figure 3. Visual feedback of molecular properties can be enhanced using a lit-spheres lighting pass.

LIVE SIMULATION AND INTERACTION

UnityMol has the ability to connect with molecular simulation software via the IMD network protocol, and supports special devices such as an haptic arm and a 3D mouse via VRPN.

Interactive Molecular Dynamics Simulation

UnityMol uses MDDriver, a light and generic framework for coupling molecular simulators with visualization software based on the IMD protocol [START_REF] Delalande | Complex molecular assemblies at hand via interactive simulations[END_REF]. Here, we focus on the coarse-grained, MDDriver-enhanced HiRE-RNA simulator. As a CG model, it provides an excellent compromise between simulation speed and biological fidelity, being in our experience more robust with respect user interactions than computations carried out at an all-atom level. Enriched by special screen-space effects such as motion blur or glowing atoms, the underlying dynamic molecular behavior is highlighted visually.

Virtual reality and Interactive Simulations

The interactive approach opens up perspectives to guide simulations via user input, within a dedicated graphical environment. Peripherals such as a 6 d.o.f. 3D mouse ease the navigation in the virtual space. The user can change the point of view intuitively, and focus on the relevant part. Haptic devices help to render accurate molecular models more real and tangible to the scientists. Hence, the user feels an immediate force feedback by a straightforward combination of molecular modeling and virtual reality leading to an intuitive understanding of the causal relationship between the theoretical model and its biologically relevant properties.

INTERACTIVE UNFOLDING OF AN RNA HAIRPIN

A typical biological application would be to stretch an RNA hairpin on two ends, controlling the successive detachment of its base pairs. When released, the structure may either progressively return to the initial hairpin state or feature a base shift. This numerical experiment provides insight into the hairpin's resilience, similarly to what can be probed experimentally with optical tweezers.

Figure 4 depicts a typical setup. The displayed screen is a montage exhibiting visualization features soon to be implemented, as described in section 5. A benchmark measuring the average FPS is displayed for systems of increasing size.

UPCOMING FEATURES AND FUTURE WORK

Our current focus is on HiRE-RNA specific features such as the detection of hydrogen bonds in realtime. They will then be rendered in 3D as well as in two-dimensional schematic graphs of RNA molecules, commonly used in this scientific community. For this purpose, we couple UnityMol with VARNA, an image generator for RNA secondary structures [START_REF] Darty | VARNA: Interactive drawing and editing of the RNA secondary structure[END_REF]. We intend to add plotting functionalities to monitor the running simulation. We may watch the evolution of the mean distance between atoms, the number of weak bonds, the energy, etc., offering a visual dashboard for investigating molecular properties of biological systems.

Figure 1 :

 1 Figure 1: UnityMol running on a display wall. The image is rendered accross a 4x4 tiled display matrix at 7680 x 4320 pixel resolution. It represents the mono-heme ferrocytochrome c protein depicted in HyperBalls representation decorated with a cracked-stone lit sphere texture in front of a clouded sky skybox as background.

Figure 2 :

 2 Figure 2: Representations of a hydrophobic surface of the ferrocytochrome c protein (A) and its helical secondary structures (B), SugarRibbons representation of the GLC-SIA-FUC-GAL-MAN pentaoligosaccharide (C).

Figure 3 :

 3 Figure 3: A coarse-grained HiRE-RNA model of a 22-Nucleotide Hairpin (1F9L). Left: HyperBall view. Middle: extension with ellipsoid base pairs. Right: textured view.

Figure 4 :

 4 Figure 4: Photo montage of a typical workspace. The user moves the camera (left hand) and manipulates the 1F9L molecule in the center with the haptic arm (right hand). To the right, a 2D representation of the RNA molecule's secondary structure and hydrogen bonds is shown. To the left, a distance time series plot and an energy surface plot are depicted. Bottom right: benchmark for five RNA systems of 146, 303, 459, 475 and 1062 coarse-grain particles, respectively.

CONCLUSION

We presented UnityMol as a platform for a full-blown virtual laboratory, featuring graphical representations of molecular data coupled with interactive visual analysis tools.


Hatem Hajri

Manon Cesaire

Lucas Schott

Sylvain Lamprier

Patrick Gallinari

Manon Césaire

Neural Adversarial Attacks with Random Noises

Keywords: deep neural network classifiers, adversarial attacks, adversarial examples

come   L'archive ouverte pluridisciplinaire

Introduction

The study of adversarial examples in machine learning has been pushed by the desire to boost the performance of deep neural networks and their use in critical applications. Most of the works related to adversarial examples have been around three categories of attacks according to the minimised distance between original and adversarial samples: L 2 (squared error) [START_REF] Szegedy | Intriguing properties of neural networks[END_REF][START_REF] Moosavi-Dezfooli | Deepfool: a simple and accurate method to fool deep neural networks[END_REF][START_REF] Carlini | Towards evaluating the robustness of neural networks[END_REF] , L ∞ (max-norm) [START_REF] Goodfellow | Explaining and harnessing adversarial examples[END_REF][START_REF] Madry | Towards deep learning models resistant to adversarial attacks[END_REF][START_REF] Kurabin | Adversarial examples in the physical world[END_REF] and L 0 (or sparse) [START_REF] Papernot | The limitations of deep learning in adversarial settings[END_REF][START_REF] Carlini | Towards evaluating the robustness of neural networks[END_REF] attacks (minimising the number of modified components).

Given a neural network classifier (NNC) F : R n → R p , the predicted label for an input x is label(x) = argmax k F k (x), where F 1 , • • • , F p are the class probabilities of F . We recall that an adversarial example to x is an item x * such that label(x * ) ̸ = label(x) (untargeted attack), or such that label(x * ) = c, with c ̸ = label(x) a specific class (targeted attack). In this paper, we mostly discuss L 0 (or sparse) attacks and L ∞ ones.

Sparse perturbations can be encountered in many situations and have been motivated in many works [START_REF] Carlini | Towards evaluating the robustness of neural networks[END_REF][START_REF] Papernot | The limitations of deep learning in adversarial settings[END_REF][START_REF] Modas | Sparsefool: a few pixels make a big difference[END_REF][START_REF] Bibi | Analytic expressions for probabilistic moments of pl-dnn with gaussian input[END_REF][START_REF] Alfadly | Network Moments: Extensions and Sparse-Smooth Attacks[END_REF][START_REF] Combey | Probabilistic jacobian-based saliency maps attacks[END_REF][START_REF] Croce | Sparse and imperceivable adversarial attacks[END_REF][START_REF] Croce | Sparse-rs: a versatile framework for query-efficient sparse blackbox adversarial attacks[END_REF][START_REF] Dong | Greedyfool: Distortion-aware sparse adversarial attack[END_REF] . For instance, they could correspond to some raindrops on traffic signs that are sufficient to fool an autonomous driver [START_REF] Modas | Sparsefool: a few pixels make a big difference[END_REF] . Understanding these special alterations is important to mitigate their effects and take a step forward trusting neural networks in real-life.

In this work, we aim at presenting a general probabilistic approach to generate L 0 attacks that rely on random noises; called throughout the paper stochastic sparse adversarial attacks (SSAA). We argue that existing deterministic attacks named XSMA (JSMA [START_REF] Papernot | The limitations of deep learning in adversarial settings[END_REF] , WJSMA, TJSMA [START_REF] Combey | Probabilistic jacobian-based saliency maps attacks[END_REF] ), which perform by sequentially applying maximal perturbations on selected components of the input, fail at reaching accurate adversarial examples on real-world large scale datasets. Fig. 1: Plots of the initial and targeted class probabilities for a one pixel version of XSMA on the left failing to converge along more than 3,000 iterations and our VFGA10 method converging efficiently in less than 400 iterations on the right.

Figure 1 (left) illustrates this failure on the ImageNet dataset [START_REF] Russakovsky | ImageNet Large Scale Visual Recognition Challenge[END_REF] for a onecomponent version of the targeted XSMA which does not succeed to affect the initial probability of the input on the Inception-v3 network [START_REF] Szegedy | Rethinking the inception architecture for computer vision[END_REF] . On the other hand, working with more than one component at a time, while more accurate, does not scale at all on datasets as ImageNet [START_REF] Carlini | Towards evaluating the robustness of neural networks[END_REF] . An alternative would be to repeatedly apply very small perturbations on components, but this would be at the cost of efficiency. These issues, well-known for the XSMA, are completely solved by SSAA through sampling from noise near the most salient component while keeping the same general iterative process followed by the XSMA. Experimental results on large scale datasets, as depicted on the same example as the failure case of XSMA in Figure 1 (on the right), show that our SSAA approaches (denoted VFGA10) succeed at efficiently producing accurate attacks in most cases.

The rest of the paper is organised as follows:

In Section 2, we introduce our best methods of SSAA called scalable SSAA. These attacks are gradient-based and scale efficiently to large datasets.

In Sections 3 and 4, we experiment scalable SSAA on deep NNC on CIFAR-10 [START_REF] Krizhevsky | Cifar-10[END_REF] and ImageNet [START_REF] Russakovsky | ImageNet Large Scale Visual Recognition Challenge[END_REF] and compare their performances with the-state-of-the-art methods SparseFool 20 , GreedyFool 10 , Brendel & Bethge L 0 attack (B&B) [START_REF] Brendel | Accurate, reliable and fast robustness evaluation[END_REF] and Sparse-RS [START_REF] Croce | Sparse-rs: a versatile framework for query-efficient sparse blackbox adversarial attacks[END_REF] . Experimental results show that our method called VFGA outperforms SparseFool while being faster. Furthermore, VFGA achieves better L 0 scores than Sparse-RS in most cases and in all cases on ImageNet; when both attacks are fully successful. Moreover, VFGA is significantly less complex than B&B and GreedyFool and achieves competitive results in many cases.

In Section 5, we introduce more examples of SSAA. These attacks are Hessianbased and so not scalable to ImageNet. However, they can be used as efficient methods on small datasets. We illustrate their usefulness through some experiments on MNIST [START_REF] Lecun | MNIST handwritten digit database[END_REF] .

In Section 6, we give two theoretical results on the expected logits of a ReLU-NNC under a one component perturbation by a Gaussian or Folded Gaussian noise. These results shed some light on the interpretation of the effect of a single pixel perturbation by Gaussian or folded Gaussian noise on the decision of the neural network.

In Section 7, we introduce L ∞ versions of SSAA and compare the L ∞ version of VFGA with the well-known projected gradient descent attack (PGD) baseline [START_REF] Madry | Towards deep learning models resistant to adversarial attacks[END_REF] on CIFAR-10. Results show the superiority of VFGA over PGD.

Finally, a nice property of SSAA is that they can be combined together in a single attack. While this combination is done only through VFGA (which combines two noises) in the present paper, we expect that the combination of more noises can significantly boost the performance of SSAA in both the L 0 and L ∞ version. This point will be discussed in the last Section 8.

Scalable stochastic sparse adversarial attacks (SSSAA)

In what follows, we introduce SSSAA by sampling from particular random noises on selected components of the input. To simplify, we mainly discuss targeted attacks and then deduce untargeted ones by applying slight modifications. The aim herein is to iteratively identify the best component to perturb and the best move for this component until the target label becomes the most probable for the NNC.

First, consider a Gaussian noise X θ ∼ N (0, θ) and denote by (e 1 , • • • , e n ) the basis of R n . Any c-targeted probability expectation of the perturbed input x + X θ e i can be expanded as follows:

E[F c (x + X θ e i )] = F c (x) + θ(G i F c )(x) + ... (1)
where

G i F c = 1 2 ∂ 2 F c ∂x 2 i
is the infinitesimal generator of X θ seen as a diffusion.

When taking the folded Gaussian noise X θ ∼ |N (0, θ)|, this expansion becomes:

E[F c (x + X θ e i )] = F c (x) + 2θ π ∂F c ∂x i (x) + θ(G i F c )(x) + ... (2)
For any given noise Y θ , we build our reasoning upon a heuristic which is to look for the input feature i that maximises

E[F c (x + Y θ e i )].
The assumption behind this heuristic is that searching for the best expectation will allow to discover the best moves according to the distribution of the noise Y θ . As a first remark, notice that the Gaussian noise X θ ∼ N (0, θ) is not a good candidate since it requires computing second derivatives of the NNC and moreover the approximations

E[F c (x + X θ e i )] ≈ F c (x) + θ(G i F c )(x) and V ar[F c (x + X θ e i )] ≈ θ ∂F c ∂x i 2 (
x) are of the same order as θ, indicating that variance should be taken into account in selecting the best components to perturb. For these reasons, we do not discuss the Gaussian noise in this section and postpone its study to Section 5.

On the other hand, considering the folded Gaussian noise, the approximation

E[F c (x + X θ e i )] ≈ F c (x) + 2θ π ∂F c ∂x i (x
) induces a negligible variance (only terms of θ p with p ≥ 1) in front of the expectation, at least when |θ| < 1 and allows us to only focus on the expected probability of the perturbed input. Note also that the approximation of the expected probability only contains first derivatives w.r.t. to the input component which is a practical advantage of the folded over the pure Gaussian noise. Taking a folded noise has another important advantage for bounded inputs. Please note that, without loss of generality, we consider inputs bounded in [0, 1] in this paper, as well as the adversarial samples which share the same support domain. In the following, we propose to automatically tune the variance parameter θ of X θ according to the distance of the input x i to these bounds. Please note that, for a given component i, x i ̸ = 0.5, the possible amplitude of move is not the same in both directions. Considering a Gaussian noise, since symmetric, would be problematic for this θ tuning. Rather, considering two folded Gaussian noises for each component, one positive (only for component increase) and one negative (only for component decrease) allows better fitted selections.

Next, we present a one-sided, only increasing perturbations, stochastic attack based on the folded Gaussian noise called FGA. Then, we introduce, in the same way another attack that relies on the uniform noise called UA. Finally, we deduce a simple both-sides attack, that considers the best choice between the increasing and decreasing FGA at each selected component, called Voting Folded Gaussian Attack (VFGA).

Folded Gaussian Attack (FGA)

For our one-side targeted attack FGA, the most relevant input feature to perturb is thus selected by the rule i = argmax j θ j ∂F c ∂x j , considering a folded Gaussian noise

|N (0, θ i )|.
Choosing the variance θ i . Since FGA only considers positive perturbations of the input, fixing the variance θ i must consider the upper-bound of the input domain. A quite natural choice could be either

θ i = 1 -x i (variance = 1 -x i ) or √ θ i = 1 -x i
(standard deviation = 1 -x i ). We choose √ θ i = 1 -x i to ensure that a generated perturbation x i +N i to x i has probability 2/3 to be inside the interval [x i , 1] (before clipping to [0, 1]) which is a more motivated choice. Our experimental results (not reported in this paper) show that this choice gives slightly more effective attacks than the second one.

After selecting the input feature i, our proposal is to simulate N S samples from |N (0, θ i )| to find an accurate move towards a close adversarial sample. The complete process is depicted in Algorithm 1 introducing the increasing FGA (and the decreasing FGA by analogy).

Algorithm 1 (Increasing) Folded Gaussian Attack (FGA) Inputs: x: input to the NNC of label l, c ̸ = l: targeted class, N S : number of samples to generate at each iteration, maxIter: maximum number of iterations. Output: x: adversarial sample to x.

Initilialise x ← x, Γ ← {1, • • • , dim(x)} \ {i : xi = 1}, iter ← 0. while Γ ̸ = ∅, label(x) ̸ = c and iter < maxIter do i 0 = argmax i∈Γ (1 -xi ) ∂F c ∂x i (x).
Generate samples S h 1≤h≤N S from |N (0, θ i0 )| where θ i0 := 1 -xi0 .

for h ∈ [[1, N S ]] do Define the input ỹ h by ỹ h j ← Clip [0,1] xj + S h if j = i 0 ỹ h j = xj otherwise. Batch compute F c ỹ h ; h ∈ [[1, N S ]] . x ← argmax ỹ h F c ỹ h , Γ ← Γ \ {i 0 } iter ← iter + 1.

return x

Choosing N S . The main hyperparameter in Algorithm 1 is the number N S . Given the way it is defined, one can expect that increasing N S will increase, up to saturation, the effectiveness of the attacks. This may, however, slow down their speeds. Thanks to batch computing, with sufficient memory, Step 8 can be performed at the cost of N S = 1 and (reasonably) augmenting N S can make Algorithm 1 converge faster as less iterations would be needed. In most of our experiments, we fix this number to N S = 10 but also address some comparisons with N S = 20, 100. We refer to the analysis of the experimental results for more discussions related to this point. Finally, we notice that batch computing used here does not often require a parallel computing effort by the user as this option is available in standard libraries.

The previous attack only applies perturbations that increase the input. Lowering the input features intensities can be as effective as increasing them. Following the same analogy, we introduce the decreasing FGA attack by taking

√ θ i = x i rather than √ θ i = 1-x i and replacing |N (0, θ i )| with -|N (0, θ i )|
in the previous algorithm. Note that FGA and XSMA are one sided attacks but XSMA apply predefined maximal perturbations and FGA explores in real time best perturbations to apply.

Uniform Attack (UA)

The UA can be introduced in the same way as FGA by sampling from the uniform noise. It follows the lines of Algorithm 1 (in its untargeted form), but Step 3 is

replaced with i 0 = argmin i∈Γ (1 -x i ) ∂F c ∂x i (x) and sampling in Step 4 is done from U([0, θ i ]), θ i = 1 -x i .

Voting attacks (VA)

In this section, we propose a two-sided attack, which both considers

E[F c (x + |X θ + i | e i )] and E[F c (x-|X θ - i | e i )]
for each feature, with X θ ∼ N (0, θ), θ + i = 1-x i and θ - i = x i . This method applies increasing and decreasing FGA at each iteration and chooses the most effective moves in both directions and is called voting FGA (VFGA) and is presented in Algorithm 2. Note that, in the same way, by considering the uniform noise it is possible to introduce more examples of voting attacks.

Untargeted attacks

The main focus for these attacks is to decrease the class probability of the input until a new class label is found. Only few modifications are required to deduce the untargeted versions of the previous Algorithms: by assuming c is the true label of x and replacing argmax with argmin in Steps 3 and 8 of Algorithm 1 and making similar slight changes in Algorithm 2.

Experiments on untargeted SSSAA

In this section, we present experiments to highlight the benefits of our untargeted SSSAA. First, we aim to showcase the relevance of FGA in comparison with UA. Second, we aim to compare our methods and more specifically VFGA with relevant state-of-the-art approaches. The codes of the present and next section are available in PyTorch at https://github.com/hhajri/stochastic-sparse-adv-attacks.

Algorithm 2 Voting Folded Gaussian Attack (VFGA) Inputs: as Algorithm 1. Output: as Algorithm 1.

Initilialise x ← x, Γ ← {1, • • • , dim(x)} \ {i : xi = 1}, iter ← 0. while Γ ̸ = ∅, label(x) ̸ = c and iter < maxIter do i + = argmax i∈Γ (1 -xi ) ∂F c ∂x i (x), i -= argmin i∈Γ xi ∂F c ∂x i (x). Generate samples S +, h 1≤h≤N S from |N (0, θ + i )| where θ + i := 1 -xi + . Generate samples S -, h 1≤h≤N S from -|N (0, θ - i )| where θ - i := xi -. for h ∈ [[1, N S ]] do Define the input ỹ h by    ỹ ±, h j ← Clip [0,1] xi ± + S ±, h if j = i ± ỹ ±, h j = xj otherwise. Batch compute F c ỹ ±, h ; h ∈ [[1, N S ]] (2N S propagations). x ← argmax ỹ ±, h F c ỹ ±, h , Γ ← Γ \ {i 0 } with i 0 = i + or i -according to the best move; iter ← iter + 1.

return x

In the experiments, we consider two popular computer vision datasets illustrating small and high dimensional data: CIFAR-10 14 (32 × 32× 3 images divided into 10 classes) and ImageNet 25 (ILSVRC2012 dataset containing 299 × 299× 3 images divided into 1,000 classes). The used neural network classifiers are described in the upcoming paragraphs.

The state-of-the-art attacks considered for comparison in this section are: SparseFool 20 . This method is fast and scalable. At each iteration, it applies Deep-Fool [START_REF] Moosavi-Dezfooli | Deepfool: a simple and accurate method to fool deep neural networks[END_REF] to estimate the minimal adversarial perturbation thanks to a linearization of a classifier. Then, it estimates the boundary point and the normal vector of the decision boundary and finally updates the input features with a linear solver.

Brendel and Bethge L 0 attack (B&B) [START_REF] Brendel | Accurate, reliable and fast robustness evaluation[END_REF] . This adversarial attack is gradientbased. It follows the boundary between the space of adversarial and non-adversarial images to find the minimum distance to the clean image. It is powerful and more efficient (but also slower and more complex) than many gradient-based approaches such as SparseFool.

GreedyFool 10 . This attack is an improvement of SparseFool. It is however more complex than the later as it needs to carefully train a distortion map which is a generative adversarial network GAN [START_REF] Goodfellow | Generative adversarial networks[END_REF] . We remark (based on one experiment on ImageNet) that it is less efficient (but also faster and less complex) than B&B.

Sparse-RS [START_REF] Croce | Sparse-rs: a versatile framework for query-efficient sparse blackbox adversarial attacks[END_REF] . This attack is fast and achieves high success rate on ImageNet outperforming many white-box attacks such as PGD 0 [START_REF] Croce | Sparse and imperceivable adversarial attacks[END_REF] . It requires fixing the maximum number of pixels to modify which is then fully exploited. In order to generate adversarial examples with minimal L 0 perturbations by Sparse-RS, one needs to run this method for several budgets before selecting an optimal budget.

A notable difference with Sparse-RS. It should be mentioned that our attacks and Sparse-RS follow different strategies. Indeed, the budget k for Sparse-RS is fixed in the pixel space. For instance, on CIFAR-10 this can go up to 32 × 32, and once k is fixed the number of modified pixels in the input space, for Sparse-RS, is near 3 × k. Our attacks compute perturbations directly in the input space. All attacks are however L 0 in the usual definition and they are compared according to the most commonly used metric which is, up to our knowledge, the L 0 distance in the input space.

All the previous attacks are experimented using the original implementations by the authors and following the recommended hyperparameters.

To compare between the different methods, we rely on the following scores: success rate (SR), Mean and Median number of changed pixels, complexity based on the number of model propagation 9 (MP). We prefer MP over the running time per image since it is not dependent on the software used when executing the codes.

More approaches. In this paper, we propose fast methods, and thus we only focus on comparisons with similar fast approaches like SparseFool and Sparse-RS. The B&B, although not fast, has been selected as a highly efficient benchmark attack. We omit comparison with Carlini&Wagner L 0 [START_REF] Carlini | Towards evaluating the robustness of neural networks[END_REF] and we believe the results would be similar to our comparison with B&B. Also, we omit comparison with CornerSearch 8 because it is less effective than Sparse-RS based on the work [START_REF] Croce | Sparse-rs: a versatile framework for query-efficient sparse blackbox adversarial attacks[END_REF] and also needs a large computational cost on ImageNet (see results in 8 ).

On CIFAR-10.

On this dataset, we use the ResNet18 13 and VGG-19 26 models. After training with PyTorch, these networks reached 95.55% and 93.87% accuracies respectively. For our attacks UA, FGA and VFGA, the hyperparameter N S is fixed to N S = 10 (the obtained attacks are denoted UA10, FGA10 and VFGA10). The effect of augmenting N S is analysed later on in this section. We notice that, otherwise stated, maxIter is put to its maximal value. The state-of-the-art approaches outlined before, except GreedyFool, are tested and compared with our methods on the correctly predicted samples among the 10, 000 CIFAR-10 test images. We refrained from comparing with GreedyFool because of the need to train the distortion map network on CIFAR-10 not provided in the code of 10 (this network has been made available for the ImageNet dataset and comparison on this dataset is considered in the next section). For Sparse-RS, several budgets of pixels k (the number of pixels to modify) have been experimented and the smallest budget giving 100% has been selected. On CIFAR-10, we choose the optimal k, i.e k ′ = k -1 does not give full success of the attack. Our intention is to show that under the condition of full success for all attacks (when possible) our VFGA method is overall more advantageous. Comments. The previous results show that the folded Gaussian noise is more advantageous in attacking than the uniform noise and that combining two folded distributions is useful not only for the SR, Mean and Median but also for the model propagation score. Concerning the comparison with the state-of-the-art methods: VFGA has less advantageous Mean and Median than B&B (near 2 times greater for Mean and the gap is reduced for Median). Nevertheless, it is up to 1 20 less complex based on the MP score. Second, all of our methods and more particularly VFGA significantly outperform SparseFool. Regarding the comparison with Sparse-RS, we remark that VFGA has notable Mean and Median advantages (up to 3 times fewer) and is also less complex given the number of experiments carried for Sparse-RS to achieve full success (with minimal Mean and Median). Notice also the difficulty to find a good k with full success for Sparse-RS as the optimal value depends on each sample and high values impact the overall performance of this attack. An advantage for us is that this parameter is set automatically and is optimal for each sample.

Attacks

A comparison between Sparse-RS and VFGA for different distortions.

As stressed before, we only focus on performances under the condition of full success which is usually reported to summarise the contribution of new methods. If we relax this condition, we remark that for small budgets k when both VFGA and Sparse-RS are not fully successful, Sparse-RS outperforms VFGA in SR but VFGA obtains better Mean and Median which are always near k for Sparse-RS. Starting from a k which approaches full success, VFGA becomes more advantageous in SR, Mean and Median. We draw in Figure 2 the SR (on the left), and Mean and Median L 0 scores (on the right) versus the number of perturbed pixels for VFGA10 and Sparse-RS. When only few pixels are perturbed, Sparse-RS performs better than VFGA10 in SR. Once the number of perturbed pixels exceeds 105, the two SR become very competitive (near and then equal to 100%). The Mean and Median L 0 scores for VFGA10 are, however, always better than those by Sparse-RS regardless the num-ber of pixels which have been modified.

For Sparse-RS, Median is a linear function of the number of perturbed pixels. For the budget of 285 modified pixels, Mean is 282.78. For VFGA10, this budget is the worst L 0 score for the 10,000 test images of CIFAR-10 and with this value, Mean is 17.03 which is about 16 times fewer than Sparse-RS. In short, even if for a small number of disturbed pixels, VFGA10 is not able to reach 100% in SR, its Mean and Median are still better than those of Sparse-RS.

Augmenting N S .
In what follows, we investigate the impact of augmenting N S on the performances of our attacks by testing UA20, FGA20 and VFGA20, which correspond to N S = 20, on the same data as Table 1. We remark that SR, Mean and Median are slightly improved but based on the MP score, the attacks become more complex. This illustrates the fact that increasing so much N S may not significantly improve the attacks but on the other hand it may slow down them. Also, we observe that despite augmenting N S from 10 to 20, the uniform attack cannot beat FGA10. This is quite remarkable since for the uniform distribution the N S generated samples are different and fall inside the domain of the input features while for the folded Gaussian distribution, due to clipping, several samples are likely to be clipped at the minimal and maximal bounds. Augmenting N S also increases more quickly MP for the uniform noise.

Attacks

On ImageNet.

In this section, we test the ability of the previously tested attacks and additionally GreedyFool to generate adversarial examples at large scale by considering models on ImageNet. Two pre-trained networks provided by PyTorch are considered for testing: Inception-v3 [START_REF] Szegedy | Rethinking the inception architecture for computer vision[END_REF] and VGG-16 [START_REF] Simonyan | Very deep convolutional networks for largescale image recognition[END_REF] whose accuracies are respectively 77.45% and 71.59%. Inputs are of size 299×299×3 and 224×224×3 for the first and second model respectively.

Again, we consider B&B as a benchmark of a highly successful attack. Greedy-Fool requires training a GAN network on ImageNet but once carefully done it is highly successful. We recover the GAN model from the available code of 10 which we complement by adding the computation of the MP score for this attack. For Sparse-RS, we again fix our objective to compare with this attack when 100% SR is achieved. This requires launching several experiments for different values of k on the whole considered set of images in order to obtain a near-optimal budget value. By this, we mean a value k giving 100% SR ; there exists k ′ < k and the performances of VFGA in Mean and Median are better than those obtained by Sparse-RS with budget k ′ . This implies in particular that VFGA gives better results than Sparse-RS when tested with the optimal value of k. To get an idea of the difference between our results and those by Sparse-RS, we always report the results for k ′ and k by Sparse-RS (in this section and next one).

The obtained results for the different attacks are reported in Table 3 and commented after. Comments. First, B&B achieves the best SR, Mean and Median scores. Greedy-Fool comes after but with the cost of training a GAN model on ImageNet. The complexity comparison between these two attacks is difficult to address and we only claim that both methods are significantly more complex than our approach (GreedyFool is complex to reproduce on new datasets). Despite this fact, we observe that on VGG-16 VFGA has a gap of Mean and Median of less than 12 pixels which is relatively small. Among the methods shown in the previous table, our attacks and SparseFool are the fastest under the full success condition and when minimising at the same time Mean and Median. Our attacks, obtain, however, overall better performances than SparseFool according to all metrics. Specifically, VFGA always significantly outperforms SparseFool. Moreover, despite the fact that we select near-optimal values of k for Sparse-RS, VFGA is still more advantageous regarding Mean and Median and also faster if the complexity of finding k is added. Finally, we notice that after finding a good k giving full success for Sparse-RS, this attack can not generate relevant adversarial examples with minimal L 0 distance 1, while due to the flexibility of our attack, several such examples can be generated. This is a further advantage of our attack.

Experiments on targeted SSSAA

Targeted attacks are more challenging than untargeted ones. The objective of this section is to compare (targeted) VFGA, our selected method, with (targeted) Sparse-RS as a fast attack outperforming several state-of-the-art methods [START_REF] Croce | Sparse-rs: a versatile framework for query-efficient sparse blackbox adversarial attacks[END_REF] recall that SparseFool is not efficient as a targeted attack. We do not report results by FGA and UA but claim that FGA is still more relevant than UA and only omit to address a similar comparison as before. We do not report results by B&B and GreedyFool as targeted attacks because of the need of the distortion map on CIFAR-10 for GreedyFool, the non ability to reproduce B&B in the targeted mode and moreover since, we consider that these approaches are complex to reproduce on new datasets. Thus, we only focus on the comparison with Sparse-RS and defend our approach as an efficient fast method. Our main conclusion in this paragraph is that, for ImageNet which is more challenging, VFGA is still more relevant than Sparse-RS regarding the same previous metrics when both attacks are fully successful and despite the fact that a near-optimal value of k is selected for Sparse-RS. On CIFAR-10, we conclude on the basis of our experiments that Sparse-RS is more advantageous in Mean and Median.

We consider the same network models as before. To simplify the experimentation, we do not consider all possible target labels but, for each test dataset, we generate a list of random labels which were fixed once for all. Each input image is then attacked to have one desired label. For Sparse-RS, we again select a nearoptimal k in all experiments. This task took several hours on CIFAR-10 and several days on ImageNet. Tables 4 and5 show the results obtained on CIFAR-10 and ImageNet (in 4: Results on the correctly-predicted test images of CIFAR-10. GS is a greedsearch to find near-optimal values of k that took several hours. The results for k ′ are not reported since Sparse-RS is here more advantageous in full success.

First, we notice that Sparse-RS obtains better Mean and Median than VFGA10. Increasing N S from 10 to 100 improves considerably VFGA but our results are still less better than Sparse-RS. Given the time needed to find the near-optimal values k, we claim that our attacks are still overall much faster than Sparse-RS to obtain full success with optimal Mean and Median.

Our interpretation of Table 5 is overall similar to that of Table 3. When attacking the Inception-v3 model, which is more challenging, VFGA10 outperforms Sparse-RS regarding all scores. When attacking VGG-16, Sparse-RS only takes a slight advantage of Mean. As for untargeted attacks, a notable advantage of our methods is the flexibility of the number of modifiable pixels allowing us to generate adversarial examples with minimum L 0 distance while being 100% successful on all samples. Table 5: Results obtained on the 5,000 firstly correctly-predicted validation images of ImageNet. GS is a greed-search to find near-optimal values of k that took several days.

Attacks

Less scalable SSAA

In this section, we focus on sparse targeted and untargeted perturbations by a distribution of the form N (aθ, b 2 θ) = aθ + bN (0, θ) where θ is a parameter to be tuned as before using the bounds of the input domain and a, b are scaling parameters which determine the contribution of the deterministic component θ and the pure Gaussian one N (0, θ). The resulting attack is denoted GA(a,b). Note that the multidimensional Gaussian noise has been used to generate effective adversarial L 2 and L ∞ attacks in [START_REF] Li | Nattack: Learning the distributions of adversarial examples for an improved black-box attack on deep neural networks[END_REF] . Our attack can be thought of as a sparse version of 18 since it follows a similar methodology based on the minimisation of the expected probability.

If X a,b θ ∼ N (aθ, b 2 θ), then E[F c (x + X a,b θ e i )] = F c (x) + θ(A a,b F c ) c,i (x) + o(θ) (3)
where

(A a,b F c ) c,i (x) := a ∂pc ∂xi (x) + b 2 2 ∂ 2 Fc ∂ 2 xi (x)
. Before explaining the minor changes to introduce G(a,b), let us remark that our GA requires computing the diagonal of the Hessian of F c with respect to the input coordinates which is costly on large datasets. For this reason, we make the following assumption under which these calculations become more feasible (at least on small datasets):

All hidden layers activation functions of the NNC are ReLU (4)

Note that state-of-the-art NNC make use of ReLU and so (4) is not a limitation as such. Under (4), the Hessian diagonal of the NNC is a simple function of the NNC Jacobian as shown here.

Fast Hessians diagonals computations for ReLU activations. Call L i the logits of the NNC. Since each L t is the composition of affine and ReLU functions which both have null second derivatives, by a simple recurrence argument ∂ 2 L t (x) ∂x i = 0. Using this fact and deriving twice the identity F t (x) × K j=1 e Gj (x) = e Gt(x) , we deduce the main relations:

• ∂Ft ∂xi = F t ( ∂Gt ∂xi -j F j ∂Gj ∂xi ) • ∂ 2 Ft ∂ 2 xi = ∂Ft ∂xi ( ∂Gt ∂xi -j F j ∂Gj ∂xi ) -F t × j ∂Fj ∂xi ∂Gj ∂xi
Experimentally, this implementation gives a very accurate approximation of the Hessian diagonal calculated with PyTorch or Tensorflow.

The (targeted) attack G(a,b) follows Algorithm 1; where Step 3 should be re-

placed with i 0 = argmax i∈Γ A a,b F c (x). For Step 4, sampling is done from N (aθ, b 2 θ i0 )
where θ i0 is to be set by the user. Next, we only experiment G(0,1) (the purely Gaussian case) and G(1,1) (the mixed balanced deterministic-Gaussian case). The parameter θ i0 for these attacks is fixed as follows:

• θ i0 = max(x i0 , 1 -xi0 ) for G(0,1). • θ i0 = 1-xi 0 2
for G(1,1).

As a complement, note that the choice (a, b) = (1, 0) is the one which has led to the XSMA attacks.

In the rest of this section, we compare the untargeted GA(0,1) and GA(1,1) with UA and FGA in attacking LeNet-5 model on MNIST [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF] . We use the same model proposed in 5 which has 99.98% accuracy on the training dataset and 99.49% accuracy on the test set. We report in Table 7 The results of this table can be ranked as follows (a<b means b outperforms a) GA(0,1)<UA<GA(1,1)<FGA. First, they confirm the overall superiority of FGA over UA yet on the MNIST dataset. GA(0,1) has overall less better results than the other attacks which is somewhat expected since it has been remarked in [START_REF] Papernot | The limitations of deep learning in adversarial settings[END_REF] that, on MNIST, increasing the input features works better than decreasing them and GA has less chance to increase the input features than the other attacks.

While the previous results could be explained by the nature of the MNIST dataset, since increasing the pixel values is known to give better results [START_REF] Papernot | The limitations of deep learning in adversarial settings[END_REF] , we do not expect similar behaviours on other datasets where there is no privileged direction for the attacker. Our main claim is that, independently of each attack performance, combining all of them (in the same way as VFGA) and possibly varying the intensity θ can lead to more efficient L 0 attacks. Indeed, experimentally, we observe many successful adversarial examples which modify a single or very few pixels by each method of Table 7 whereas, at the same time, the other methods fail in doing so. Two illustrations are given through Figure 3.

Fig. 3: Top: left-right: an input of label 3 and its perturbation by GA(0,1) (interpreted as 5) by modifying one pixel. The L 0 scores for FGA, GA(1,1) and UA are respectively: 2, 5, 2. Down: left-right: an input of label 9 and its perturbation by GA(1,1) (interpreted as 7) while modifying one pixel. The L 0 scores for GA(0,1), FGA and UA are respectively: 20, 6, 8.

A fine analysis of one-pixel perturbations by the Gaussian and folded Gaussian noise for ReLU-networks

In this section, we want to interpret the behaviour of a ReLU NNC under one pixel perturbation by a Gaussian or Folded Gaussian distribution. We focus on the output logits (G i ) i as they fully determine the decision of the NNC.

The main result is the next Theorem which gives a decomposition of the expected logit E[G c (x + X θ e i )], where X θ ∼ N (0, θ) or X θ ∼ |N (0, θ)|. Theorem 6.1. Denote by (G i ) 1≤i≤p the output logits of a ReLU-neural network classifier. Then for every c ∈ [1, p], we have

(i) E[G c (x + N (0, θ)e i )] = G c (x) + 1 2 j∈D Γ j θ (Gr j+ -Gr j-) (ii) E[G c (x + |N (0, θ)|e i )] = G c (x) + Γ θ ∂Gc ∂xi (x + 0 + e i ) + 1 2 j∈D + Γ j,+ θ (Gr j+ - Gr j-)
where Gr u = ∂ u G c (x + ue i ) and D is the set of discontinuities (or jumps) of Gr,

Γ j θ = E[|N (0, θ) -j|] -|j|; D + is the set of jumps of Gr u , u > 0 in [0, θ], Γ θ = 2θ π and Γ j,+ θ = E[||N (0, θ)| -|j||] -|j|, ∂Gc ∂xi (x + 0 + e i ) = lim h→0+ ∂Gc ∂xi (x + he i ).
Proof. Let us give the main steps to prove (i). By applying Itô-Tanaka formula (see [START_REF] Revuz | Continuous martingales and brownian motion[END_REF] ) for the function g(u) = G c (x + ue i ) (which can be written as the difference of two convex functions and so satisfies the assumption of the formula), we have

G c (x + W θ e i ) = G c (x) + M θ + 1 2 L a θ T ′′ g (da)
where M θ is a martingale; in particular E[M θ ] = 0; L a θ is the local time of W at level a and time θ and T ′′ g is the second derivative in the sense of distributions of g. Next, we use the decomposition of distributions:

(T g ) ′′ = T g ′′ + a∈D ′ g (g ′ (a+) -g ′ (a-))δ a
where D ′ g is the set of discontinuities of g ′ . Note that T g ′′ = 0 since g is a piecewise constant function. Thus, with the notation of the theorem, we have

E[G c (x + W θ e i )] = G c (x) + 1 2 j∈D Γ j θ (Gr j+ -Gr j-)
where Γ j θ = E[L j θ ] which also coincides with the expression given in the theorem. (ii) can be proved in a similar way.

Let us give some remarks regarding the previous theorem:

• In the particular case of one hidden layer, Theorem 6.1 (i) allows to retrieve the expression found in Theorem 1 in 2 . In this case, the set D is reduced to one point and there is only one jump of Gr. The expression given here is much more interpretable than 2 and general for any ReLu-based network. Note also that the number of points contained in D is bounded by the number of hidden layers.

• The constants Γ j θ , Γ θ , Γ j,+ θ in Theorem 6.1 are easy and fast to approximate with high precision. Thus the main complexity in the expansions given above is the computation of the jumps of Gr. Once these jumps are computed efficiently, they can be used once for all to approximate all E[G c (x + X θ e i )], X θ ∼ N (0, θ) or |N (0, θ)| for every possible θ. In short, we claim that Theorem 1 reduces the computation of E[G c (x + X θ e i )] for all θ to the computation of the jumps of the trajectory of the gradient of G c (x) with respect to the i-th coordinate.

In some situations, when simple Gaussian sparse attacks or also our attacks such as GA(0,1) and FGA are successful in one iteration, the user or designer of the NNC would be interested in understanding such a failure. The main tool given in this section suggests to compare the quantities given in front of the expectations in Theorem 6.1 for different c, which as seen before, depends on the jumps of the gradient paths of the logits according to the chosen component.

L ∞ bounded attacks

In this section, we focus on L ∞ attacks. This category seems to be more studied in the literature than L 0 and L 2 ones. One reason is that L ∞ attacks have reached the-state-of-the-art performances in Adversarial Training [START_REF] Madry | Towards deep learning models resistant to adversarial attacks[END_REF] . Indeed, [START_REF] Madry | Towards deep learning models resistant to adversarial attacks[END_REF] introduced a strong L ∞ baseline, the PGD attack, to evaluate NNC and adversarially train them. Several improvements of PGD have been proposed but they all either follow the main idea of PGD or combine many PGD varieties at once (see [START_REF] Croce | Robustbench: a standardized adversarial robustness benchmark[END_REF] ).

Our attacks presented previously can be applied as L ∞ ones by slight modifications. Algorithm 1 below introduces the L ∞ ε-bounded version of VFGA where the main idea is to sample perturbations from N (0, ε). Notice that N S is the only hyperparameter for the attack. Its interpretation is almost the same as before. Our main claim is that, despite the iterative scheme of VFGA which makes the attack slow in comparison with PGD, this attack can outperform PGD. We propose an experiment to support this by attacking the previous CIFAR-10 Resnet model on 1000 correctly predicted input test images. The applied attacks to this model are VFGA, PGD and PGD with random restart denoted PGD-RR [START_REF] Madry | Towards deep learning models resistant to adversarial attacks[END_REF] for two values of ε = 4/255, 2/255. The results are given in Table 7.

We remark that for ε = 4/255, the performances are very competitive between our attack and the two PGD varieties. However, for for ε = 4/255 the results are significantly better for VFGA.

Conclusion

This paper introduced noise-based attacks to generate sparse adversarial samples to inputs of deep neural network classifiers. A first advantage of our methods is that they work as both untargeted and targeted attacks. Moreover, they are very simple to put in place and require fixing only one parameter whose interpretation x is intuitive (the bigger the best up to saturation in performance). Our attacks are faster to apply on new models and datasets than existing approaches (SparseFool, GreedyFool) while ensuring full success. They are much less complex than thestate-of-the-art method B&B relying on the model propagation score (near 1 20 on CIFAR-10 and 1 5 on ImageNet) and achieve competitive results in some cases. In comparison with Sparse-RS, our attacks are flexible allowing to find an optimal budget of pixels for each input image and achieve full success with minimal L 0 scores. They are applicable as L ∞ attacks and can outperform the PGD baselines. We believe that our results can be improved in a crucial way by combining different kind of noises as explained in the paper.

Fig. 2 :

 2 Fig. 2: Results over the number of perturbed pixels for untargeted attacks (VFGA10 and Sparse-RS) on CIFAR-10 for the ResNet-18 model. On the left, the success rate SR and on the right, the Mean and Median L 0 scores.

Algorithm 3 L←

 3 ∞ VFGA with intensity ε Inputs: as Algorithm 1. Output: as Algorithm 1.Initilialise x ← x, Γ ← {1, • • • , dim(x)} \ {i : xi = 1}, iter ← 0. while Γ ̸ = ∅, label(x) ̸ = c and iter < maxIter do i + = argmax i∈Γ ∂F c ∂x i (x), i -= argmin i∈Γ ∂F c ∂x i (x). Generate samples S +, h 1≤h≤N S from Clip [0,ε] (|N (0, ε)|); Generate samples S -, h 1≤h≤N S from -Clip [0,ε] (|N (0, ε)|); for h ∈ [[1, N S ]] do Define the input ỹ h by Clip [0,1] xi ± + S ±, h if j = i ± ỹ ±, h j = xj otherwise. Batch compute F c ỹ ±, h ; h ∈ [[1, N S ]] (2N S propagations).

Table 1 :

 1 Results on the correctly predicted samples among the 10, 000 test images of CIFAR-10. SR is the success rate of the attack, Mean, Median are the average and median number of modified pixels on successful samples and MP is the number of model propagations. GS is a greed-search to find optimal values of k giving full success that took several hours. The highlighted results of our VFGA in comparison with the fast methods SparseFool and Sparse-RS are in bold.

		SR	Mean Median	MP
			ResNet18
	B&B	100	8.33	8.0	1927
	SparseFool	99.31	36.48	9.0	520
	Sparse-RS (k = 10)	100	29.79	30	GS + 49
	UA10	100	30.94	20.0	363
	FGA10	100	29.70	20.0	134
	VFGA10	100	17.03	11.0	99
			VGG-19
	B&B	100	5.30	6.0	1483
	SparseFool	97.98	67.71	8.0	686
	Sparse-RS (k = 7)	100	20.82	21	GS + 55
	UA10	100	22.16	11.0	281
	FGA10	100	19.67	11.0	103
	VFGA10	100	11.40	7.0	80

Table 2 :

 2 Comparison between our attacks for N S = 20.

		SR	Mean	Median	MP
			ResNet18
	UA20	100	30.93	20.0	684
	FGA20	100	30.23	19.0	193
	VFGA20	100	16.76	11.0	131
			VGG-19
	UA20	100	20.60	10.0	457
	FGA20	100	19.71	10.0	134
	VFGA20	100	11.22	7.0	113

Table 3 :

 3 Results on the firstly 6, 000 correctly predicted validation images of Ima-geNet. GS is a greed-search to find near-optimal values of k that took several days. GAN is a generative network trained on ImageNet. Our results in comparison with the fast methods SparseFool and Sparse-RS when fully successful are highlighted in bold.

	. We

Table 4

 4 VFGA100 is simply VFGA with N S = 100).

	Attacks	SR	Mean	Median	MP
			ResNet18
	Sparse-RS (k = 30)	100 89.43	90.0	GS + 1678
	VFGA10	100 641.49	174.0	13427
	VFGA100	100 154.43	105.0	20619
			VGG-19
	Sparse-RS (k = 25)	100 73.17	75.0	GS + 1123
	VFGA10	100 551.27	150.0	12213
	VFGA100	100 174.93	97.0	21417

Table

Table 6 :

 6 the results on 1, 000 images. Results obtained on 1, 000 images of the MNIST test set.

	Attacks	SR	Mean	Median
	UA	100%	19.53	13.0
	FGA	100%	17.69	13.0
	GA(0,1)	100%	22.45	16.0
	GA(1,1)	100%	17.68	12.0

Table 7 :

 7 Results obtained on the first 1000 correctly classified images of the CIFAR-10 test set.

	return	x
		Attacks	ε	SR
		PGD	2/255	63.5%
		PGD-RR	2/255	67.5%
		VFGA	2/255	86.4%
		PGD	4/255	97.6%
		PGD-RR	4/255	98.4%
		VFGA	4/255	97.9%

← argmax ỹ ±, h F c ỹ ±, h , Γ ← Γ \ {i 0 } with i 0 = i + or i -according

to the best move; iter ← iter + 1.

Acknowledgements. This research work has been carried out in the framework of IRT SystemX, Paris-Saclay, France, and therefore granted with public funds within the scope of the French Program Investissements d'Avenir. This work is a part of the project EPI project (EPI for "Evaluation des Performances de l'Intelligence artificielle"). The project is supervised by IRT systemX and its partners which are Apsys, Expleo France, Naval Group and Stellantis. Also, we would like to thank Jérôme Rony for useful discussions and for his short implementation of our attack in 24 .