Guillaume Carlier
email: carlier@ceremade.dauphine.fr

Maxime Laborde
email: laborde@ceremade.dauphine.fr.

A splitting method for nonlinear diffusions with nonlocal, nonpotential drifts

Keywords: Wasserstein gradient flows, Jordan-Kinderlehrer-Otto scheme, splitting, nonlocal drift, nonlinear diffusions, Helmholtz decomposition. MS Classification: 35K15, 35K40

We prove an existence result for nonlinear diffusion equations in the presence of a nonlocal densitydependent drift which is not necessarily potential. The proof is constructive and based on the Helmholtz decomposition of the drift and a splitting scheme. The splitting scheme combines transport steps by the divergence-free part of the drift and semi-implicit minimization steps à la Jordan-Kinderlherer Otto to deal with the potential part.

Introduction

In their seminal paper [START_REF] Jordan | The Variational Formulation of the Fokker-Plank Equation[END_REF], Jordan-Kinderlehrer and Otto, showed that the Fokker-Planck equation

∂ t ρ -∆ρ -div(ρ∇V ) = 0, ρ |t=0 = ρ 0
with a probability density ρ 0 as Cauchy datum can be viewed as the gradient flow for the Wasserstein metric of the relative entropy with respect to the equilibrium measure ρ ∞ := e -V , E(ρ|ρ ∞ ) := ´ρ log(ρ/ρ ∞ ). They also introduced an implicit Euler-scheme (nowadays referred to as the JKO scheme): given a time step h > 0 and starting from ρ 0 h = ρ 0 , construct inductively a sequence ρ k h by:

ρ k+1 h ∈ argmin 1 2h W 2 2 (ρ, ρ k h ) + E(ρ|ρ ∞ )
where W 2 denotes the 2-Wasserstein distance (see section 2) and proved its convergence to the solution of the Fokker-Planck equation as h → 0. The theory of Wasserstein gradient flows has developed rapidly in the last twenty years with many applications for instance to porous medium equations [START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF] or aggregation equations [START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF].

The reference textbook of Ambrosio, Gigli and Savaré [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF] gives a very detailed account of this powerful theory, which enables to study general nonlinear diffusion equations of the form ∂ t ρ -∆P (ρ) -div(ρ∇V ) = 0 where P (ρ) := ρF ′ (ρ) -F (ρ),

as Wasserstein gradient flows for the energy ´(F (ρ) + V ρ).

The purpose of the present paper is to present a splitting transport-JKO scheme to study nonlinear diffusion equations (or more generally, systems) with a general density-dependent drift: ∂ t ρ(t, x) -∆P (ρ(t, x)) -div(ρ(t, x)U [ρ(t, .)](x)) = 0, t ≥ 0, x ∈ Ω, (∇P (ρ) + U [ρ]) • ν = 0 on ∂Ω, ρ |t=0 = ρ 0 , where for every probability density ρ, U [ρ] is a-not necessarily potential-vector field, for instance, in even dimensions, it can mix a gradient and Hamiltonian structures i.e. be of the form ∇V [ρ] + J∇H[ρ] (where J is the usual symplectic matrix). The potential case where U [ρ] = ∇V [ρ] can be studied by means of a semiimplicit JKO scheme introduced by Di Francesco and Fagioli [START_REF] Francesco | Measure solutions for nonlocal interaction PDEs with two species[END_REF] in the nondiffusive case and further developed by Laborde [START_REF] Laborde | On some non linear evolution systems which are perturbations of Wasserstein gradient flows[END_REF] for the case of a non linear diffusion. The idea of our splitting scheme is natural and consists in performing a Helmholtz decomposition of U [ρ]. We then treat the divergence-free part purely by (continuous in time) transport and the potential part by the semi-implicit JKO scheme. For the transport steps of the splitting scheme, we essentially need the divergence-free part to have some Sobolev regularity in x so as to be able to apply DiPerna-Lions theory, we will also need both the potential and divergence-free part of ρ → U [ρ] to satisfy some Lipschitz continuity condition with respect to the Wasserstein distance. In our recent work [START_REF] Carlier | Remarks on continuity equations with nonlinear diffusion and nonlocal drifts[END_REF], we studied the same type of equations or systems (in the periodic in space setting) by quite different arguments (approximation by uniformly parabolic equations). One advantage of the constructive splitting method presented here is that the transport steps by a divergence-free vector field preserve the internal energy, this is one way to overcome some difficulties discussed in [START_REF] Mészáros | Advection-diffusion equations with density constraints[END_REF] (section 5, variant 3).

The paper is organized as follows. Section 2 recalls some results from optimal transport and DiPerna Lions theory. Section 3 lists the various assumptions, explains the splitting scheme and gives the main result. Section 4 gives estimates on the discrete sequences of measures obtained by the splitting scheme. Convergence of the scheme as the time step goes to 0 to a solution of the PDE is proved in section 5. In the concluding section 6, we briefly discuss extension to systems and uniqueness issues.

Preliminaries

Wasserstein space

We recall some results from optimal transport theory that we will use in the sequel, we refer the reader to the textbooks of Villani [START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Villani | Optimal transport: Old and New[END_REF], Ambrosio, Gigli and Savaré [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF] or Santambrogio [START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF] for a detailed exposition. Let Ω be an open subset of R n , we denote by P(Ω) the set of Borel probability measures on Ω, P 2 (Ω) := {ρ ∈ P(Ω) : M (ρ) := ´Ω |x| 2 dρ(x) < +∞} and P ac 2 (Ω) the elements of P 2 (Ω) which in addition are absolutely continous with respect to the n-dimensional Lebesgue measure. Given ρ 0 and ρ 1 in P 2 (R n ) and denoting by Π(ρ 0 , ρ 1 ) the set of transport plans between ρ 0 and ρ 1 i.e. the set of Borel probability measures on R n × R n having ρ 0 and ρ 1 as marginals, the 2-Wasserstein distance between ρ 0 and ρ 1 , W 2 (ρ 0 , ρ 1 ) is defined as the value of the optimal transport problem

W 2 2 (ρ 0 , ρ 1 ) = inf γ∈Π(ρ0,ρ1) ˆRn × R n |x -y| 2 dγ(x, y).
Since this is a linear programming problem, it admits a dual formulation, which reads

1 2 W 2 2 (ρ 0 , ρ 1 ) = sup ˆRn ϕdρ 0 + ˆRn ψdρ 1 : ϕ(x) + ψ(y) ≤ 1 2 |x -y| 2 .
Optimal potentials in the problem above are called Kantorovich potentials between ρ 0 and ρ 1 , they can be taken semi-concave and their existence is well-known (see [START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Villani | Optimal transport: Old and New[END_REF][START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF][START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF]). If ρ 0 ∈ P ac 2 (Ω) a celebrated result of Brenier [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF] states that there is a unique optimal plan γ between ρ 0 and ρ 1 and it is induced by a transport map T i.e. is of the form γ = (id, T ) # ρ 0 and optimality is characterized by the fact that T = ∇u with u convex, moreover T (x) = x -∇ϕ(x) where ϕ is the Kantorovich potential between ρ 0 and ρ 1 . This in particular gives

W 2 2 (ρ 0 , ρ 1 ) = ˆRn |∇ϕ(x)| 2 ρ 0 (x)dx.
Another important result is the Benamou-Brenier formula [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF] which gives a dynamic formulation of W 2 2 (ρ 0 , ρ 1 ) and expresses it as the infimum of the kinetic energy:

ˆ1 0 ˆRn |v t (x)| 2 dρ t (x)dt
among solutions of the continuity equation

∂ t ρ + div(ρv) = 0, ρ |t=0 = ρ 0 , ρ |t=1 = ρ 1 .

Flows of weakly differentiable vector fields

We will also need to apply the DiPerna Lions theory [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF] in the special case of divergence-free vector fields. Let W ∈ W 1,1 loc (R n , R n ) be divergence-free div(W ) = 0 and with at most linear growth

|W (x)| C(1 + |x|).
Then there exists a unique flow map X : R

+ × R n → R n , X ∈ C(R + , L 1 loc (R n )) such that
• X(0, .) = id and for a.e. x, t ∈ R + → X(., x) is a solution of the ODE Ẋ = W (X) i.e.:

X(t, x) = x + ˆt 0 W (X(s, x))ds, t ≥ 0,
• X satisfies the group property X(t, X(s, x)) = X(t + s, x) for a.e. x and every t, s ≥ 0,

• for every X(t, .) preserves the n-dimensional Lebesgue measure.

Moreover given ρ 0 ∈ P ac 2 (R n ), ρ(t, .) := X(t, .) # ρ 0 = ρ 0 (X(t, .) -1 ) is the unique weak solution of the continuity equation

∂ t ρ + div(ρW ) = 0, ρ |t=0 = ρ 0 , (2.1)
which, since W is divergence-free, can also be rewritten as the transport equation

∂ t ρ + ∇ρ • W = 0. If we are
given an open subset Ω of R n with a smooth boundary and which is tangential (in the sense of traces) to ∂Ω, then the DiPerna-Lions flow X leaves Ω invariant so that if ρ 0 ∈ P ac 2 (Ω) (extended by 0 outside Ω, say), the solution ρ t = X(t, .) # ρ 0 of (2.1) remains supported in Ω hence may be viewed as a curve with values in P ac 2 (Ω).

Assumptions and main result

Given a suitable convex nonlinearity F and its associated pressure P (ρ) = ρF ′ (ρ) -F (ρ) as well as a nonlocal drift ρ → U [ρ], our goal is to solve

∂ t ρ -∆P (ρ) -div(ρU [ρ]) = 0, ρ |t=0 = ρ 0 , (3.1)
on (0, +∞) × Ω, where Ω is a smooth domain of R n (not necessary bounded), in case Ω has a boundary, the previous equation is supplemented with the no-flux boundary condition (ν denotes the outer unit normal to ∂Ω):

(∇P (ρ) + U [ρ]ρ) • ν = 0 on ∂Ω. (3.2)
For every ρ ∈ P(Ω), we assume that the Helmholtz decomposition of the vector field U [ρ]:

U [ρ] = -W [ρ] + ∇V [ρ], (3.3)
with

∇ • W [ρ] = 0, W [ρ] • ν = 0 on ∂Ω,
satisfies the following assumptions.

Assumptions on the potential part V :

• ∇V [ρ] ∈ L ∞ loc uniformly in ρ i.e for all K ⊂⊂ Ω, there exists C > 0 such that for all ρ ∈ P(Ω),

∇V [ρ] ∞,K C. (3.4)
Note that by Rademacher's condition, this condition implies that V [ρ] is differentiable a.e.,

• V [ρ] is semi-convex uniformly in ρ i.e there exists C such that for all ρ ∈ P(Ω), for every y ∈ Ω and every x ∈ Ω, point of differentiability of V [ρ]:

V [ρ](y) V [ρ](x) + ∇V [ρ](x), y -x - C 2 |y -x| 2 , (3.5)
• there exists C 0 such that for all ρ ∈ P(Ω), for all x ∈ Ω:

V [ρ](x) -C(1 + |x|), (3.6)
• ∇V [ρ] ∈ L 2 (ρ) uniformly in ρ, i.e there exists C > 0 such that for all ρ ∈ P(Ω),

ˆΩ |∇V [ρ]| 2 dρ C, (3.7)
• There exists C > 0 such that for all ρ, µ ∈ P(Ω),

ˆΩ |∇V [ρ] -∇V [µ]| 2 dρ CW 2 2 (ρ, µ), (3.8)
Assumptions on the divergence-free part W :

• there exists C > 0 such that for all ρ ∈ P(Ω) and

W [ρ] ∈ W 1,1 loc (R n ) and |W [ρ](x)| C(1 + |x|), for all x ∈ R n (3.9)
• There exists C > 0 such that for all ρ, µ ∈ P(Ω),

ˆΩ |W [ρ] -W [µ]| 2 dρ CW 2 2 (ρ, µ), (3.10)
Assumptions on the internal energy F and the associated pressure P :

The nonlinear diffusion term is given by a continuous strictly convex superlinear (i.e.

F (ρ)/ρ → +∞ as ρ → +∞) function F : R + → R of class C 2 ((0, +∞)) which satisfies F (0) = 0, and P (ρ) C(ρ + F (ρ)). (3.11)
where P (ρ) := ρF ′ (ρ) -F (ρ) is the pressure associated to F . Moreover, we define F : P(Ω) → R by

F (ρ) := ´Ω F (ρ(x)) dx if ρ ≪ L n , +∞ otherwise.
And we assume that

F (ρ) -C(1 + M (ρ)) α , for all ρ ∈ P(Ω), (3.12)
where α ∈ (0, 1) and M (ρ) := ´Ω |x| 2 dρ(x) is the second moment of ρ.

The typical examples of energies we have in mind are F (ρ) := ρ log(ρ), which gives a linear diffusion driven by the laplacian, and F (ρ) := ρ m (m > 1) which corresponds to the porous medium equation.

A weak solution of (3.1)

-(3.2) is a curve ρ : t ∈ (0, +∞) → ρ(t, •) ∈ P ac 2 (Ω) such that ∇P (ρ) ∈ M n ([0, +∞)× Ω) and ˆ∞ 0 ˆΩ(∂ t φρ -∇φ • U [ρ]ρ)dx -ˆΩ ∇φ • d∇P (ρ) dt = -ˆΩ φ(0, x)ρ 0 (x) dx, (3.13)
for every φ ∈ C ∞ c ([0, +∞) × R n ).
Our main result is the following:

Theorem 3.1. Assume ρ 0 ∈ P ac 2 (Ω) such that F (ρ 0 ) < +∞, (3.14)
then (3.1) admits at least one weak solution.

The proof of this theorem is given in the next sections and is based on the following splitting scheme that combines pure transport steps by the divergence-free part of the drift U and Wasserstein gradient flow steps taking into account the potential V in a semi-implicit way. More precisely, given a time step h > 0, we construct by induction a sequence ρ k h ∈ P ac 2 (Ω) by setting ρ 0 h = ρ 0 and given ρ k h we find ρ k+1 h using the following scheme:

• pure transport phase: we introduce an intermediate measure, ρk+1 h (with ρ0 h = ρ 0 ) defined by

ρk+1 h = X k h (h, •) # ρ k h , (3.15)
where X k h is solution of

∂ t X k h = W [ρ k h ] • X k h , X k h (0, •) = id . (3.16)
Since W [ρ k h ] satisfies (3.9), as recalled in section 2, DiPerna-Lions theory [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF] implies that X k h is well defined. Moreover, since W [ρ k h ] is divergence-free then X k h preserves the Lebesgue measure and leaves the domain Ω invariant thanks to the fact that W [ρ] is tangential to ∂Ω. Therefore ρk+1

h = ρ k h (X k h -1
) which implies the conservation of the internal energy:

F (ρ k+1 h ) = ˆΩ F (ρ k h (X k h -1 (x))) dx = ˆΩ F (ρ k h (x)) dx = F (ρ k h ) (3.17)
In addition, we can see ρk+1

h
is the value at time h of the solution µ of the continuity equation

∂ t µ + div(µW [ρ k h ]) = 0, µ |t=0 = ρ k h . (3.18)
Thanks to these observations, we can easily control the W 2 -distance between ρk+1 h and ρ k h . Indeed, using Benamou-Brenier formula and (3.9), we obtain

W 2 2 (ρ k+1 h , ρ k h ) h ˆh 0 ˆΩ |W [ρ k h ]| 2 dµ t dt Ch ˆh 0 ˆΩ(1 + |x| 2 ) dµ t dt Ch ˆh 0 (1 + M (µ t )) dt Moreover, d dt M (µ t ) = ˆΩ |x| 2 ∂ t µ t = -ˆΩ |x| 2 div(W [ρ k h ]µ t ) = 2 ˆΩ x • W [ρ k h ]µ t C(M (µ t ) + 1).
We obtain the last line using (3.9), Cauchy-Schwarz inequality and Young's inequality. Then,

M (µ t ) C(t + 1)e Ct 2Ce C for t 1, which implies W 2 2 (ρ k+1 h , ρ k h ) Ch 2 . (3.19)
• semi-implicit JKO scheme: In the second step we use a semi-implict version of the Jordan-Kinderlehrer-Otto scheme [START_REF] Jordan | The Variational Formulation of the Fokker-Plank Equation[END_REF], introduced by Di Francesco and Fagioli in [START_REF] Francesco | Measure solutions for nonlocal interaction PDEs with two species[END_REF] and used in [START_REF] Laborde | On some non linear evolution systems which are perturbations of Wasserstein gradient flows[END_REF], with ρk+1 h being the measure defined in the previous step. More precisely, we select ρ k+1 h as a solution of

inf ρ∈P ac 2 (Ω) E h (ρ|ρ k+1 h ) := W 2 2 (ρ, ρk+1 h ) + 2h F (ρ) + V(ρ|ρ k+1 h ) , (3.20)
where

V(ρ|µ) := ˆΩ V [µ] dρ.
By standard compactness and lower semicontinuity argument, (3.20) admits at least one solution (see for example [START_REF] Jordan | The Variational Formulation of the Fokker-Plank Equation[END_REF][START_REF] Laborde | On some non linear evolution systems which are perturbations of Wasserstein gradient flows[END_REF]) so the sequence ρ k h is well defined (it is even actually unique by strict convexity of E h (.|ρ k+1 h )). To summarize, given a time step h > 0, we construct by induction two sequences ρ k h and ρk h with the following splitting scheme: ρ 0 h = ρ0 h = ρ 0 and for all k 0,

ρk+1 h = X k h (h, •) # ρ k h , ρ k+1 h ∈ argmin ρ∈P ac 2 (Ω) W 2 2 (ρ, ρk+1 h ) + 2h F (ρ) + V(ρ|ρ k+1 h ) .
(3.21)

We finally introduce three different interpolations:

• We denote ρ h the usual piecewise constant interpolation of the sequence

ρ k h ρ h (t, •) := ρ k+1 h if t ∈ (hk, h(k + 1)], (3.22)
• similarly, we interpolate in a piecewise constant way the sequence ρk h :

ρ1 h (t, •) := ρk+1 h if t ∈ (hk, h(k + 1)], (3.23)
• finally, we denote by ρ2 h the continuous interpolation of ρk

h ρ2 h (t, •) := X k h (t -hk, •) # ρ k h if t ∈ (hk, h(k + 1)]. (3.24)
We remark that on (hk, h(k + 1)], ρ2 h is the solution on (0, h) of the continuity equation (3.18). The next two sections are devoted to the proof of theorem 3.1. In section 4, we derive various estimates on the sequences generated by the splitting scheme above, in particular thanks to the Euler-Lagrange equation of the semi-implicit JKO steps. This enables us to pass to the limit as the time step goes to 0 (the difficult term being of course the nonlinear pressure term) and thus to conclude the existence proof, this is done in section 5.

Estimates

Basic a priori estimates

Using the semi-implicit JKO scheme we first obtain the following a priori estimates on ρ h , ρ1 h and ρ2 h . Proposition 4.1. There exists h 0 > 0, such that for T > 0, there exists C > 0 such that, for all h, k, with h ∈ (0, h 0 ) and hk < T , N = ⌈ T h ⌉, we have

M (ρ k h ) C, (4.1)
F (ρ k h ) C, (4.2)
N -1 k=0 W 2 2 (ρ k+1 h , ρ k+1 h ) Ch. (4.3)
Proof. Using ρk+1 h as a competitor of ρ k+1 h in (3.20), we obtain

1 2h W 2 2 (ρ k+1 h , ρk+1 h ) F (ρ k+1 h ) -F (ρ k+1 h ) + ˆΩ V [ρ k+1 h ] ρk+1 h -ρ k+1 h . (4.4)
Let γ be the optimal transport plan between ρk+1 h and ρ k+1 h . Then we have

ˆΩ V [ρ k+1 h ] ρk+1 h -ρ k+1 h = ˆΩ V [ρ k+1 h ](x) -V [ρ k+1 h ](y) dγ(x, y) = ˆΩ V [ρ k+1 h ](x) -V [ρ k+1 h ](y) + ∇V [ρ k+1 h ](x) • (y -x) dγ(x, y) -ˆΩ ∇V [ρ k+1 h ](x) • (y -x) dγ(x, y).
Using (3.5) for the first part and Cauchy-Schwarz inequality and (3.7) for the second part of the right hand side, we find

ˆΩ V [ρ k+1 h ] ρk+1 h -ρ k+1 h C 2 ˆΩ |x -y| 2 γ(x, y) + C ˆΩ |x -y| 2 γ(x, y) 1 2 = C 2 W 2 2 (ρ k+1 h , ρ k+1 h ) + CW 2 (ρ k+1 h , ρ k+1 h ).
Choosing h ≤ h 0 ≤ 1 2C and using Young's inequality

W 2 (ρ k+1 h , ρ k+1 h ) ≤ 1 8hC W 2 2 (ρ k+1 h , ρ k+1 h ) + 2Ch, (4.4) becomes 1 8h W 2 2 (ρ k+1 h , ρk+1 h ) F (ρ k+1 h ) -F (ρ k+1 h ) + Ch.
Now using (3.17), to recover a telescopic sum, and summing over k, we obtain

N -1 k=0 W 2 2 (ρ k+1 h , ρk+1 h ) 8h F (ρ 0 ) -F (ρ N h ) + CT ,
this inequality and (3.14) imply (4.2). In addition, since the lower bound of F is controlled by the second moment,

N -1 k=0 W 2 2 (ρ k+1 h , ρk+1 h ) 8h F (ρ 0 ) + C(1 + M (ρ N h )) α + CT . (4.5)
But, with (4.5) and by standard arguments (see [START_REF] Jordan | The Variational Formulation of the Fokker-Plank Equation[END_REF][START_REF] Laborde | On some non linear evolution systems which are perturbations of Wasserstein gradient flows[END_REF]), we deduce that M (ρ k h ) satisfies (4.1) and then (4.5), (3.14) 3), we also have

N -1 k=0 W 2 2 (ρ k h , ρ k+1 h ) Ch and N -1 k=0 W 2 2 (ρ k h , ρk+1 h ) Ch.
Moreover, using (3.17), we have for all t ∈ [0, T ],

F (ρ 1 h (t)), F (ρ 2 h (t)) C and M (ρ 1 h (t)), M (ρ 2 h (t)) C.

Discrete Euler-Lagrange equation and stronger estimates

Let us start with the Euler-Lagrange equation of (3.20).

Proposition 4.3. For all k 0, we have

P (ρ k+1 h ) ∈ W 1,1 (Ω) and h ∇V [ρ k+1 h ]ρ k+1 h + ∇P (ρ k+1 h ) = -∇ϕ k+1 h ρ k+1 h a.e, (4.6)
where ϕ k+1 h is a Kantorovich potential from ρ k+1 h to ρk+1 h (so that its gradient is unique ρ k+1 h -a.e.) for W 2 . Proof. The proof is the same as in [START_REF] Agueh | Existence of solutions to degenerate parabolic equations via Monge-Kantorovich theory[END_REF][START_REF] Laborde | On some non linear evolution systems which are perturbations of Wasserstein gradient flows[END_REF] for example. We start by taking the first variation in the semi-implicit JKO scheme along the flow of a smooth vector field. Let ξ ∈ C ∞ c (Ω; R n ) be given and Φ τ the corresponding flow defined by

∂ τ Φ τ = ξ • Φ τ , Φ 0 = id .
We define a pertubation of ρ k+1 h by ρ τ := Φ τ # ρ k+1 h . Then we get

1 τ E h (ρ τ |ρ k+1 h ) -E h (ρ k+1 h |ρ k+1 h ) 0. (4.7)
By standard computations, we have

lim sup τ ց0 1 τ (W 2 2 (ρ τ , ρk+1 h ) -W 2 2 (ρ k+1 h , ρk+1 h )) ˆΩ×Ω (x -y) • ξ(x) dγ k+1 h (x,
1 τ (F (ρ τ ) -F (ρ k+1 h )) -ˆΩ P (ρ k+1 h (x)) div(ξ(x)) dx. (4.9)
Finally,

lim sup τ ց0 1 τ (V(ρ τ |ρ k+1 h ) -V(ρ k+1 i,h |ρ k+1 h )) ˆΩ ∇V [ρ k+1 h ] • ξρ k+1 h dx, (4.10)
If we combine (4.7), (4.8), (4.9) and (4.10), and if we replace ξ by -ξ, we find that, for all

ξ ∈ C ∞ c (Ω; R n ), ˆΩ ∇ϕ k+1 h • ξρ k+1 h -h ˆΩ P (ρ k+1 h ) div(ξ) + h ˆΩ ∇V [ρ k+1 h ] • ξρ k+1 h = 0. (4.11)
Now we claim that P (ρ k+1 h ) ∈ W 1,1 (Ω). Indeed, since P is controlled by F thanks to assumption (3.11), (4.2) gives P (ρ k+1 h ) ∈ L 1 (Ω). Moreover, using (4.11), we obtain

ˆΩ P (ρ k+1 h ) div(ξ) ˆΩ |∇ϕ k h (y)| h ρ k+1 h + ˆΩ |∇V [ρ k+1 h ]|ρ k+1 h ξ L ∞ (Ω) .
But using (3.8), (3.7), (4.3) and Cauchy-Schwarz inequality we get

ˆΩ |∇V [ρ k+1 h ]|ρ k+1 h ˆΩ |∇V [ρ k+1 h ] -∇V [ρ k+1 h ]|ρ k+1 h + ˆΩ |∇V [ρ k+1 h ]|ρ k+1 h ˆΩ |∇V [ρ k+1 h ] -∇V [ρ k+1 h ]| 2 ρ k+1 h 1/2 + ˆΩ |∇V [ρ k+1 h ]| 2 ρ k+1 h 1/2 C W 2 (ρ k+1 h , ρ k+1 h ) + 1 C.
We thus have

ˆΩ P (ρ k+1 h ) div(ξ) W 2 (ρ k+1 h , ρ k+1 h ) h + C ξ L ∞ (Ω) .
This implies P (ρ k+1 h ) ∈ BV (Ω) and ∇P (ρ

k+1 h ) = -∇V [ρ k+1 h ]ρ k+1 h - ∇ϕ k+1 h h ρ k+1 h in M n (Ω). In fact, P (ρ k+1 h ) is in W 1,1 (Ω) because ∇V [ρ k+1 h ]ρ k+1 h + ∇ϕ k h h ρ k+1 h ∈ L 1
(Ω) and we have proved (4.6).

We immediately deduce an L 1 ((0, T ), BV (Ω)) estimate for P (ρ h ):

Corollary 4.4. For all T > 0, we have

P (ρ h ) L 1 ((0,T );W 1,1 (Ω)) CT. (4.12)
Proof. If we integrate (4.6), we obtain

h ˆΩ |∇P (ρ k+1 h )| W 2 (ρ k+1 h , ρk+1 h ) + Ch,
Then we sum from k = 0 to N -1 and thanks to (4.3), we have

ˆT 0 ˆΩ |∇P (ρ h )| CT.
We conclude thanks to (3.11) and (4.2).

Proposition 4.5. Let h > 0, N ∈ N * , T := N h, t k := hk, for k = 0, • • • , N , then, for every φ ∈ C ∞ c ([0, T )×R n ) ˆT 0 ˆΩ ρ2 h (t, x)(∂ t φ(t, x) + W [ρ h (t -h)] • ∇φ(t, x)) dxdt = h N -1 k=0 ˆΩ ∇P (ρ k+1 h (x)) • ∇φ(t k , x) dx + h N -1 k=0 ˆΩ ∇V [ρ k+1 h ] • ∇φ(t k , x)ρ k+1 h dx + N -1 k=0 ˆΩ×Ω R[φ(t k , •)](x, y)dγ k+1 h (x, y) -ˆΩ ρ 0 (x)φ(0, x) dx, with, for all φ ∈ C ∞ c ([0, T ) × R n ), |R[φ](x, y)| 1 2 D 2 φ L ∞ ([0,T )×Ω) |x -y| 2 ,
and γ k+1 h is the optimal transport plan in

Π(ρ k+1 h , ρk+1 h ). Proof. Let ϕ ∈ C ∞ c (R n ), multiplying (4.6
) by ∇ϕ and integrating on Ω, we obtain

-ˆΩ ∇ϕ k+1 h • ∇ϕρ k+1 h = h ˆΩ ∇P (ρ k+1 h ) • ∇ϕ + ˆΩ ∇V [ρ k+1 h ] • ∇ϕρ k+1 h .
But, we can rewrite the left hand side by

-ˆΩ ∇ϕ k+1 h • ∇ϕρ k+1 h = ˆΩ×Ω (y -x) • ∇ϕ(x) dγ k+1 h (x, y).
A second-order Taylor-Lagrange formula then gives

ˆΩ×Ω (y -x) • ∇ϕ(x) dγ k+1 h = ˆΩ×Ω (ϕ(y) -ϕ(x)) dγ k+1 h (x, y) - ˆΩ×Ω R[ϕ](x, y)dγ k+1 h (x, y) = ˆΩ ϕ(ρ k+1 h -ρ k+1 h ) - ˆΩ×Ω R[ϕ](x, y)dγ k+1 h (x, y). (4.13) Now let φ ∈ C ∞ c ([0, T ) × R n ), we have ˆT 0 ˆΩ ρ2 h (t, x)(∂ t φ(t, x) + W [ρ h (t -h)] • ∇φ(t, x)) dxdt = N -1 k=0 ˆtk+1 t k ˆΩ ρ k h (x)(∂ t φ + W [ρ k h ] • ∇φ)(t, X k h (t -t k , x)) dxdt But, on [t k , t k+1 ], d dt [φ(t, X k h (t -t k , x))] = (∂ t φ + W [ρ k h ] • ∇φ)(t, X k h (t -t k , x)).
Then,

N -1 k=0 ˆtk+1 t k ˆΩ ρ k h (x)(∂ t φ + W [ρ k h ] • ∇φ)(t, X k h (t -t k , x)) dxdt = N -1 k=0 ˆtk+1 t k ˆΩ ρ k h (x) d dt [φ(t, X k h (t -t k , x))] dxdt = N -1 k=0 ˆΩ ρ k h (x) φ(t k+1 , X k h (h, x)) -φ(t k , x) dx = N -1 k=0 ˆΩ φ(t k+1 , x)ρ k+1 h (x) -φ(t k , x)ρ k h (x) dx = N -1 k=0 ˆΩ φ(t k+1 , x) ρk+1 h (x) -ρ k+1 h (x) -ˆΩ φ(0, x)ρ 0 (x) dx.
Then the proof is complete by applying (4.13) with ϕ = φ(t k+1 , •).

5 Convergence and proof of Theorem 3.1

Weak and strong convergences

Using a refined version of Ascoli theorem (see [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF]), estimate (4.3) and remark 4.2 and taking subsequences, if necessary, we have that, for every T < +∞, ρ h , ρ1 h and ρ2 h converge in L ∞ ((0, T ), W 2 ) to some respective limits ρ, ρ1 and ρ2 :

sup t∈[0,T ] max(W 2 (ρ h (t, .), ρ(t, .)), W 2 (ρ 1 h (t, .), ρ1 (t, .)), W 2 (ρ 2 h (t, .), ρ2 (t, .))) → 0 as h → 0.
In fact, these three sequences have to converge to the same limit, ρ. Indeed, for all ϕ ∈ C ∞ c ((0, T ) × Ω),

ˆT 0 ˆΩ ϕ(ρ h -ρ1 h ) = N -1 k=0 ˆtk+1 t k ˆΩ ϕ(ρ k+1 h -ρk+1 h ) Ch N -1 k=0 W 2 (ρ k+1 h , ρk+1 h ) ChN 1/2 N -1 k=0 W 2 2 (ρ k+1 h , ρk+1 h ) 1/2 CT 1/2 h,
because of (4.3). With a similar computation, we find that ρ h and ρ2 h converge to the same limit. We thus have

sup t∈[0,T ] max(W 2 (ρ h (t, .), ρ(t, .)), W 2 (ρ 1 h (t, .), ρ(t, .)), W 2 (ρ 2 h (t, .), ρ(t, .))) → 0 as h → 0. (5.1)
Moreover it is classical to deduce from (4.3) and remark 4.2 an Hölder-like estimate of the form W 2 (ρ h (t, .), ρ h (s, .)) ≤ C |t -s| + h from which one deduces that the limit curve ρ actually belongs to C 1/2 ((0, T ), W 2 ). This kind of convergence will be enough to pass to the limit in ∇V [ρ 1 h ]ρ h and W [ρ h ]ρ 2 h , because of assumptions (3.8) and (3.10), but we will need a stronger convergence to deal with the nonlinear diffusion term P (ρ h ). Fo this purpose, we will use an extension of the Aubin-Lions Lemma due to Rossi and Savaré in [START_REF] Rossi | tightness, integral equicontinuity and compactness for evolution problem in Banach spaces[END_REF]: Theorem 5.1 (th. 2 in [START_REF] Rossi | tightness, integral equicontinuity and compactness for evolution problem in Banach spaces[END_REF]). On a Banach space B, let be given • a normal coercive integrand G : B → R + , i.e, G is l.s.c and its sublevels are relatively compact in B,

• a pseudo-distance g : B × B → [0, +∞], i.e, g is l.s.c and [g(ρ, µ) = 0, ρ, µ ∈ B with G(ρ), G(µ) < ∞] ⇒ ρ = µ.
Let T > 0 and U be a set of measurable functions u : (0, T ) → B. Under the hypotheses that

sup u∈U ˆT 0 G(u(t)) dt < +∞ and lim hց0 sup u∈U ˆT -h 0 g(u(t + h), u(t)) dt = 0, (5.2)
U contains a subsequence (u n ) n∈N which converges (strongly in B) in measure with respect to t ∈ (0, T ) to a limit u ⋆ : (0, T ) → B.

We now apply this theorem to B = L 1 (Ω), U = {ρ h } h , g defined by

g(ρ, µ) := W 2 (ρ, µ) if ρ, µ ∈ P 2 (Ω), +∞ otherwise,
and G by

G(ρ) := F (ρ) + P (ρ) BV (Ω) + M (ρ) if ρ ∈ P ac 2 (Ω), P (ρ) ∈ BV (Ω) and F (ρ) ∈ L 1 (Ω), +∞ otherwise.
Lemma 5.2. G is l.s.c on L 1 (Ω) and its sublevels are relatively compact in L 1 (Ω).

Proof. Let us start by proving that sublevels of G are relatively compact in L 1 (Ω). Let

A c := ρ ∈ L 1 (Ω) : G(ρ) c
and (ρ k ) be a sequence in A c then P (ρ k ) is bounded in BV (Ω) thus, up to a subsequence P (ρ k ) converges to some Φ in L 1 loc (Ω) and a.e.. Since P is continuous, one to one and its inverse is continuous, ρ k converges to ρ := P -1 (Φ) a.e.; and, since G(ρ k ) c and F is superlinear, ρ k is uniformly integrable, using Vitali's convergence theorem, we obtain that ρ k converges to ρ in L 1 (K ∩ Ω) for every compact K. To conclude that there is convergence in L 1 (Ω), we use the fact that the second momentum of ρ k is uniformly bounded:

ˆΩ |ρ k -ρ| ˆΩ\BR |x| 2 R 2 |ρ k -ρ| + ˆBR∩Ω |ρ k -ρ| 2c R 2 + ˆBR∩Ω |ρ k -ρ|.
The first term in the right hand can be made arbitrary small by choosing R large enough and the second term converges to zero by L 1 (B R ∩ Ω)-convergence. Now we have to show the lower semi-continuity of G on L 1 (Ω). Let (ρ k ) be a sequence which converges strongly to ρ in L 1 (Ω) with (without loss of generality) sup k G(ρ k ) C < +∞. Without loss of generality, we can assume that ρ k converges to ρ a.e. Since sup k G(ρ k ) C, P (ρ k ) is uniformly bounded in BV (Ω) so P (ρ k ) converges weakly to µ in BV (Ω). Moreover, P (ρ k ) converges strongly to µ in L 1 loc (Ω). We can conclude that µ = P (ρ) and by lower semi-continuity of F , M and the BV -norm we have

G(ρ) lim inf kր+∞ G(ρ k ).
Thanks to lemma 5.2, to apply theorem 5.1, it remains to verify (5.2). The first condition of (5.2) is satisfied because of the estimate on the momentum, (4.1), on the internal energy F , (4.2) and on the gradient of P (ρ h ) (4.12). The second condition of (5.2) comes from the estimate on the distance (4.3) and remark 4.2 (see for example [START_REF] Francesco | curves of steepest descent are entropy solutions for a class of degenerate convection-diffusion equations[END_REF] for a detailed proof). Then theorem 5.1 implies that ρ h converges in measure with respect to t in L 1 (Ω) to ρ. Since convergence in measure implies a.e convergence, up to a subsequence, we may also assume that ρ h (t, .) converges strongly in L 1 (Ω) to ρ(t, .) for a.e. t. Then Lebesgue's dominated convergence theorem implies that ρ h converges strongly in L 1 ((0, T ) × Ω) to ρ.

Thanks to (3.11) and (4.2) P (ρ h ) is uniformly bounded in L ∞ ((0, T ), L 1 (Ω)). In addition, by corollary 4.4, P (ρ h ) is uniformly bounded in L 1 ((0, T ), W 1,1 (Ω)). Thanks to the Sobolev embedding, we deduce that P (ρ h ) is uniformly bounded in L ∞ ((0, T ), L 1 (Ω)) ∩ L 1 ((0, T ), L n/n-1 (Ω)). To have uniform integrability of P (ρ h ) both in the time and space variables, the following will be useful: Lemma 5.3. Let p > 1, q := 2p-1 p and f ∈ L ∞ ((0, T ), L 1 (Ω)) ∩ L 1 ((0, T ), L p (Ω)) then f ∈ L q ((0, T ) × Ω) and we have f q L q ((0,T )×Ω)

f q-1 L ∞ t (L 1 x ) f L 1 t (L p x ) . Proof. Writing 1 q = p 2p -1 = θ + 1 -θ p , with θ = p -1 2p -1 ∈ (0, 1)
and observing that (1 -θ)q = 1 and θq = q -1, the interpolation inequality yields that

f q L q x f q-1 L 1 x f L p x , since f ∈ L ∞ ((0, T ), L 1 (Ω)) ∩ L 1 ((0, T ), L p (Ω)) this implies that f ∈ L q ((0, T ) × Ω) and f q L q ((0,T )×Ω) = ˆT 0 f q L q x = ˆT 0 f q-1 L 1 x f L p x f q-1 L ∞ t (L 1 x ) f L 1 t (L p x ) .
Applying lemma 5.3 we deduce that P (ρ h ) is uniformly bounded in L (n+1)/n ((0, T ) × Ω). This implies that P (ρ h ) is uniformly integrable and since we know that it converges a.e. to P (ρ), we can deduce from Vitali's convergence theorem that P (ρ h ) converges strongly to P (ρ) in L 1 loc ((0, T ) × Ω). Thanks to Corollary 4.4 we deduce that ∇P (ρ h ) converges vaguely to ∇P (ρ) in M n loc ((0, T ) × Ω). In fact, we have ∇P (ρ) in M n ((0, T ) × Ω) and narrow convergence of ∇P (ρ h ) to ∇P (ρ), thanks to Prokhorov Theorem and the following tightness estimate: Lemma 5.4. The family ∇P (ρ h ), viewed as vector-valued measures on [0, T ] × Ω, is tight, more precisely, for every h and every A measurable, A ⊂ Ω

ˆT 0 ˆA |∇P (ρ h )| ≤ C 1 + √ h ˆT 0 ˆA ρ h (t, x)dxdt 1/2 . ( 5.3)
Proof. Integrating (4.6) on (0, T ) × A together with Cauchy Schwarz inequality and (4.3), we get (taking

N = ⌈ T h ⌉ + 1, say) ˆT 0 ˆA |∇P (ρ h )| N k=0 ˆA |∇ϕ k+1 h |ρ k+1 h + ˆT 0 ˆA |∇V [ρ 1 h ]|ρ h N k=0 ˆΩ |∇ϕ k+1 h | 2 ρ k+1 h 1/2 ˆA ρ k+1 h 1/2 + ˆT 0 ˆA |∇V [ρ 1 h ]|ρ h N k=0 W 2 2 (ρ k+1 h , ρ k+1 h ) 1/2 ˆT 0 ˆA ρ h (t, x)dxdt 1/2 + ˆT 0 ˆA |∇V [ρ 1 h ]|ρ h C √ h ˆT 0 ˆA ρ h (t, x)dxdt 1/2 + ˆT 0 ˆA |∇V [ρ 1 h ]|ρ h .
Moreover, with Cauchy Schwarz inequality, (3.7) and (3.8), we also have

ˆT 0 ˆA |∇V [ρ 1 h ]|ρ h ˆT 0 ˆA |∇V [ρ h ]|ρ h + ˆT 0 ˆA |∇V [ρ 1 h ] -∇V [ρ h ]|ρ h C 1 + √ h ˆT 0 ˆA ρ h (t, x)dxdt 1/2
which proves (5.3). The tightness of ∇P (ρ h ) therefore immediately follows from that of ρ h and (5.3).

We can summarize all of this in the next result:

Theorem 5.5. Up to a subsequence ρ h converges strongly in L 1 ((0, T ) × Ω), P (ρ h ) converges strongly to P (ρ) in L 1 loc ((0, T ) × Ω) and ∇P (ρ h ) converges to ∇P (ρ) narrowly in M((0, T ) × Ω).

End of the proof of theorem 3.1

In this section we finish the proof of theorem 3.1. We have to pass to the limit in all terms in proposition 4.5. The linear term (with time derivative) and the diffusion term converge to the desired result because ρ2 h converges to ρ in L ∞ ([0, T ], W 2 ) and ∇P (ρ h ) converges to ∇P (ρ) narrowly in M n ([0, T ] × Ω). Remainder term goes to 0 when h goes to 0 because of (4.3). So we just have to check the convergence in the transport terms, in what follows the test-function φ belongs again to

C ∞ c ([0, T ) × R n ). • term in W : We have to show: ˆT 0 ˆΩ ρ2 h (t, x)W [ρ h (t -h)](x) • ∇φ(t, x) dxdt → ˆT 0 ˆΩ ρ(t, x)W [ρ(t, •)](x) • φ(t, x) dxdt.
We first have

ˆT 0 ˆΩ ρ2 h (t, x)W [ρ h (t -h)](x) • ∇φ(t, x) dxdt = ˆT 0 ˆΩ ρ2 h (t, x)W [ρ 2 h (t)](x) • ∇φ(t, x) dxdt + ˆT 0 ˆΩ ρ2 h (t, x)(W [ρ h (t -h)](x) -W [ρ 2 h (t)](x)) • ∇φ(t, x) dxdt.
The second term in the right hand side goes to zero when h goes to 0. Indeed, the first term in the right hand-side converges to 0 because of (5.1). As for the second one, it also converges to 0, because W [ρ] • ∇φ belongs to L ∞ ((0, T ) × Ω) and ρ2 h is uniformly integrable by remark 4.2 and the superlinearity of F , hence, up to a subsequence it converges to ρ weakly in L 1 ((0, T ) × Ω).

ˆT 0 ˆΩ ρ2 h (t, x)(W [ρ h (t -h)](x) -W [ρ 2 h (t)](x)) • ∇φ(t, x) dxdt C ˆT 0 ˆΩ ρ2 h (t, x)|W [ρ h (t -h)](x) -W [ρ 2 h (t)(x)]| 2 dx 1/2 dt, then using (3.10), ˆT 0 ˆΩ ρ2 h (t, x)(W [ρ h (t -h)](x) -W [ρ 2 h (t)](x)) • ∇φ(t, x) dxdt C ˆT 0 W 2 (ρ h (t -h), ρ2 h (t)) dt C N -1 k=0 ˆtk+1 t k W 2 (ρ k h , X k h (t -t k ) # ρ k h ) dt CT h
• term in ∇V : We claim that The proof is the same as the previous one for W , using (3.8), (3.4) and the convergence of ρ h to ρ in L 1 ((0, T ) × Ω)) ∩ L ∞ ((0, T ), W 2 ).

On extension to systems and uniqueness

The splitting transport-JKO scheme described above, can easily be adapted, under suitable assumptions to the case of systems for the evolution of N species coupled by nonlocal drifts:

∂ t ρ i -∆P i (ρ i ) -div(ρ i U i [ρ 1 , • • • , ρ N ]) = 0, ρ i (0, .) = ρ i,0 , i = 1, • • • , N, (6.1)
where P i (s) = sF ′ i (s) -F i (s) is the pressure associated to a strictly convex superlinear function F i with corresponding internal energy F i (ρ i ) := ´Fi (ρ i (x))dx. Decomposing each drift

U i [ρ 1 , • • • , ρ N ] = ∇V i [ρ 1 , • • • , ρ N ] - W i [ρ 1 , • • • , ρ N ] with div(W i [ρ 1 , • • • , ρ N ]
) = 0 and under similar assumptions as in paragraph 3, one can show, by similar arguments as above, convergence as h → 0 to a solution of (6.1) of the following splitting scheme.

Starting form ρ 0 i,h = ρ i,0 and given ρ k h = (ρ k 1,h , • • • ρ k N,h ) we find ρ k+1 h = (ρ k+1 1,h , • • • ρ k+1 N,h ) by:

• setting ρk+1 i,h = X k i,h (h, .) # ρ k i,h where

∂ t X k i,h = W i [ρ k h ] • X k i,h , X k i,h (0, •) = id,
• defining ρk+1

h = (ρ k+1 1,h , • • • , ρk+1 N,h ), ρ k+1 h = (ρ k+1 1,h , • • • ρ k+1 N,h
) is obtained by the semi-implicit JKO scheme:

ρ k+1 i,h = argmin ρi∈P ac 2 (Ω) W 2 2 (ρ i , ρk+1 i,h ) + 2h F i (ρ i ) + ˆΩ V i [ρ k+1 h ]ρ i .
Finally, let us say a few words on uniqueness which we have not addressed here, but which can be obtained at least in two ways: either by assuming some displacement semiconvexity of the internal energy and proving some exponential in time contraction estimate on the W 2 distance between two solutions (see [START_REF] Francesco | Measure solutions for nonlocal interaction PDEs with two species[END_REF][START_REF] Laborde | On some non linear evolution systems which are perturbations of Wasserstein gradient flows[END_REF]), or by assuming some nondegeneracy of the diffusion and establishing some H -1 contraction estimate (see [START_REF] Carlier | Remarks on continuity equations with nonlinear diffusion and nonlocal drifts[END_REF]).

h N - 1 k=0ˆΩ

 1 ∇V [ρ k+1 h ](x) • ∇φ(t k , x)ρ k+1 h dx → ˆT 0 ˆΩ ∇V [ρ(t, •)](x) • ∇φ(t, x)ρ(t, x) dxdt.

  and (4.1) give (4.3).

	Remark 4.2. Using estimate (3.19) between ρk+1

h and ρ k h , and (4.


Cristóbal Bertoglio

Rodolfo Núñez

Felipe Galarce

David Nordsletten

Axel Osses

Relative pressure estimation from velocity measurements in blood flows: state-of-the-art and new approaches

The pressure gradient across stenotic blood vessels is an important clinical index for the diagnosis of many cardiovascular diseases. While clinical gold standard for its measurement (in terms of accuracy) is invasive catheterization, Phase-Contrast MRimaging has emerged as a promising tool for enabling a non-invasive quantification, by linking (highly spatially resolved) velocity measurements with relative pressures via Navier-Stokes equations. In this work we provide a review of current methods for relative pressure estimation, propose new ones and compare them on numerical examples using subsampled and noisy synthetic data. We verify that the newly proposed approaches are more robust with respect to data perturbations and therefore more precise.

Introduction

Pressure difference across stenotic blood vessels is a standard clinical measurement that serves to assess the severity of the pathology. It may be performed several times during diagnosis and follow up, for example, in congenital heart diseases [START_REF] Warnes | ACC/AHA 2008 guidelines for the management of adults with congenital heart disease[END_REF]. The most reliable procedure in clinical practice to obtain pressure gradients involves invasive X-ray guided catheterization, which should be avoided unless non-invasive evaluations (like ultrasound) are inconclusive or discordant with clinical findings [START_REF] Vahanian | Guidelines on the management of valvular heart disease (version 2012)[END_REF].

Alternatively, Phase-Contrast Magnetic Resonance Imaging (PC-MRI) allows to obtain three-dimensional, time-resolved measurements of velocity fields, allowing visualization of complex blood flow patterns in large vessels and the heart [START_REF] Markl | 4D flow MRI[END_REF]. Typical PC-MRI spatial and temporal resolutions go from 1.5 to 3 mm 3 and 20 -40 ms, respectively [START_REF] Dyverfeldt | 4D flow cardiovascular magnetic resonance consensus statement[END_REF].

Previous works have postulated methods to estimate relative pressures using these measured velocities and the Navier-Stokes equations. To the best of the authors knowledge, the oldest method in use is the so called Pressure Poisson Estimator (PPE), see for example [START_REF] Ebbers | Estimation of relative cardiovascular pressures using time-resolved three-dimensional phase contrast MRI[END_REF][START_REF] Ebbers | Improving computation of cardiovascular relative pressure fields from velocity MRI[END_REF][START_REF] Sebastian Bs Krittian | A finite-element approach to the direct computation of relative cardiovascular pressure from time-resolved MR velocity data[END_REF]. It consists in reconstructing the relative pressure field by directly plugging the velocity measurements into the Navier-Stokes equations (NSE), and solving a Poisson equation for the pressure.

As a new approach, recently in [START_REF] Švihlová | Determination of pressure data from velocity data with a view toward its application in cardiovascular mechanics. Part 1. Theoretical considerations[END_REF] the authors proposed to perturb the Navier-Stokes residual with an auxiliary divergence-free vector, resulting in a mixed problem for this vector and the pressure to be estimated. Therefore it was denoted as Stokes Estimator (STE). In the tests presented in [START_REF] Švihlová | Determination of pressure data from velocity data with a view toward its application in cardiovascular mechanics. Part 1. Theoretical considerations[END_REF], STE turned out to outperform the PPE in terms of precision.

Also very recently in [START_REF] Donati | Non-invasive pressure difference estimation from PC-MRI using the workenergy equation[END_REF], relative pressure and velocity measurements were proposed to be related by using the classical integral energy balance of the incompressible Navier-Stokes equations, denoted as WERP. This has the advantage that only few integrals have to be evaluated, making it computationally faster than PPE and STE. However, the authors indicated some potential drawbacks: (i) the velocity measurement terms are squared, amplifying the noise, in particular for small in-/outflows, and (ii) it cannot be applied for simultaneous estimation of relative pressure in multiple outlet geometries. Nevertheless, WERP also outperformed PPE for peak relative pressure estimates.

The purpose of this manuscript is twofold. First, to compare the methods recently proposed (STE and WERP). Secondly, to propose improvements of such methods. For the STE, we show that a simple integration by parts of the terms containing the velocity considerably improves the robustness regarding spatial subsampling. Then, inspired from the WERP, we formulate an integral approach based on the momentum balance for estimating relative pressures. We also derived formulae of the statistical bias of all methods, what to the best of the authors knowledge has not yet been reported. We confirm numerically that all methods, with exception of the WERP, have a negligible bias.

The rest of this paper is organized as follows. In Section 2 we set up the pressure estimation problem. In Section 3 we revisit the state-of-the-art methods mentioned above. In Section 4 we propose new pressure estimation approaches. In Section 5 we perform a theoretical bias analysis for all the methods. In Section 6 we test and compare all methods using synthetic data, including sensitivity to noise, spatial and temporal subsampling. Finally in Section 7 we discuss the results and draw some conclusions and perspectives. the pressure p(t) : Ω → R, t ∈ [0, T ]:

Γ i Γ w Γ w Γ o
ρ∂ t u + ρ(u • ∇)u -µ∆u + ∇p = 0 in Ω × [0, T ], ∇ • u = 0 in Ω × [0, T ] ( 1
)
where ρ is the density and µ the dynamic viscosity. For our purposes, Ω will be the domain around the stenosis. We also assume that the boundary of Ω is given by ∂Ω

:= Γ = Γ i ∪ Γ o ∪ Γ w , with Γ i ∩ Γ o = ∅. Γ i is the inlet boundary (i.
e. upstream of the obstruction), Γ o the outlet boundary (downstream of the obstruction) and Γ w is the arterial wall boundary, see example in Figure 1. The task of the relative pressure estimation is then to compute:

δp = 1 |Γ i | Γ i p - 1 |Γ o | Γo p (2)
given (perturbed) measurements of u in Ω.

As is obtained in Phase-Contrast MRI, we will assume that the velocity measurements are available at h-spaced discrete points of the domain Ω, but we will be able to reconstruct them everywhere by: (a) computing a mesh T h of the domain Ω using the measurement locations and a segmented geometry from the images, and then (b) interpolating the velocity components at every point in Ω using simple linear piecewise finite element basis functions. Additionally, we will asume that these space-discrete measurements will be available at N discrete measurement times t 1 , ..., t N at a constant time interval τ . We will denote the measurements as u 1 m , ..., u N m ∈ [P 1,h ] d , with P 1,h the usual piecewise continuous linear polynomial finite element space.

Moreover, for further analysis we assume that the measurements include a random, additive perturbation with respect to the ground truth velocity, i.e. u n m = u n h + ε n , with u n h ∈ [P 1,h ] d a spatial subsampling of the true field u into [P 1,h ] d at observation time t n , and ε n ∈ [P 1,h ] d a discrete noise field, with all degrees-of-freedom independently identically distributed (i.i.d.) following a normal distribution N (0, σ 2 ).

Notice that we assume only the reference field u to satisfy Problem (1) (in practice we will compute one from a highly resolved numerical simulation). The discrete measured velocity fields u n h are the result of a subsampling operator, and therefore they do not fulfill neither [START_REF] Bertoglio | A Stokes-residual backflow stabilization method applied to physiological flows[END_REF] nor a discretized version of it. Naturally, the same applies to u n m which includes additional random perturbations.

State of the art of estimation methods

In this section, we overview the mathematical formulations of the currently available methods to estimate relative pressure. We start with the Poisson Pressure Estimator (PPE), which is the oldest and most popular. We continue with the Stokes estimator (STE) and then with the Work-Energy Relative pressure estimator (WERP), both recently introduced. We remark that the later two approaches has been independently reported to have a superior performance than the PPE, but not have been yet compared with each other.

The Poisson Pressure Estimator (PPE)

The PPE, in the context of finite elements, is based on the assumption that the pressure gradient ∇p satisfies a weak version of Equation [START_REF] Bertoglio | A Stokes-residual backflow stabilization method applied to physiological flows[END_REF]. Hence, the estimated relative pressure at time t n+1/2 := (t n + t n+1 )/2 can be found as a solution of the following problem: Find

p n+1/2 ppe ∈ P 1,h such that Ω ∇p n+1/2 ppe • ∇q = - ρ τ Ω (u n+1 m -u n m ) • ∇q -ρ Ω (u n+1/2 m • ∇u n+1/2 m ) • ∇q (3)
for all q ∈ P 1,h , and p ppe = q = 0 on Γ o . Note that we leave out the term ∆u n+1/2 m since it was reported in [START_REF] Lamata | Aortic relative pressure components derived from four-dimensional flow cardiovascular magnetic resonance[END_REF] that the viscous part of the pressure gradient in large vessels is negligible. Moreover, since u n m ∈ P 1,h , the data does not possess enough regularity to compute its laplacian.

Then, the pressure drop at time t n+1/2 between the outlet and the inlet is evaluated with Equation (2) using p = p n+1/2 ppe . Note also that the midpoint time evaluation scheme allows a second order approximation with respect to time and to reduce the noise variance of the measurements u n+1/2 m by a factor of 2 with respect to the original variance of the measurements u n m .

The Stokes estimator (STE)

The STE formulated in [START_REF] Švihlová | Determination of pressure data from velocity data with a view toward its application in cardiovascular mechanics. Part 1. Theoretical considerations[END_REF] consists in perturbing the Navier-Stokes equation with the Laplacian of an auxiliary divergence-free velocity field w. This leads to a mixed problem for the pressure, in terms of the measured velocity, instead of a Poisson problem as in the PPE.

Defining P b 1,h as the finite element space with linear piecewise polynomials with an additional bubble function at the centroid of the element, the STE is formulated as follows for the time step t n+1/2 : Find w ∈ [P b 1,h ] d and p n+1/2 ste ∈ P 1,h , such that

Ω ∇w : ∇ ŵ - Ω p n+1/2 ste (∇ • ŵ) + Ω (∇ • w)q = - ρ τ Ω (u n+1 m -u n m ) • ŵ -ρ Ω (u n+1/2 m • ∇u n+1/2 m ) • ŵ (4)
for all ŵ ∈ [P b 1,h ] d and q ∈ P 1,h , with w = ŵ = 0 on Γ. Then, the pressure drop at time t n+1/2 between the outlet and the inlet is evaluated with Formula (2) using p = p n+1/2 ste . We remark that no additional information on the pressure field p n+1/2 ste is needed for ensuring solvability of Problem [START_REF] Dyverfeldt | 4D flow cardiovascular magnetic resonance consensus statement[END_REF]. Notice also that, as for the PPE, the viscous term is neglected.

The work-energy relative pressure estimator (WERP)

The starting point of the formulation of the WERP method introduced in [START_REF] Donati | Non-invasive pressure difference estimation from PC-MRI using the workenergy equation[END_REF] is the classical energy relation of an incompressible Newtonian fluid, which can be obtained either by writing the conservation law for the energy density ρu 2 /2 directly or multiplying (1) by u and integrating over Ω. This relation reads

ρ 2 Ω ∂ t (|u| 2 ) + ρ 2 Γ (u • n)|u| 2 + µ Ω |∇u| 2 + Γ p (u • n) -µ Γ (∇u • n) • u = 0 , (5)
with | • | denoting the Euclidian norm.

The WERP is then formulated by assuming that: (i) the measurements u n m satisfy relation ( 5), (ii) u n m • n ≈ 0 on Γ w (i.e. the vessel walls nearly don't move), (iii) the viscous forces on Γ are negligible, (iv) and the pressure is nearly constant on Γ i and Γ o . Doing so, the WERP pressure estimator at time t n+1/2 is written as:

δp n+1/2 werp (u n,n+1 m ) = -1 Λ(u n+1/2 m ) E kin (u n+1 m ) -E kin (u n m ) + E conv (u n+1/2 m ) + E visc (u n+1/2 m ) (6)
with

E kin (w) = ρ 2τ Ω |w| 2 (7)
E conv (w) = ρ 2 Γ i ∪Γo (w • n)|w| 2 (8)
E visc (w) = µ Ω |∇w| 2 (9) Λ(w) = Γ i w • n. (10)
Note that the superscript n, n + 1 means that both velocity at the measured time steps n and n + 1 are considered.

As it is well known from numerical analysis, the use of the mid-point scheme in (6) does not induce any perturbation compared with the time-continuous energy relation [START_REF] Ebbers | Improving computation of cardiovascular relative pressure fields from velocity MRI[END_REF], independently on the time step τ .

We can see that this formulation may be unstable at small in-/outflows due to the division by Λ(u n+1/2 m

). In other words, when the inflow is small, perturbations (e.g. due to the presence of measurement noise) can induce large errors in Λ(u n+1/2 m ) leading to an unphysical amplification of the estimated relative pressure. Fortunately, low flow regimes are not relevant in clinical practice since only peak relative pressures are of interest. Another important issue is that the WERP induces systematic shifts in the pressure curve in the presence of noise, as we will see in the next section. Moreover, the WERP can only by construction estimate relative pressure in geometries with only one inflow and one outflow, due to the assumption u n m • n ≈ 0 on Γ w needed to eliminate the pressure field on Γ w out of the formulation.

New estimation methods

In this section we present new pressure estimation methods, inspired from the state-ofthe-art approaches, and we will see later in the numerical examples that they can lead to improvements in the estimation results.

The integrated STE (STEint)

We propose to modify the STE method originally proposed in [START_REF] Švihlová | Determination of pressure data from velocity data with a view toward its application in cardiovascular mechanics. Part 1. Theoretical considerations[END_REF] and formulated in Equation ( 4) by taking advantage of the regularity of the auxiliary velocity field and integrating the convective and viscous terms by parts.

The STEint is formulated as follows for the time step t n+1/2 . Find w

∈ [P b 1,h ] d and p n+1/2 stei ∈ P 1,h , such that Ω ∇w : ∇ ŵ - Ω p n+1/2 stei (∇ • ŵ) + Ω (∇ • w)q = - ρ τ Ω (u n+1 m -u n m ) • ŵ +ρ Ω (u n+1/2 m • ∇ ŵ) • u n+1/2 m -µ Ω ∇u n+1/2 m : ∇ ŵ ( 11
)
for all ŵ ∈ [P b 1,h ] d and q ∈ P 1,h , with w = ŵ = 0 on Γ. Notice that the boundary terms of the integration by parts do not appear since ŵ = 0. Then, the pressure drop at time t n+1/2 between the outlet and the inlet is evaluated with Formula (2) using p = p n+1/2 stei . The advantage of reducing the derivative order of the measurements will be noticeable in the robustness of the pressure estimation results when spatially subsampling the data. Unfortunately, this integration by parts is not straightforward for the PPE method due to the lack of enough regularity of ∇q in Equation (3) when standard piecewise polynomial finite element spaces are used (like P 1,h ).

The Darcy estimator (DAE)

Note that, a similar alternative for the STE would be to use a Darcy formulation for w instead of Stokes. Let RT k h be the usual Raviart-Thomas finite element space of degree k. Then, find w ∈ RT 1 h and p n+1/2 dae ∈ P 1,h , such that

Ω w • ŵ - Ω p n+1/2 dae (∇ • ŵ) + Ω (∇ • w)q = - ρ τ Ω (u n+1 m -u n m ) • ŵ -ρ Ω (u n+1/2 m • ∇u n+1/2 m ) • ŵ (12)
for all ŵ ∈ RT 1 h and q ∈ P 1,h , with w • n = ŵ • n = 0 on Γ. This approach has not yet been reported, to the best of the authors knowledge. Note also that here the integration by parts of the convective and viscous terms of the right-hand side, as by the STEint, is for formal reasons not possible since ŵ ∈ H div (Ω) and it can be easily shown that H 1 (Ω) ⊂ H div (Ω).

Integral momentum relative pressure estimator (IMRP)

In this section we derive a new formulation for the integral relative pressure estimator based on a general class of test functions.

Assume that we have computed some specific function v ∈ H 1 (Ω) satisfying

∇ • v ≈ 0 , v • n = 0 on Γ w (13)
If we multiply Equation (1) 1 by v and we integrate over Ω

ρ Ω ∂ t u • v I kin (∂ t u) + ρ Ω (u • ∇u) • v I conv (u) + Ω ∇p • v I pres -µ Ω ∆u • v I visc (u) = 0 . ( 14
)
we obtain for each of the terms

I conv (u) = -ρ Ω (u • ∇v) • u + ρ Γ (u • n)(u • v) (15)
I pres = - Ω p(∇ • v) + Γ p (v • n) = Γ i ∪Γo p (v • n) (16)
I visc (u) = µ Ω ∇u : ∇v -µ Γ (∇u • n) • v (17)
using standard integration by parts. The integral momentum relative pressure estimator (IMRP) is then formulated by assuming that: (i) the measurements u n m satisfy relation ( 14), (ii) the pressure is nearly constant on Γ i and Γ o . Doing so, the IMRP estimator at time t n+1/2 is formulated as:

δp n+1/2 imrp (u n,n+1 m ) = - 1 Λ(v) I kin ( u n+1 m -u n m τ ) + I conv (u n+1/2 m ) + I visc (u n+1/2 m ) . (18)
Note that the requirements on the function in [START_REF] Warnes | ACC/AHA 2008 Guidelines for the Management of Adults with Congenital Heart Disease: a report of the American College of Cardiology/American Heart Association Task Force on Practice Guidelines (writing committee to develop guidelines on the management of adults with congenital heart disease)[END_REF] serve only so that the pressure field vanishes (up to the relative pressure) in ( 16).

Notice also that the integral terms I kin , I conv , I visc are multiplied by 1/Λ(v) in the estimator formulation, which is independent of the measurements. This is an advantage compared to the WERP for the potential spurious amplifications for low flows. Another advantage will be clear after the bias analysis. We also point out that we did not make any assumption on the velocity measurements on the boundary Γ w , as it was done for the WERP. Hence it is potentially better at capturing the pressure gradient changes due to the wall deformation.

We propose the test function v as the solution of the following Poiseuille flow, i.e. a Stokes flow solution with unitary Neumann load on Γ

i : Find v ∈ [P b 1,h ] d and z ∈ P 1,h , Ω ∇v : ∇ψ + Ω z∇ • ψ + Ω q∇ • v + Γ i ψ • n = 0 , (19)
for all ψ ∈ [P b 1,h ] d and q ∈ P 1,h , and with v = ψ = 0 on Γ w . It is not necessary to introduce any coefficient in the viscous term in (19) since the IMRP estimator is independent of any constant scaling of the velocity field, see Equation ( 18). The same applies to the unitary Neumann load in (19). Note also that in the specific choice of this auxiliary function, since v = ψ = 0 on Γ w all terms in that boundary in ( 15)-(17) vanish.

Remark 1. The final relative pressure estimation seems to not considerably depend on the order of the finite element space used for the test function in the numerical experiments. Differences between the different spaces can only be appreciated for coarse mesh spacing of the measurements, i.e. larger than 3 mm. The same applies to the discretization of the STE, STEint and DAE mixed problems.

Bias analysis

We remind the reader that we deal with additive random perturbations of the measured (subsampled) velocity. In this section we detail then how to compute the bias of the estimation methods introduced previously, what has not been done in all original works (i.e for PPE, STE and WERP), to the best of the authors knowledge.

Generalities

The computation of the bias implies in practice to calculate the following quantity:

E δp n+1/2 (u n,n+1 h + ε n,n+1 ) -δp n+1/2 (u n,n+1 h ) ( 20
)
with E the usual expected value operator. Notice that it does not take into account the estimation errors due to the spatial or temporal subsampling. In simple words, bias is how much the mean of the estimators when including noise in the measurements differ from the estimator with noise-free measurements.

Since ε n ∈ [P 1,h ] d is a discrete noise field we can write it (using the Einstein summation convention from now on) as ε n i (x) = N j (x)e n i,j , with N j (x) the finite element shape function for the j-th degree-of-freedom of P 1,h and i denotes the spatial direction.

Remember also that the degrees-of-freedom e n i,j are independently identically distributed (i.i.d.) following a normal distribution N (0, σ 2 ) for all n, i, j. For any linear, differential operator with deterministic coefficients D : H 1 (Ω) → L 2 (Ω) we have therefore the following Identity:

E Ω D(ε n i ) = Ω E e n i,j D(N j (x)) = 0 (21)

PPE

To compute then the bias (20) for the PPE method we proceed as follows. Note that the bias function b ppe := E(p

n+1/2 ppe (u n,n+1 m )) -p n+1/2 ppe (u n,n+1 h ) satisfies the problem: Find b ppe ∈ P 1,h such that Ω ∇b ppe • ∇q = E - ρ τ Ω (u n+1 m -u n m ) • ∇q + ρ τ Ω (u n+1 h -u n h ) • ∇q -E ρ Ω (u n+1/2 m • ∇u n+1/2 m ) • ∇q + ρ Ω (u n+1/2 h • ∇u n+1/2 h ) • ∇q (22)
for all q ∈ P 1,h , and b ppe = q = 0 on Γ o . Since u n m = u n h + ε n , the kinetic terms of the right-hand-side of Problem (22) are zero due to Identity (21). For the convective terms, the only ones surviving are (since the expected value of the cross terms is zero due to Identity (21))

ρ Ω E (ε n+1/2 • ∇ε n+1/2 ) • ∇q = ρ Ω E ε n+1/2 k ∂ k ε n+1/2 i ∂ i q = ρ Ω E N j e n+1/2 k,j ∂ k (N e n+1/2 i, ) ∂ i q = ρ Ω E e n+1/2 k,j e n+1/2 i, N j ∂ k N ∂ i q = ρ Ω σ 2 2 N j ∂ k N j ∂ k q = ρσ 2 4 Ω ∂ k j (N j ) 2 ∂ k q = ρσ 2 4 Ω ∇α • ∇q with α = j (N j ) 2 in P 2,h
. Therefore, we can compute the bias spatial distribution b ppe by solving:

Ω ∇b ppe • ∇q = - ρσ 2 4 Ω ∇α • ∇q . ( 23
)
with b ppe = q = 0 on Γ o . Finally, the bias for the relative pressure δp can be computed by evaluating Formula (2) with p = b ppe . We want to remark that this expression provides a closed form for the bias, depending on known constants and the geometry itself, and it is independent of the real velocity field and hence also of time. Therefore, it can be always evaluated a priori and in the numerical examples we will verify that it is absolutely negligible with respect to the relative pressure of interest. This will be also the case with the rest of the methods, except WERP.

STE and DAE

Analogous to the PPE case, since they all share the same form for the right-hand-side, the bias of STE and DAE can be computed by solving the following problems.

For the STE: Find w ∈ [P b 1,h ] d and b ste ∈ P 1,h , such that

Ω ∇w : ∇ ŵ - Ω b ste (∇ • ŵ) + Ω (∇ • w)q = - ρσ 2 4 Ω ∇α • ŵ (24)
for all ŵ ∈ [P b 1,h ] d and q ∈ P 1,h , with w = ŵ = 0 on Γ. And now for the DAE: Find w ∈ RT 1 h and b dae ∈ P 1,h , such that

Ω w • ŵ - Ω b dae (∇ • ŵ) + Ω (∇ • w)q = - ρσ 2 4 Ω ∇α • ŵ (25)
for all ŵ ∈ RT 1 h and q ∈ P 1,h , with w • n = ŵ • n = 0 on Γ.

STEint

Due to the integration by parts in the right-hand-side, the STEint case has to be handled slightly different. Again, due to Identity (21) the kinetic, and now also the viscous terms in the right-hand-side of the bias equation, vanish. The only surviving term in the convective part has the form:

ρ Ω E (ε n+1/2 • ∇ ŵ) • ε n+1/2 = ρ Ω E ε n+1/2 i ∂ i ŵj ε n+1/2 j = ρ Ω E e n+1/2 i,k e n+1/2 j, N k N ∂ i ŵj = ρσ 2 2 Ω α∇ • ŵ (26)
Therefore, for the STEint the bias equation reads:

Find w ∈ [P b 1,h ] d and b stei ∈ P 1,h , such that Ω ∇w : ∇ ŵ - Ω b stei (∇ • ŵ) + Ω (∇ • w)q = ρσ 2 2 Ω α∇ • ŵ (27)
for all ŵ ∈ [P b 1,h ] d and q ∈ P 1,h , with w = ŵ = 0 on Γ.

IMRP

The bias of the IMRP can be directly defined as

b irmp = E(δp n+1/2 imrp (u n,n+1 m )) -δp n+1/2 imrp (u n,n+1 h )
It can be easily verified that due to Identity (21) the bias expression reduces to convective terms only:

b irmp = - 1 Λ(v) -ρ Ω E ε n+1/2 i ∂ j v i ε n+1/2 j + ρ Γ E(ε n+1/2 i n i ε n+1/2 j v j ) = ρσ 2 2Λ(v) Ω α∇ • v - ρ Λ(v) Γ N k N n i v j E(e n+1/2 i,k e n+1/2 j, ) = ρσ 2 2Λ(v) Ω α∇ • v - ρσ 2 2Λ(v) Γ αv • n .

WERP

We will now end the bias analysis with the WERP estimator. We will need for the analysis the following additional identities:

E ρ Ω |ε n | 2 = σ 2 tr(M Ω ), E µ Ω |∇ε n+1/2 | 2 = σ 2 2 tr(K Ω ) , (28)
with M Ω and K Ω the classical mass and stiffness finite element matrices for the Stokes problem in P 1,h , respectively. Then, let us point out that the following identities hold:

The goal is now to compute:

b werp = E(δp n+1/2 werp (u n,n+1 m )) -δp n+1/2 werp (u n,n+1 h )
In order to perform such an analysis, we include an additional assumption that Λ(u

n+1/2 m ) ≈ Λ(u n+1/2 h
), what is reasonable if we are interested in estimating the peak systolic relative pressure, which typically is simultaneous to the peak flow where the signal-to-noise ratio is best.

We proceed therefore as follows:

E(δp n+1/2 werp (u n,n+1 m )) ≈ -1 Λ(u n+1/2 m ) E E kin (u n+1 m ) -E kin (u n m ) + E conv (u n+1/2 m ) + E visc (u n+1/2 m
) .

(29) We now compute separately each term. First the kinetic part

E (E kin (u n m )) = E ρ 2τ Ω |u n h + ε n | 2 = ρ 2τ Ω E |u n h | 2 + 2u n h • ε n + |ε n | 2 = E kin (u n h ) + σ 2 2τ tr(M Ω ) .
Then, the convective part can be reduced using Identity (21):

E (E conv (u n m )) = ρ 2 Γ i ∪Γo E ((u n h + ε n ) • n)|u n h + ε n | 2 = ρ 2 Γ i ∪Γo E (u n h • n) |u n h | 2 + 2u n h • ε n + |ε n | 2 +(ε n • n) |u n h | 2 + 2u n h • ε n + |ε n | 2 = E conv (u n h ) + ρ 2 Γ i ∪Γo E (u n h • n)|ε n | 2 + 2(ε n • n)(u n h • ε n )
where it can be easily shown that the very last (cubic) term in the second row vanishes due to the standard result that the cube of a normal variable has also zero expected value. Now we continue using the fact that

ε n i = N j e n i,j E (E conv (u n m )) = E conv (u n h ) + ρ 2 Γ i ∪Γo (u n h • n) i E (N j e n i,j ) 2 + 2E ((ε n • n)(u n h • ε n )) = E conv (u n h ) + ρ 2 Γ i ∪Γo (u n h • n)(dσ 2 )α + 2E (N k e n i,k )(N e n j, ) n i u n h,j = E conv (u n h ) + ρ 2 Γ i ∪Γo (u n h • n)(dσ 2 )α + 2σ 2 αu n h • n = E conv (u n h ) + B conv (u n h , σ 2 ) , B conv (u n h , σ 2 ) = ρσ 2 (d + 2) 2 Γ i ∪Γo (u n h • n)α
with d the spatial dimension (in our examples later d = 3). And at last, the viscous part

E E visc (u n+1/2 m ) = E visc (u n+1/2 h ) + µ Ω 2∇u n+1/2 h : ∇ε n+1/2 + |∇ε n+1/2 | 2 = E visc (u n+1/2 h ) + σ 2 2 tr(K Ω ) .
Now, we inserted the computed expected values into Equation (29) and we then obtain the following expression for the bias of the WERP estimator

b werp ≈ -1 Λ(u n+1/2 m ) 1 2 B conv (u n+1/2 h , σ 2 ) + σ 2 2 tr(K Ω ) . (30)
Note that B conv depends on both real (unknown) spatially subsampled velocity u n+1/2 h and noise variance σ 2 , while the second term of the bias depends on σ 2 and the geometry Ω only. However, we will show in the examples that this term is very small B conv compared with the viscous part of the bias.

Remark 2. It is not straightforward to perform similar computations for the variance of the estimators. For the PPE, STE, STEint and DAE, no expression or equation (to our knowledge) can be formulated. For the WERP and IMRP, even though the expressions for the pressure drop are explicit, in the variance the expected value of products of integrals containing the noise is needed. This can not be computed as before since these integrals are not i.i.d.

Numerical examples 6.1 Forward simulations

The synthetic measurements are obtained from a Navier-Stokes simulation in a threedimensional geometry representing a stenotic blood vessel, see Figure 1. The radii of inlet and outlet are 1 cm and the stenosis has a radius of 0.4 cm, i.e. 60 % of coarctation. The physical parameters were chosen as µ = 0.035 Poise and ρ = 1.0 gr/cm 3 . An homogeneous Neumann load on Γ o was imposed, while an homogeneous Dirichlet condition on Γ w and the following Dirichlet condition on Γ i are considered:

u = 60(x 2 + z 2 -1) sin 5π 2 t n on Γ w (31)
where x, z are the coordinates of the plane perpendicular to the main flow direction (the vessel centerline lies on x = z = 0). This input data leads to a peak inflow rate and pressure drop of 94 cm 3 /s and 21 mmHg, respectively. This pressure drop is typical from transitional coarctations from mild to severe and represents a typical threshold for clinical decisions [START_REF] Warnes | ACC/AHA 2008 Guidelines for the Management of Adults with Congenital Heart Disease: a report of the American College of Cardiology/American Heart Association Task Force on Practice Guidelines (writing committee to develop guidelines on the management of adults with congenital heart disease)[END_REF].

The numerical solution of the reference problem was performed using a monolithic velocitypressure coupling. The space-semi-discretization was performed with P b 1,h /P 1,h finite elements for the velocity and pressure, respectively, with a mesh spacing of h = 0.06 cm. A backward Euler scheme was used for the time-semi-discretization with time-step of 0.002 s, with a semi-implicit treatment of the convective term, and including a Temam stabilization term to ensure Lyapunov stability of the solution in time [START_REF] Temam | Une méthode d'approximation de la solution des équations de navierstokes[END_REF]. At that spatial refinement level, no volume or backflow convective instabilities were observed, so no volume convective or backflow stabilization terms were included in order to not perturb the original equation, e.g. SUPG [START_REF] Alexander | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF] or backflow stabilization [START_REF] Bertoglio | A Stokes-residual backflow stabilization method applied to physiological flows[END_REF].

Synthetic measurements

Synthetic measurements (see Figure 2) were generated considering the following realistic perturbations of the data:

• First, the highly resolved reference solution was subsampled on grids with element spacing of 0.1, 0.15, 0.2 and 0.25 cm.

• Then, it was subsampled in time at 0.02 s.

• Finally, an additive Gaussian noise was added to the subsampled velocity at each mesh node and for each time. The standard deviation chosen as 10% of the peak velocity for all velocity components as it is usually done in clinical Phase-Contrast imaging adjusting the expected peak velocity (called VENC parameter) [START_REF] Dyverfeldt | 4D flow cardiovascular magnetic resonance consensus statement[END_REF].

Weighting functions for IMRP

We will use a Poiseuille flow as test functions v for the IMRP method, which can be appreciated in Figure 3. We remark that this is not the only possible choice for v. Alternatives may be Brinkman or Darcy flows. Nevertheless, In several numerical test that we performed no advantage over a simple Poiseuille flow was established in the pressure gradient results. Specifically, the results of the Brinkman flow were dependent on the weighting of the mass term, possibly leading to deteriorated results if this weight was chosen either too large or too small. In the case of a Darcy flow, this is in fact a naive choice due to the lack of regularity of the solutions gradient, which is need in the evaluation of the IMRP in Equation (18). A further option would be to include convective terms for v, but this leads to considerably worse approximation results of the IMRP. Therefore, we do not present any results for these alternative test functions for the sake of clarity and conciseness.

Bias of the estimators

As we remarked in Section 5, we can compute a priori the bias for all methods, which are shown in Table 1. Only in the case of the WERP, the bias coming from the convective term needs to be computed from the ground truth. However, we show that for small mesh spacing the bias of the WERP is dominated by the viscous part, which is again independent on the ground truth velocity field, see Equation (30).

For this example, the bias for the STE, STEint and DAE are of the order of the machine precision. The bias of the PPE is also very small. The IMRP turns to give a negligible bias.

Estimation results: noise-free measurements

We first study the sensitivity of the relative pressure estimation to spatial subsampling of the measurements, without any random noise. The results are shown in Figure 4, and the error convergence curves with respect to the data resolution h are shown in Figure 5. Note that the best methods turn out finally to be the IMRP and STEint, with the former giving a slightly better precision. They turn out to be more robust when subsampling the measurements compared with the state of the art methods, i.e. PPE, STE and WERP. Note that the PPE, which is the oldest and the most spread approach, was consistently the worst method, in line with the findings of [START_REF] Švihlová | Determination of pressure data from velocity data with a view toward its application in cardiovascular mechanics. Part 1. Theoretical considerations[END_REF][START_REF] Donati | Non-invasive pressure difference estimation from PC-MRI using the workenergy equation[END_REF]. We also remark that the temporal subsampling at this realistic level does not have a considerable impact on the estimation precision, compared with the space subsampling, and hence only the time-subsampled results are shown for the sake of conciseness.

Estimation results: noisy measurements

Now we summarize the results with noisy measurements, for which we consider 100 random realizations of the noise ε n for each n. The standard deviations of the relative pressure estimation and the error of the mean with respect to the ground truth (both at peak systole) are presented in Figure 6 and7. Note that both diminish with the spatial discretization size h for all methods. In addition, by comparing Figures 4 and6 that the bias is imperceptible, consistently with the values computed a priori in Table 1. Also consistent with Table 1, the only exception is the WERP.

As expected due to the "squaring" of the measurements, the method with highest variance is the WERP, while the methods delivering sharper standard deviations are PPE and DAE. However, the differences in standard deviations among all methods are around 1.0 mmHg, meaning around 2.0 mmHg for a 95% confidence interval.

In summary, these results confirm that the best performing methods (also when including noise and subsampling) are the IMRP and the STEint.

Conclusions

In this work we reviewed, improved, analyzed (both theoretically and numerically) and compared old and recent approaches for estimating relatives pressures from 4D flow data. Consistent with recent findings, the worst performance was found the be of the widely spread PPE method in terms of its robustness the spatial subsampling.

The WERP method recently introduced in [3], while being the computationally cheapest, it is more sensitive to noise, delivering biased estimations and larger variances than the other methods. The bias of the WERP at higher data resolutions is controlled by the viscous term. However, in our specific numerical test, this term is anyway negligible with respect to the convective one, so it could be neglected from the WERP formulation, avoiding biasing.

We also performed computations for the STE method also recently introduced in [START_REF] Švihlová | Determination of pressure data from velocity data with a view toward its application in cardiovascular mechanics. Part 1. Theoretical considerations[END_REF], confirming that it performs better than the PPE. However, it also presented an important sensitivity when the measurements are subsampled in space. We found out that just a simple integration by parts, in particular of the convective term, considerably reduces this sensitivity, and we called it STEint. This was not shown in the original article [START_REF] Švihlová | Determination of pressure data from velocity data with a view toward its application in cardiovascular mechanics. Part 1. Theoretical considerations[END_REF]. In this context, we investigated an alternative to the STE method based on the Darcy flow, which we called DAE, obtaining slightly better results than the PPE but worse than the rest of the approaches in terms of estimated mean.

Finally, inspired from the integral approach introduced by the WERP, we proposed a new method based on an integral momentum balance weighted by a simple Poiseuille solution. This turns out to be the most accurate method together with the STEint, with a low sensitivity to spatial subsampling in current clinical ranges of voxel size used for 4D flow acquisitions.

Summarizing, based on the numerical experiments presented in this work, the methods of choice would be the IMRP and/or the STEint, both newly introduced in our contribution. While the STEint has the advantage of given the full spatial distribution of the pressure (relative to some point), the IMRP is computationally cheaper when the relative pressure needs to be computed at many time frames. Therefore, future work will consist on further comparing and validating these approaches using real 4D flow data and catheter measurements.

Figure 1 :

 1 Figure 1: Stenosed geometry with 60% coarctation.

Figure 2 :

 2 Figure 2: (a)Reference simulated velocity field (with h = 0.06 cm and mini-elements). (b)-(c) Spatially subsampled measurements at time t = 0.25 s (noise free). (d)-(e) Spatially subsampled and noisy (n) measurements at time t = 0.25.

  (a) h = 0.15 cm (b) h = 0.25 cm

Figure 3 :

 3 Figure 3: Poiseuille test functions for IMRP for different spatial resolutions.

Figure 4 :

 4 Figure 4: Reference ground-truth (GT) and noise-free relative pressure estimations for two different spatial subsampling resolutions including temporal subsampling of 0.02 s.

Figure 5 :

 5 Figure 5: Error of estimation with respect to the ground truth at peak systole for noise-free but spatial subsampled measurements.

Figure 6 :

 6 Figure 6: Mean pressure estimations with noisy measurements.

Figure 7 :

 7 Figure 7: Estimation results at peak systole with noisy measurements.

Table 1 :

 1 Computation of the estimators bias [mmHg] for 10% noise in the velocity.

	h [cm]	0.1	0.15	0.2	0.25
	PPE	5.7e -13	2.1e -13	1.4e -13	1.0e -13
	STE	-1.6e -16 -1.0e -16 -1.6e -16 -1.6e -16
	STEint	1.0e -17	8.6e -17	7.8e -17 -2.4e -17
	DAE	3.3e -17	1.7e -16	7.0e -17 -9.6e -18
	WERP	2.3	0.9	0.6	0.5
	WERP(visc)	2.3	1.0	0.3	1.1
	IMRP	4.2e -06 -5.1e -06 -7.9e -06 -4.1e -08

The relative pressure estimation problemIn this work we deal with the estimation of relative pressures from velocity image data, so we first detail in this section the assumptions about the flow physics and the data.Let us consider an incompressible, Newtonian fluid in a bounded domain Ω ⊂ R d , d =

3, modeled by the incompressible Navier-Stokes equations with the velocity u(t) : Ω → R d and

Acknowledgements

We sincerely thank discussions with Alfonso Caiazzo and Alexander Linke from WIAS Berlin about finite element discretizations, and with Sergio Uribe and Joaquín Mura from the Center for Biomedical Imaging of the Universidad Católica de Chile about Phase-Contrast MRI. D. Nordsletten acknowledges funding from BHF New Horizons program (NH/11/5/29058) and the Engineering and Physical Sciences Research Council (EP/N011554/1).

The rest of the authors sincerely thank the Conicyt Basal Program PFB-03 for its support.