About the entropic structure of detailed balanced multi-species cross-diffusion equations
Esther S. Daus Laurent Desvillettes Helge Dietert


Macroscopic limit]About the entropic structure of detailed balanced multi-species cross-diffusion equations

E. S. Daus]Esther S. Daus
Institute for Analysis and Scientific Computing, Vienna University of
	Technology, Wiedner Hauptstraße 8–10, 1040 Wien, Austria
esther.daus@tuwien.ac.at

L. Desvillettes]Laurent Desvillettes
Université Paris Diderot, Sorbonne Université, CNRS, Institut de Mathématiques de Jussieu-Paris Rive Gauche, IMJ-PRG, F-75013, Paris, France.
desvillettes@math.univ-paris-diderot.fr

H. Dietert]Helge Dietert
Université Paris Diderot, Sorbonne Université, CNRS, Institut de Mathématiques de Jussieu-Paris Rive Gauche, IMJ-PRG, F-75013, Paris, France.
dietert@math.univ-paris-diderot.fr



The first author acknowledges partial support from the
  Austrian Science Fund (FWF), grants P27352 and P30000. The
  research leading to this paper was also partially funded by the
  French “ANR blanche” project Kibord: ANR-13-BS01-0004, by the GDRI
  of the CNRS ReaDiNet, Reaction-Diffusion Network in Biomedecine, and
  by Université Sorbonne Paris Cité, in the framework of the
  “Investissements d'Avenir”, convention ANR-11-IDEX-0005. The last
  author acknowledges partial support from the People Programme (Marie
  Curie Actions) of the European Union’s Seventh Framework Programme
  (FP7/2007-2013) under REA grant agreement n.
  PCOFUND-GA-2013-609102, through the PRESTIGE programme coordinated
  by Campus France.


  This paper links at the formal level the entropy structure of a
  multi-species cross-diffusion system of Shigesada-Kawasaki-Teramoto
  (SKT) type (cf. <cit.>) satisfying the detailed balance
  condition with the entropy structure of a reversible microscopic
  many-particle Markov process on a discretised space. The link is
  established by first performing a mean-field limit to a master
  equation over discretised space. Then the spatial discretisation
  limit is performed in a completely rigorous way. This by itself
  provides a novel strategy for proving global existence of weak
  solutions to a class of cross-diffusion systems.




[2000]35K55, 35K57, 35Q92, 60J28, 82C22, 92D25

[
    [
    December 30, 2023
=====================





§ INTRODUCTION



We consider the population dynamics cross-diffusion system model
coming out of the classical paper by Shigesada, Kawasaki and Teramoto
<cit.> (SKT model) for n≥ 2 species without reaction
term. For clarity, we suppose that the species live on the torus
= [0,1) with periodic boundary conditions. Thus, the density
u_i := u_i(t,x) of species i=1,… ,n evolves as

    ∂_t u_i = Δ(D_i u_i + ∑_j=1^n A_ij u_j u_i )

with diffusion constants D_i ≥ 0, self-diffusion coefficients A_ii>0 and cross-diffusion coefficients A_ij≥ 0 for i,j=1,…,n.

For this system (<ref>), Chen, Daus and Jüngel showed in <cit.> that

    (u) := ∫_∑_i=1^nπ_i [
        u_i log (u_i(x)) - u_i(x) + 1
      ]  x

with positive constants π_i>0 for i=1,… ,n is an entropy (Lyapunov)
functional if the following condition holds

    π_i A_ij = π_j A_ji    for i,j=1,…,n,

which for n ≥ 3 gives a constraint on the cross-diffusion
coefficients A_ij.  Under this condition (they called it
detailed balance condition), the authors were then able to
construct global weak solutions to (<ref>) for an
arbitrary number of population species with the help of the gradient
estimates coming from the entropy production of the entropy
(<ref>).

The motivation of this work is to understand the origin of the entropy
(<ref>) under the condition
(<ref>). In particular, we wanted to link the
condition (<ref>) to the detailed balance
equation of finite-state Markov chains, where the detailed balance
equation has been identified as necessary and sufficient condition for
the existence of a gradient flow structure with respect to the
relative entropy
<cit.>.

In this work, we establish the formal link between the entropy
structure of (<ref>) and the entropy structure of a
microscopic many-particle Markov process on a discrete space. The link
is established in two steps. In the first step, we perform a formal
mean-field limit keeping the spatial discretisation fixed. The
resulting system is a quadratic master equation. In the second step,
we then refine the spatial discretisation and arrive at the cross-diffusion system
(<ref>).

In this many-particle derivation, the condition
(<ref>) enters as a natural necessary condition
for the construction of a reversible Markov process and the
constants π_i can be interpreted as relative portions in the
many-particle model.

In both limits, the entropy structure is preserved and, in particular,
the master equation on the discretised space has the corresponding
entropy structure. This allows us to perform the spatial
discretisation limit in a rigorous way, which is an interesting result
by itself (which will be discussed after the statement of our main
Theorem <ref> in <ref>).

Note that the transfer of the entropy structure from a microscopic model towards a mesoscopic model has been extensively studied for equations belonging to other classes. The spatially homogeneous Boltzmann equation is for example a model in which many results have been proven (cf. <cit.>). It shares some features with the SKT model (quadraticity of course, but also diffusive properties when the angular cutoff of Grad is not performed).

For the second rigorous limit from the space discretised master
equation to the SKT cross-diffusion system, similar discrete in space
approximation schemes for the SKT model were studied in <cit.>, but they are not necessarily entropy preserving. Very
recently, an entropy preserving numerical scheme was proposed in
<cit.>, though not for the SKT model, but for a volume-filling
type cross-diffusion system.

Other approaches have been proposed for obtaining cross-diffusion
equations of SKT type out of microscopic models. First (stochastic) approaches from particle models to reaction-diffusion systems trace back to Oelschläger <cit.> in the late 1980s. Recently, Fontbona and Méléard <cit.> managed to prove the convergence from realistic individual-based models in a suitable
limit towards non-local (convoluted w.r.t. space) SKT-type
systems. Note that because of the lacking evidence of existence and
uniqueness of strong solutions to the limiting model, it looks
difficult to provide a rigorous proof of passage to the limit towards
the full (i.e. nontriangular) local multi-species SKT system when one
starts with a microscopic model (whether on a discrete set of
positions or on a continuous set of positions, using a nonlocality
which disappears in the limit). Note, however, that very recently, Moussa [21] manged to prove the convergence in the
  case of strictly triangular limiting local SKT model with bounded
  coefficients starting on a continuous set of positions with a
  nonlocality which disappears in the limit by using duality
  techniques (introduced for instance by Pierre and Schmitt in <cit.>).
  
  
  
  
  
  




§ FORMAL MEAN-FIELD LIMIT

For the microscopic derivation, we first consider a many-particle
system on a fixed spatial discretisation. The spatial discretisation
consists of M positions given by

    Ω_M ={x_k : k=0,…,M-1}    with    
      x_k = k/M = kh,

which is understood in the periodic setting, and where  we set h=M^-1.

Given the relative fractions π_1,…,π_n with
  π_i >0 between the species, we consider the
many-particle system with π_i N particles
of species i=1,…,n, where π_i N
denotes the largest integer smaller than π_i N. The aim of this
section is to obtain a suitable master equation when N→∞.

The microscopic configuration is given by

    x := (x_1^1,…,x_1^π_1 N,x_2^1,…,x_2^π_2
        N,……,x_n^1,…,x_n^π_n N)
      ∈Ω_M^⊗(π_1N+…+π_nN) =: Ω_M^N

and this configuration is set to evolve in time as a time-continuous
Markov chain.

The distribution over the microscopic configurations at time t is
given by a density μ^N_t ∈(Ω^N_M). In terms of
statistical physics, this means that we consider an ensemble over the
microscopic configurations.

We assume that the particles within a species are
indistinguishable. The class of such measures is denoted by
_s(Ω_M^N) and defined as follows:


  A measure μ∈(Ω^N_M) is in _s(Ω^N_M)
  if and only if for all permutations σ_1,…,σ_n of
  resp. {1,..,π_1N}, …, {1,..,π_nN} and configurations
  x∈Ω_M^N it holds that
  
    μ^N(x_1^σ_1(1),…,x_1^σ_1(π_1 N),
          x_2^σ_2(1),…,x_2^σ_2(π_2 N),
          ……, x_n^σ_n(1),…,x_n^σ_n(π_n N)) 
       = μ^N(x_1^1,…,x_1^π_1 N,x_2^1,…,x_2^π_2
            N,……,x_n^1,…,x_n^π_n N).



By the indistinguishability, the distribution of a typical particle is
given by the marginal distribution. For this, we first introduce the
following notation for projections.


  Let p=(p_1,…,p_n) and N be such that
  p_i ≤π_i N for i=1,…,n.  We define the
  projection
  
    ^N;(p) : (Ω_M^N)
        ↦(Ω_M^⊗ (p_1+… +p_n))

  by
  
    (^N;(p)μ^N)(x)
          := ∑_x_1^p_1+1∈Ω_M…∑_x_1^π_1N∈Ω_M   ∑_x_2^p_2+1∈Ω_M…∑_x_2^π_2N∈Ω_M……∑_x_n^p_n + 1∈Ω_M…∑_x_n^π_nN∈Ω_M
       μ^N(x_1^1,…,x_1^π_1 N,x_2^1,…,x_2^π_2
            N,……,x_n^1,…,x_n^π_n N)

  for
  
    x := (x_1^1,…,x_1^p_1,x_2^1,…,x_2^p_2,……,x_n^1,…, x_n^p_n).

  For μ^N ∈(Ω_M^N), we denote the marginal by
  
    μ^N;(p) = ^N;(p)μ^N∈(Ω_M^⊗ (p_1+…+p_n)).



We then expect to recover the master equation from the first marginals

    u_i := μ^N;(e_i),

where e_i = (0,…,0,1,0,…,0) denotes the unit vector with n
components, where the 1 is at the i-th component.

The quadratic terms are expected to originate from a binary
interaction in the particle model. The diffusion is the result of a
random walk and the nonlinearity given in the SKT model
(<ref>) is expected to come out from jumps of

particles interacting at the same position.

The entropy structure is expected to be linked to the reversibility of
the Markov chain. We therefore introduce a binary interaction which
happens in a reversible way. This can be realised by imposing the same
jump for the interacting particles.

This leads us to consider the following class of particle models,
describing the evolution of the microscopic configuration.


  
  Let D_i and D_ij be nonnegative constants such that D_ij = D_ji for
  i,j=1,…,n. For a fixed N, define the time-continuous Markov
  chain on Ω_M^N by the transitions
  
    .
          x→x + e_i^a + e_j^b 
    x→x - e_i^a - e_j^b
          }    with rate    δ_(i,a)≠(j,b)δ_x_i^a=x_j^bD_ij/N
       .
          x→x + e_i^a 
    x→x - e_i^a
          }    with rate    
                             D_i

  for i,j=1,…,n and a=1,…,π_iN, b=1,…,π_jN,
  where e_i^a is the vector with components of value zero at all places, except for the ath particle of species
  i, where the value is h = 1/M. The Markov chain is defined to
  have no other transitions.



  The transition rates are well-defined if and only if
  D_ij = D_ji, which will lead to condition
  (<ref>).


From the construction, we directly see the reversibility.

  The Markov chain given in <ref> is
  reversible and the stationary distribution is the homogeneous
  distribution, where each x∈Ω^N_M has the same
  probability |Ω^N_M|^-1 = M^-(π_1 N + …+ π_n N).


  This can be obtained by a direct computation.


We now suppose that the microscopic configuration evolves according to
the Markov chain. Then the distribution μ^N solves the following linear ODE:

    / tμ^N(x)
           = ∑_i=1^n∑_a=1^π_iN
        D_i
        [μ^N(x+e_i^a)+μ^N(x-e_i^a)-2μ^N(x)] 
       + 1/2∑_i=1^n ∑_a=1^π_iN∑_j=1^n ∑_b=1^π_jNδ_(i,a)≠(j,b)δ_x_i^a=x_j^bD_ij/N[μ^N(x+e_i^a+ e_j^b)+μ^N(x- e_i^a-e_j^b)-2μ^N(x)].

Here we used that x_i^a=x_j^b holds after the pairwise interaction
if and only if it holds before, so that we can factor it out
(that is, δ_x_i^a=x_j^b =
δ_(x+e_i^a)_i^a=(x+e_j^b)_j^b =
δ_(x-e_i^a)_i^a=(x-e_j^b)_j^b).

We further suppose that the particles are indistinguishable, which is
propagated in time.

  Suppose that μ^N is the distribution for the Markov chain given
  in <ref>. If μ^N ∈_s(Ω_M^N) initially holds, then it also holds at all later times.


  It follows directly from the definition of the transition rates, which
  respect the indistinguishability.


 We now write explicitly the formula emphasizing the entropy structure of our reversible Markov process, we recall that due to <cit.>, the time-reversible many-particle continuous time Markov chain is a gradient flow of the relative entropy with respect to its stationary distribution.

  We assume that D_i≥ 0, and
   D_ij=D_ji≥ 0 for all 1≤ i,j≤ n. We also assume that μ^N is initially strictly positive.
   Then, the entropy functional defined by
   
    ℋ̃(μ^N) : =∑_xμ^N(x)log(μ^N(x)/M^(π_1 N+ ⋯ + π_n N))

   is decreasing with respect to time, i.e.
   
    / tℋ̃(μ^N) ≤ 0   for all t>0

   along the flow of (<ref>).
 

  Note first that the strict positivity of μ^N is maintained in
  the evolution of the process, so that the logarithm of μ^N is
  always well defined.

  The proof works in a totally analogous way as the proof of
  (<ref>), where the entropy decay is shown on the macroscopic
  level. For completeness, we sketch the proof also here, by using the
  following notation for any function f: Ω^N_M → (0,∞):
  
    Δ_e_i^a(f(x))   := f(x + e_i^a) + f(x - e_i^a) - 2f(x), 
    Δ_(e_i^a + e_j^b)(f(x))   := f(x + e_i^a + e_j^b) + f(x - e_i^a - e_j^b) - 2f(x), 
    ∇^+_(e_i^a)(f(x))   := f(x + e_i^a) - f(x), 
    ∇^+_(e_i^a + e_j^b)(f(x))   := f(x + e_i^a + e_j^b)- f(x).

  Here e_i^a is the vector with components of value zero at
  all places, except for the ath particle of species i,
  where the value is h = 1/M.  This coincides with the
  notation of a discrete Laplacian and discrete gradient up to
  positive scaling constants. We will introduce them more rigorously
  on the level of the master equation in Section
  <ref>.  Thanks to the periodicity of
  the domain and to a discrete integration by parts (see detailed
  formulas at the beginning of Section <ref>),
  it holds that
  
    / tℋ̃(μ^N)    = ∑_x(logμ^N(x) +1 )d/dt(μ^N(x)) 
        = ∑_x∑_i=1^n ∑_a=1^π_i ND_i (logμ^N(x) +1 ) Δ_e_i^a(μ^N(x)) 
          +1/2∑_x∑_i,j=1^n ∑_a=1^π_i N∑_b=1^π_j Nδ_(i,a)≠(j,b)δ_x_i^a=x_j^bD_ij/N(logμ^N(x) +1 )Δ_(e_i^a + e_j^b)(μ^N(x)) 
       = -∑_x∑_i=1^n ∑_a=1^π_i ND_i ∇^+_(e_i^a)(logμ^N(x))·∇^+_(e_i^a)(μ^N(x)) 
          -1/2∑_x∑_i,j=1^n ∑_a=1^π_i N∑_b=1^π_j Nδ_(i,a)≠(j,b)δ_x_i^a=x_j^bD_ij/N∇^+_(e_i^a + e_j^b)(logμ^N(x))·∇^+_(e_i^a + e_j^b)(μ^N(x)) 
       ≤ 0,

   thanks to the monotonicity of x ↦log x.



The evolution of the marginals is given by the BBGKY hierarchy.

  
  Suppose that μ^N ∈_s(Ω_M^N) is the density
  evolving according to (<ref>). Then the
  marginals evolve as
  
    / tμ^N;(p)(x)
        = I + II + III,

  where
  
    I    =  D_i∑_i=1^n∑_a=1^p_i
            [μ^N;(p)(x+e_i^a)+μ^N;(p)(x-e_i^a)-2μ^N;(p)(x)],
        
    
        II    = 1/2∑_i=1^n ∑_a=1^p_i∑_j=1^n ∑_b=1^p_jδ_(i,a)≠(j,b)δ_x_i^a=x_j^bD_ij/N
             [μ^N;(p)(x+e_i^a+e_j^b)+μ^N;(p)(x-e_i^a-e_j^b)-2μ^N;(p)(x)],

  with e_i^a defined as the vector of size p_1 + … + p_n with all coordinates with value
  0, except the coordinate of index
  p_1 + … + p_i-1 +a, which value is 1/M=h. Finally,
  
    III
             =
            ∑_i=1^n ∑_a=1^p_i∑_j=1^n
            ∑_x_j^p_j+1∈Ω_Mδ_x_i^a=x_j^p_j+1 D_ijπ_jN-p_j/N
         [μ^N;(p+e_j)((x# x_j^p_j+1)+ẽ_i^a+ẽ_j^p_j + 1)
            +μ^N;(p+e_j)((x# x_j^p_j+1)-ẽ_i^a-ẽ_j^p_j +1)
            -2μ^N;(p+e_j)((x# x_j^p_j+1))]

  with
  
    (x# x_j^p_j+1)
        = (x_1^1,…,x_1^p_1,……,x_j^1,…,x_j^p_j,x_j^p_j+1,x_j+1^1,…, x_j+1^p_j+1,……,x_n^1,…,x_n^p_n),

  i.e. x with x_j^p_j+1 added between x_j^p_j and x_j+1^1, and where ẽ_i^a
is defined  as the vector of size p_1 + … + (p_j +1) + … + p_n with all coordinates with value 0, except the coordinate of index p_1 + … + p_j +1,
  which value  is 1/M=h.

The term I is the standard linear diffusion. The term II is the
quadratic interaction between the considered particles, which should
be negligible as N →∞. The term III is the interaction
between the considered particles and the averaged particles, which
leads to the quadratic term. The interaction between the averaged
particles does not appear in the projection.


  Take the projection ^N;(p). The terms I and II follow
  directly. The third term appears as
  
    III =
        ∑_i=1^n ∑_a=1^p_i∑_j=1^n ∑_b=p_j+1^π_jNδ_x_i^a=x_j^bD_ij/N^N;(p)[μ^N(x+e_i^a+e_j^b)+μ^N(x-e_i^a-e_j^b)-2μ^N(x)],

  where we ordered the pair (i,a) and (j,b) so that the factor
  1/2 is not appearing there. Thanks to the indistinguishability, this
  takes the claimed form.


Thus, as usually in the BBGKY hierarchy (tracing back to <cit.>),  in order to compute the evolution of the one-particle marginals
u_i, we need the knowledge of the two-particle marginals, whose
evolution in turn requires the three-particle marginals.

In order to close an equation on u_i, we thus need an additional
assumption. This assumption has been identified as chaos by Kac <cit.>, in
a mathematical setting following the famous Stoßzahlansatz by
Boltzmann (first suggested by J. Clerk Maxwell in  <cit.>).
It states that in the limit N →∞, the
different particles are becoming independent. In terms of the measure
μ^N, it means that

    μ^N(x_1^1,…,x_1^π_1 N,x_2^1,…,x_2^π_2
          N,……,x_n^1,…,x_n^π_n N) 
       ≈ u_1(x_1^1)… u_1(x_1^π_1 N)   u_2(x_2^1)… u_2(x_2^π_2 N)…… u_n(x_n^1)…
        u_n(x_n^π_nN),

as N→∞.  The formal idea is that the interaction
between two particles is scaled as N^-1 so that the correlation
between two particles should also be scaled as N^-1. Therefore, as
N→∞, the particles become independent in the limit. In a way
for a rigorous mathematical treatment, Kac suggested to initially
assume the factorisation and then prove that this is preserved in time
with an error going to zero as N→∞ (propagation of chaos).

 
  Formally, as N →∞ we have under the Stoßzahlansatz
  that the marginals u_i evolve as
  
    / t u_i(x)    = D_i [u_i (x+h) + u_i (x-h) - 2 u_i (x)]
          +∑_j=1^n D_i jπ_j
        [u_j(x+h)u_i (x+h) + u_j(x-h)u_i (x-h) - 2 u_j(x)u_i (x)].



  This follows from <ref>, where in the term III, we see
  that  π_jN-p_j/N→π_j, and we reduce the two-marginal
  density as a product by (<ref>).


This gives the desired quadratic master equation with the final rate
A_ij := D_ijπ_j. This is equivalent to the detailed balance
equation (<ref>), which follows from the
symmetry D_ij=D_ji in the following way:

    π_i A_ij = π_i D_ijπ_j = π_j D_jiπ_i = π_j A_ji    1≤ i,j ≤ n.


 Thanks to the chaos assumption, we can relate the relative entropy of μ^N to the relative entropy of the
 u_i using the following proposition.

  Assume
  
    μ^N(x_1^1,…,x_1^π_1 N,x_2^1,…,x_2^π_2
            N,……,x_n^1,…,x_n^π_n N) 
       = u_1(x_1^1)… u_1(x_1^π_1 N)   u_2(x_2^1)… u_2(x_2^π_2 N)…… u_n(x_n^1)…
          u_n(x_n^π_nN),

  then
  
    1/Nℋ̃(μ^N) = ∑_xμ^N(x)
        log( μ^N(x)/M^(π_1 N+…+π_nN))
        = 1/N ∑_i=1^n ∑_ℓ=0^M-1π_i N u_i(x_ℓ)
        log(u_i(x_ℓ)/M)



  This follows from expanding the logarithm as product and using that
 the  u_i are probability distributions.


When N →∞, this last quantity converges towards ∑_i=1^n π_i  ∑_ℓ=0^M-1 u_i(x_ℓ) log(u_i(x_ℓ)/M).
 Because of Lemma <ref> and Proposition <ref>, we expect that this quantity decreases along the flow of eq. (<ref>). We shall indeed prove this in the next section, thus establishing the link between the
entropy structure for eq. (<ref>) (and its limit when the discretisation step h tends to 0) and the  classical
relative entropy of Markov chains.




§ RIGOROUS DERIVATION TO THE CROSS-DIFFUSION SYSTEM


Starting from the Markov chain defined in
Definition <ref>, we showed in  the last section how performing the mean-field limit on the formal
level leads to the spatial discretisation
(<ref>) of the SKT system (<ref>).  In fact, we
expect this particular discretisation to preserve the entropy
structure of the Markov chain. In this section we shall check this
property and use it to pass rigorously to the limit when the
discretisation step h tends to 0, thus recovering the existence of
weak solutions for the SKT model.

For this, we recall the discretisation
Ω_M = {x_k = kh : k=0,…,M-1} with h=M^-1 from
(<ref>). Moreover, we introduce the discrete derivatives and
discrete Laplacian by

    (∇^+_h f) (x)     := f(x+h) - f(x)/h,
                            (∇^-_h f) (x)  := f(x) - f(x-h)/h, 
    
      (Δ_h f)(x)     := [∇^-_h(∇^+_h f)](x)
                        = [∇^+_h(∇^-_h f)](x)
                        = f(x+h) + f(x-h) - 2 f(x)/h^2.

We now rewrite (<ref>) together with its initial boundary
conditions (and A_ij := π_j  D_ij), after a suitable rescaling in time (such that _t is replaced by h^2 _t). This yields


    {_t u_i(t,x_k) = D_i [Δ_h u_i(t,·)](x_k) + [Δ_h (u_i(t,·)∑_j=1^n A_iju_j(t,·))] (x_k) ,    k=0,…,M-1,
    
        u_i(0,x_k) =u_i^0(x_k)≥ 0,      k=0,…,M-1,    i=1,…,n,
    
        u_i(t,x_0) = u_i(t,x_M),     u_i(t,x_-1) = u_i(t,x_M-1),     ∀ t,    i=1,…,n.
      .


Given the values (w(x_k))_k=0,…,M-1 over Ω_M, let
w̃ : ↦ be the linear interpolant (P_1 discretisation), which can be
defined by

    w̃(x) := ∑_k=0^M-1 w(x_k)  T(x - x_k) + w(x_0)  T(x - x_M),

where

    T(x) = (1 - |x|/h) _[ |x|≤ h ].

With this we can state our main theorem.

  
  Let D_i ≥ 0 and A_ij≥ 0 be coefficients satisfying
  
  
  
  * A_ii > 0 (strict positivity of self-diffusion),
  
  * π_i A_ij = π_j A_ji for some constants π_i >0
    (detailed balance equation).
  
  We also assume continuous positive initial data u_i^0:= u_i^0(x) > 0 on
  𝕋 for all 1 ≤ i ≤ n.

  Then for all M ∈ℕ, there exists a unique global
  solution u_i := u_i(t,x_k) > 0 of class C^∞ to the
  discrete system (<ref>) with h := 1/M.

  Denoting [ũ_i]^M the interpolant obtained from
  u_i(t,x_k) by formula (<ref>), then there exists a
  subsequence such that the following holds:
  [ũ_i]^M →_M →∞ u_i in
  L^4 -([0,T] ×𝕋) for all T>0 and >0,
  where u_i ∈ L^4([0,T]×) ∩ L^2([0,T], H^1()) is a
  weak solution to the SKT system
  _t u_i = Δ (D_i   u_i + ∑_j=1^n A_ij u_i u_j)
  (with initial data u_i^0 and periodic boundary conditions), in the
  following sense: For all
  φ∈ C^2_c([0,∞) ×𝕋) and for all
  i=1,…,n,
  
    - ∫_𝕋 u_i(0,x)  φ(0,x)   x -
          ∫_0^∞∫_𝕋 u_i(t,x) _t φ(t,x)   x  t 
       = ∫_0^∞∫_𝕋[ D_i u_i(t,x) + ∑_j=1^n
          A_ij u_i(t,x)  u_j(t,x)] Δφ(t,x)   x  t.



The existence of weak solutions to the SKT model has been known for a
long time in the case of two equations, and it has been studied more
recently in the case of more than two species, under the detailed
balance condition (cf. <cit.>). We do not go
beyond the existing theory of existence in this paper. We, however,
present an approximation procedure which is extremely simple (it
consists only in discretizing w.r.t. the space variable)
compared to most previous procedures (cf. for
example <cit.>).

Though this approximation procedure is presented here only in the
specific case of the quadratic SKT model without reaction terms under
the assumption of detailed balance in dimension 1 and in the
presence of self-diffusion, our feeling is that it can be easily
extended to more general cases. First, one can introduce (not too quickly increasing) reaction terms. Second, one can go to higher space dimensions d≥ 1 keeping periodic boundary conditions. In a third step, by introducing a reasonable grid, one can expect that the same procedure  works for any reasonably smooth domain with Neumann boundary conditions.
Moreover, one can also think of non quadratic cases, provided that a good Lyapunov
functional is known, or of the quadratic case without self-diffusion when the standard
diffusion term or the reaction terms are sufficient to guarantee the
equintegrability. It is less clear if duality arguments (cf. for
example <cit.>) are compatible with this approximation (as they
are when time discretisation is performed, cf. <cit.>): this
issue will be investigated in future works.

The possibility of extending the formal results of the first part to
more general systems (thus giving a microscopic background for an
entropy structure which is known to exist at the macroscopic level)
will also be studied further, especially in the direction of non
quadratic systems, and systems presenting exclusion processes, see for
instance <cit.>. 

We now begin the


Proof of Theorem <ref>:

We first observe that there exists a unique global solution
t∈_+ ↦ (u_1(t, x_k),…,u_n(t,x_k)) with
u_i(t,x_k)≥ 0 to the ODE system (<ref>). We briefly sketch the proof of this result, which uses standard theorems for ODEs.
We denote by T_1>0 the maximal time of existence for the equation (obtained thanks to Cauchy-Lipschitz theorem), and by T_2 ∈ ]0, T_1[ the maximal time for which  u_i(t,x_k)>0 for all i,k and t∈ [0,T_2[. Note that T_2>0 because all initial data are assumed to be strictly positive.
On the interval [0,T_2[, we use the conservation of
the total number of individuals (of each species)

    / t( ∑_k=0^M-1 u_i(t, x_k)) = 0,

for i=1,…,n, and get that

    ∀ t∈ ]0,T_2[,   i=1,..,n,   k=0,..,M-1,      0 ≤  u_i(t, x_k) ≤ C,

where C : = sup_i=1,..,n[ ∑_k=0^M-1 u_i(0, x_k) ].
Then, we observe that on the interval [0,T_2[,


    h^2  ∂_t u_i(t,x_k)    =  D_i [u_i(t,x_k+h) + u_i(t, x_k-h) - 2 u_i(t, x_k)] 
          +  ∑_j=1^n A_ij[(u_i u_j)(t,x_k+h) + (u_i u_j)(t, x_k-h) - 2 (u_iu_j)(t, x_k)] 
       ≥ -2 D_i  u_i(t,x_k) -2 ∑_j=1^n A_ij(u_iu_j)(t,x_k) 
       ≥ -2 (D_i + C  ∑_j=1^n A_ij)  u_i(t,x_k),

and consequently

    u_i(t,x_k) ≥ u_i(0,x_k)exp(- 2/h^2 [D_i + C  ∑_j=1^n A_ij]  T_2 )>0     t ∈ [0,T_2].

Then T_2=T_1, and finally thanks to estimate eq. (<ref>), T_1=∞.



 §.§ Discrete system


We start by studying the discrete system and establishing the main
a priori estimates, which are uniform with respect to the
spatial discretisation M. Those estimates are a direct consequence
of the entropy structure of our models.

  Under the same assumptions
  on the coefficients and initial data as in <ref>, the
  unique solution to the system (<ref>) satisfies the following a
  priori estimates, for some constants C_T>0 depending only on T,
  the initial data, and the coefficients π_i, A_ij and D_i:
  
    ∑_i=1^n∫_0^T     h ∑_k=0^M-1 |(∇^+_h u_i) (t,x_k)|^2
         t ≤ C_T,

  and
  
    / t^h(u(t,·)) ≤ 0,     sup_t ∈ [0,T]^h(u(t,·))
        ≤ C_T,

  where
  
    ^h(u) := ∑_i=1^n h ∑_k=0^M-1π_i  [u_i(x_k) log(u_i(x_k)) - u_i(x_k) +1 ] .



 The normalised entropy described at the end of Section 2
 
 differs from ^h in (<ref>) only by terms which are constant with respect to time t (due to mass conservation), thus the entropy dissipation of those terms is 0.



For the proof, we rely on the following elementary properties for the
discrete derivatives:



  * Discrete integration by parts: For all 1-periodic
  functions p,q: →,
 
    ∑_k=0^M-1 (∇^+_h p) (x_k)  q(x_k) =
       - ∑_k=0^M-1  p(x_k)  (∇^-_h q)(x_k) ;


  * Discrete product rule: For all functions p,q: →,
  
    (∇^+_h(pq))(x) = p(x+h)  (∇^+_h q)(x) +
        (∇^+_hp)(x)  q(x).




  By using the abbreviation ã_ij:=π_i A_ij (and
  therefore assuming that ã_ij=ã_ji for all
  i,j ∈{1,…,n}), we can compute
  
    / t^h(u) = ∑_i=1^n ∑_k=0^M-1 h π_i _t u_i(x_k) log u_i(x_k) 
       =∑_i=1^n ∑_k=0^M-1 hπ_i (Δ_h[D_i u_i + u_i∑_j=1^n A_i,ju_j])(x_k)log u_i(x_k) 
       = h ∑_i=1^n ∑_k=0^M-1  D_iπ_i log u_i(x_k)  (Δ_h u_i)(x_k)  + h∑_i,j=1^n ∑_k=0^M-1ã_ij log u_i(x_k)  (Δ_h (u_i u_j))(x_k)
       =  -h∑_i=1^n ∑_k=0^M-1 D_iπ_i ( ∇^+_h (log u_i))(x_k)   (∇^+_h(u_i))(x_k) -h∑_i,j=1^n ∑_k=0^M-1ã_ij( ∇^+_h (log u_i))(x_k)   (∇^+_h(u_i u_j))(x_k)  
        =  -h∑_i=1^n ∑_k=0^M-1 D_iπ_i ( ∇^+_h (log u_i))(x_k)   (∇^+_h(u_i))(x_k) - h/2 ∑_i,j=1^n ∑_k=0^M-1ã_ij( ∇^+_h (log (u_i u_j)))(x_k)   (∇^+_h(u_i u_j))(x_k)
       ≤ -4h∑_i=1^n ∑_k=0^M-1 D_iπ_i  |∇^+_h(√(u_i))(x_k)|^2   - 2h ∑_i,j=1^n ∑_k=0^M-1ã_ij | ∇^+_h (√(u_i u_j)))(x_k)|^2  ≤ 0 .

  We used above the elementary inequality
  (x-y)(log x- log y)≥ 4 (√(x)-√(y))^2 for all
  x>0,y>0.

  We end up the proof of estimate (<ref>) by noticing that all
  terms above are nonpositive and by integrating between 0 and
  T. Estimate (<ref>) is obtained by using only the self
  diffusion terms (that is, the ones corresponding to
  ã_ij, for i=j) and also by integrating between 0 and
  T.


Next, we introduce for 1≤ p<∞ the discrete norm for
(w(x_k))_k=0,…,M-1 by defining

    w^p_h,p := h ∑_k=0^M-1 |w(x_k)|^p.


With the help of the following lemma, we can switch between the
discrete norm and the norm of the continuous linear interpolant
w̃ of w defined in (<ref>):


 For 1≤ p < ∞ and w(x_k)≥ 0, k=0,..,M-1, it holds that
 
    w̃^p_L^p()≤w^p_h,p≤p+1/2w̃^p_L^p(),

    
    ∇w̃^p_L^p() = ∇^+_h w_h,p^p,

 where w_h,p is the discrete norm defined in
 (<ref>), and w̃ is the linear interpolant defined
 in (<ref>).


  The factor 2/(p+1) is necessary, as  can be seen from the case when
  w(x_k)=1 for k=1 and w(x_k)=0 otherwise.




  Note that the linear interpolant w̃ can also be rewritten as
  
    w̃(x)=∑_k=0^M-1(α_k(x)w(x_k) + (1-α_k(x))w(x_k+1))_[x_k,x_k+1)(x),

  where w(x_M) := w(x_0) and
  
    α_k(x)=x_k+1-x/h.


  For x ∈ [x_k,x_k+1) with k=0,…,M-1, we know thanks to
  (<ref>) that
  
    w̃(x)=α_k(x)w(x_k)+(1-α_k(x))w(x_k+1).

  Since x ↦ x^p is convex, we see that
  
    |w̃(x)|^p ≤α_k(x)  |w(x_k)|^p  + (1-α_k(x)) |w(x_k+1)|^p,

  so that integrating between x_k and x_k+1, we get
  
    ∫_x_k^x_k+1|w̃(x)|^p  x ≤h/2|w(x_k)|^p + h/2|w(x_k+1)|^p,

  which shows the first part of (<ref>).

  In the other direction, we find that
  
    ∫_x_k^x_k+1|w̃(x)|^p  x    = h ∫_β=0^1 |β w(x_k) + (1-β)w(x_k+1)|^p β
       = h/p+1[w(x_k+1)]^p+1-[w(x_k)]^p+1/w(x_k+1) - w(x_k)
       ≥h/p+1[[w(x_k+1)]^p + [w(x_k)]^p],

  where we used the elementary inequality
  
    A^p+1-B^p+1/A-B≥ A^p + B^p     A,B ≥ 0.

  This elementary inequality is easily proved (by considering
  A/B). This finishes the proof of (<ref>).

  For (<ref>), we see that for k=0,…,M-1,
  
    ∇w̃(x) = (∇^+_h w)(x_k)  for x ∈ (x_k,x_k+1).

  This implies
  
    ∇w̃^p_L^p()
        = ∫_𝕋|∇w̃(x)|^p  x
        =∑_k=0^M-1∫_x_k^x_k+1|∇w̃(y)|^p  y
        = ∑_k=0^M-1 h |(∇^+_h w)(x_k)|^p = ∇^+_h w_h,p^p.





 §.§ Uniform a priori estimates for the linear interpolant


From now on, when we interpolate functions which depend on t, we
systematically write w̃(t,x) instead of
w(t,·)(x). We also use the notation C_T for any
constant depending on the time T, on the inital data and the
parameters π, A_i,j and D_i of the problem, but not on the
discretisation parameter h= 1/M.

Combining <ref> and <ref>, we obtain
for ũ_i with i=1,…,n that

    sup_t ∈ [0,T]∫_𝕋ũ_i(t,x)  x ≤ C_T, 
    ∫_0^T ∫_𝕋|∇ũ_i(t,x)|^2  x  t ≤ C_T.

Using the Gagliardo-Nirenberg inequality, this implies for
p ∈ [1,4] that

    ∫_0^T ∫_𝕋|ũ_i(t,x)|^p   x  t ≤ C_T.

Note that in the estimate above, when the dimension d=1 is replaced by a more
general dimension d, the maximal value 4 of p is replaced by
2+2/d.

Indeed, the Gagliardo-Nirenberg interpolation allows to estimate
ũ_i_L^p() by
∇ũ_i_L^2()^θũ_i_L^1()^1-θ. Choosing θ p=2, means
that θ=(2d(p-1))/((d+2)p) ∈ [0,1] so that p=2+ 2/d. With
this choice, we find

    ũ_i^p_L^p([0,T], L^p())   = ∫_0^T ũ_i^p_L^p()  t 
       ≤ C ∫_0^T ∇ũ_i^θ p_L^2()ũ_i^(1-θ)p_L^1()  t 
       ≤ C ũ_i^(1-θ)p_L^∞([0,T],
        L^1())∫_0^T ∇ũ_i^θ p_L^2()  t 
       ≤ C ũ_i_L^∞([0,T], L^1())^(1-θ)p∇ũ_i^2_L^2([0,T] ×).


Using <ref>, we can relate the estimate back to the discrete
system as

    ∫_0^T h ∑_k=0^M-1|u_i(t,x_k)|^p  t ≤ C_T,
        for all p∈ [1,4].


We now show that for all ϕ∈ W^1,∞(),

    ∫_0^T | _t ∫_ũ_̃ĩ(t,x) ϕ(x)  x
      |   t ≤ C_T ϕ_W^1,∞().


Indeed, performing a discrete integration by parts in x_k
(cf. (<ref>)) and performing the translation
x ↦ x+h inside the integral over , we get that

    _t ∫_ũ_̃ĩ(t,x) ϕ(x)  x
       = ∑_k=0^M-1[Δ_h (D_i  u_i(t,·) +
        u_i(t,·)∑_j=1^n A_ij u_j(t,·))] (x_k)  ∫_  T(x - x_k)  ϕ(x)  x 
       =  - ∑_k=0^M-1[∇^+_h (D_i  u_i(t,·) +
        ∑_j=1^n A_ij  u_i(t,·)  u_j(t,·))]
        (x_k)  ∫_  [(∇^+_h T)(x - · )](x_k)  ϕ(x)  x 
       = ∑_k=0^M-1[∇^+_h (D_i  u_i(t,·) +  ∑_j=1^n A_ij  u_i(t,·)  u_j(t,·))] (x_k)  ∫_T(x-x_k) - T(x-x_k+1)/h ϕ(x)  x 
       = - ∑_k=0^M-1[∇^+_h (D_i  u_i(t,·) +  ∑_j=1^n A_ij  u_i(t,·)  u_j(t,·))] (x_k)  ∫_  T(x-x_k)   [∇^+_hϕ](x)  x.

By the discrete product rule (<ref>), this can be
estimated as

    ∫_0^T | _t ∫_ũ_̃ĩ(t,x) ϕ(x)  x |   t 
       ≤∫_0^T  ∑_k=0^M-1| [∇^+_h (D_i  u_i(t,·) +  ∑_j=1^n A_ij u_i(t,·)  u_j(t,·))] (x_k) |  ∫_  T(x-x_k)   |[∇^+_hϕ](x)|  x   t 
       ≤∇^+_hϕ_∞ h ∫_0^T  ∑_k=0^M-1| [∇^+_h (D_i  u_i(t,·) +  ∑_j=1^n A_ij  u_i(t,·)  u_j(t,·)) ] (x_k) |   t 
       ≤∇^+_hϕ_∞ ∫_0^T   h ∑_k=0^M-1∑_j=1^n A_ij  |u_i(t,x_k+h)| 
        | [∇^+_h u_j(t,·)](x_k) |  t 
          + ∇^+_hϕ_∞ ∫_0^T   h ∑_k=0^M-1∑_j=1^n A_ij  |u_j(t,x_k)| 
        | [∇^+_h u_i(t,·)](x_k) |  t 
          + ∇^+_hϕ_∞ ∫_0^T   h ∑_k=0^M-1 D_i  
        | [∇^+_h u_i(t,·)](x_k) |  t 
       ≤∇^+_hϕ_∞∑_j=1^n A_ij( ∫_0^T   h ∑_k=0^M-1  |u_i(t,x_k+h)|^2  t )^1/2  ( ∫_0^T   h ∑_k=0^M-1| [∇^+_h u_j(t,·)](x_k) |^2  t )^1/2
          + ∇^+_hϕ_∞∑_j=1^n A_ij( ∫_0^T   h ∑_k=0^M-1  |u_j(t,x_k)|^2  t )^1/2  ( ∫_0^T   h ∑_k=0^M-1| [∇^+_h u_i(t,·)](x_k) |^2  t )^1/2
          + ∇^+_hϕ_∞ D_i ( ∫_0^T   h  ( ∑_k=0^M-1 1 )  t )^1/2  ( ∫_0^T   h ∑_k=0^M-1| [∇^+_h u_i(t,·)](x_k) |^2  t )^1/2
       ≤ C_T   [n max_i,j A_ij + max_i D_i] ϕ_W^1,∞() [   T^1/2 + max_j ( ∫_0^T   h ∑_k=0^M-1  |u_j(t,x_k)|^2  t )^1/2 ] 
       ≤ C_T,

where we used estimates (<ref>) and (<ref>).



 §.§ Compactness


In order to stress the dependence w.r.t. the spatial discretisation, we denote by
[ũ_i]^M the interpolant associated to the discrete system on
Ω_M.

The classical Aubin-Lions lemma shows with the estimates
(<ref>), (<ref>) and (<ref>) that there
exists a subsequence such that

    [ũ_i]^M →_M →∞ u_i   strongly in L^4-([0,T] ×)

for some u_i ∈ L^4([0,T]×) ∩ L^2([0,T], H^1()).



 §.§ Passing to the limit


We now show that the limit is a solution in the weak formulation
stated in <ref>.

We first find that the interpolation [ũ_i]^M of the discrete
solution on Ω_M satisfies for all test functions φ :=
φ(t,x) ∈ C_c([0,+∞)×) and i=1,…,N that

    -  ∫_[ũ_i]^M(0,x)  φ(0,x)  x -
        ∫_0^∞∫_[ũ_i]^M(t,x) _t φ(t,x)  x  t 
        = D_i ∫_0^∞∫_ [ũ_̃ĩ]^M(t,x) Δ_h φ(t,x)  x  t + ∑_j=1^n A_ij∫_0^∞∫_[u_iu_j]^M(t,x)
         Δ_h φ(t,x)  x  t,

where [u_iu_j]^M is the interpolant (in the sense of
(<ref>)) of u_iu_j from the values on Ω_M.

Indeed, differentiating the interpolant [ũ]^M in time, we
find from the ODE system (<ref>) that

    _t [ũ_i]^M(t,x)
      = ∑_k=0^M-1[Δ_h( D_i u_i(t,·) + ∑_j=1^n A_ij u_i(t,·)  u_j(t,·))](x_k)  T(x-x_k).

Multiplying with the compact test function φ and integrating shows that

    -  ∫_[ũ_i]^M(0,x)  φ(0,x) dx-∫_0^∞∫_[ũ_i]^M_t φ  x  t 
       = ∫_0^∞∫_∑_k=0^M-1[Δ_h(D_i u_i(t,·) + ∑_j=1^n A_ij u_i(t,·)  u_j(t,·))](x_k)T(x-x_k)φ(t,x)  x  t
       =∫_0^∞∫_∑_k=0^M-1[D_i u_i(t,·) + ∑_j=1^n A_ij u_i(t,·) u_j(t,·)](x_k)Δ_h[T(x- ·)](x_k) φ(t,x)  x  t
       =∫_0^∞∑_k=0^M-1[D_i u_i(t,·) + ∑_j=1^n A_ij u_i(t,·) u_j(t,·)]
        (x_k)∫_T(x-x_k)(Δ_hφ)(x)  x  t
       = D_i ∫_0^∞∫_ [ũ_̃ĩ]^M(t,x)Δ_h φ(t,x)  x  t + ∑_j=1^n A_ij∫_0^T∫_[u_iu_j]^M(t,x)Δ_h φ(t,x)  x  t,

so that (<ref>) holds.


Next, we use the following result explaining how the linear interpolation behaves on products:


  Under the assumptions of <ref>, the following estimate holds:
  
    ∫_0^T ∫_|[u_iu_j]^M(t,x)-[ũ_i]^M(t,x)[ũ_j]^M(t,x)|  x  t ≤ C_T   h.



  For x ∈ [x_k,x_k+1) the representation formula (<ref>)
  shows for i=1,…,n that
  
    [ũ]^M_i(t,x)    = α_k(x) u_i(t,x_k) + (1-α_k(x))u_i(t,x_k+1),
    
        [u_iu_j]^M(t,x)    = α_k(x) [(u_iu_j)(t,x_k)] + (1-α_k(x))[(u_iu_j)(t,x_k+1)],

  where we recall that α_k(x)=(x_k+1-x)/h. Then, for
  x ∈ [x_k,x_k+1):
  
    |[u_iu_j]^M(t,x) -
          [ũ_i]^M(t,x)[ũ_j]^M(t,x)| 
       =|α_k(x) [(u_iu_j)(t,x_k)] + (1-α_k(x))[(u_iu_j)(t,x_k+1)] 
          -( [α_k(x)u_i(t,x_k) +  (1-α_k(x))u_i(t,x_k+1) ] [α_k(x) u_j(t,x_k) + (1-α_k(x))u_j(t,x_k+1)])|
       ≤α_k(x)(1-α_k(x))(u_i(t,x_k) + u_i(t,x_k+1) ) |u_j(t,x_k+1) - u_j(t,x_k)|.

  Consequently, we get that
  
    ∫_0^T ∫_|[u_iu_j]^M-[ũ_i]^M[ũ_j]^M |  x  t 
       ≤∫_0^T (∑_k=0^M-1∫_x_k^x_k+1α_k(x)(1-α_k(x))  x (|u_i(t,x_k)| + |u_i(t,x_k+1)|) |u_j(t,x_k+1) - u_j(t,x_k)| )  t 
       =1/6∫_0^T∑_k=0^M-1 h  (|u_i(t,x_k)| + |u_i(t,x_k+1)|)  |u_j(t,x_k+1) - u_j(t,x_k)|  t 
       ≤1/3∫_0^T (∑_k=0^M-1 h |u_i(t,x_k)|^2)^1/2(∑_k=0^M-1 h |u_j(t,x_k+1)-u_j(t,x_k)|^2)^1/2  t
       =1/3  h (∫_0^T ∑_k=0^M-1 h |u_i(t,x_k)|^2  t)^1/2(∫_0^T∑_k=0^M-1 h |(∇^+_h u_j(t,·))(x_k)|^2  t )^1/2.

  This conclude the proof with the estimates (<ref>) and
  (<ref>).


Thanks to the lemma above, the weak formulation (<ref>) of the
discretised system implies that for all
φ :=φ(t,x) ∈ C_c([0,T)×) and i∈ 1,…,n,

    - ∫_ [ũ_i]^M(0,·) φ(0,·)  x
      - ∫_0^∞∫_([ũ_i]^M _t φ + [ D_i  [ũ_i]^M + ∑_j=1^n A_ij [ũ_i]^M [ũ_j]^M ] Δ_h φ)  x  t = 𝒪(h).


We end up the proof of <ref> by observing that the limit
satisfies

    - ∫_ u_i(0,·) φ(0,·)  x
      - ∫_0^∞∫_(u_i_t φ + [D_i  u_i + ∑_j=1^n A_iju_iu_j] Δφ)  x  t = 0.

Indeed,

    | - 	∫_ u_i(0,·) φ(0,·) - ∫_0^∞∫_(u_i_t φ + [D_i   u_i + ∑_j=1^n A_iju_iu_j] Δφ)  x  t|
       ≤∫_|u_i(0,·)-[ũ_i]^M(0,
        ·)|  |φ(0, ·)| dx
        + ∫_0^∞∫_|u_i-[ũ_i]^M||_tφ|  x  t 
    
         +∑_j=1^n A_ij∫_0^∞∫_|[ũ_i]^M [ũ_j]^M-u_iu_j|   |Δ_h φ|   x  t
       +D_i ∫_0^∞∫_|[ũ_i]^M-u_i||Δ_hφ|   x  t
       +∫_0^∞∫_| D_i u_i + ∑_j=1^n A_ij u_iu_j||Δ_hφ - Δφ|  x  t +  𝒪(h).

The first integral tends to 0 because u_i(0,·) is continuous
on . The second and fourth integrals converge to 0 because
[ũ_̃ĩ]^M→ u_i strongly in L^1,

 the third integral converges
to 0 because [ũ_̃ĩ]^M→ u_i strongly in L^2, and the last
integral converges to 0 since φ is smooth.


This concludes the proof of <ref>.

elsarticle-num


