Convergence of a Fully Discrete and Energy-Dissipating Finite-Volume Scheme for Aggregation-Diffusion Equations

We study an implicit finite-volume scheme for non-linear, non-local aggregation-diffusion equations which exhibit a gradient-flow structure, recently introduced by Bailo, Carrillo, and Hu (2020). Crucially, this scheme keeps the dissipation property of an associated fully discrete energy, and does so unconditionally with respect to the time step. Our main contribution in this work is to show the convergence of the method under suitable assumptions on the diffusion functions and potentials involved.


Introduction
This work is concerned with the analysis of a finite-volume scheme for the family of aggregationdiffusion equations on Ω = (−L, L), for some L > 0, and for t > 0, together with a given non-negative initial datum ρ(0, x) = ρ 0 (x). (1.1b)

INTRODUCTION
The equation is equipped with no-flux boundary conditions: at x = ±L. The unknown, ρ, typically models the density of particles of a given species. The term H is called the internal energy density, and models linear or non-linear diffusion effects within the species ρ. The potential V models the drift of the species; it is sometimes referred to as the external potential, since it provides an energetic landscape that drives the density towards favourable regions. Finally, the potential W models the interactions of particles with each other; we call it the internal potential or the interaction potential.
It is important to highlight the generality of Eq. (1.1), since it encompasses many popular equations. In the absence of the potential terms, it becomes the well-known filtration equation. Linear diffusion can be obtained by choosing Boltzmann's entropy as the internal energy, i.e. H(ρ) = ρ log(ρ) − ρ; moreover, the linear Fokker-Planck equation is recovered by letting V (x) = x 2 /2. The porous medium equation [40] can be found instead by considering a power law for the internal energy density: H(ρ) = ρ m /(m − 1), for some m > 1. Typically, the non-linearity accounts for a density-dependent diffusion, as observed, for example, in porous medium flow [5], in heat transfer [42], or in population dynamics [28]. The interplay between diffusion and potential effects has a myriad of applications in life sciences; see [12] for a recent survey.
The choice of the potential W models certain phenomena, such as attraction, repulsion, and combinations thereof. Pairwise interactions often depend solely on the distance between agents; as a result, many applications assume the interaction potential to be radially symmetric, i.e., W (x) = w(|x|), for some w : R → R. This radial potential is called attractive if w > 0 and repulsive if w < 0. Typical choices include power laws in the context of granular media [6,39], or combinations of power laws (W (x) = |x| a /a − |x| b /b), which feature short-range repulsion and long-range attraction, in swarming and population dynamics [31,27,30,7,17]; these effects may appear simple at first glance, but they result in complex and subtle dynamics [13,4,25,30]. Other choices include local sensitivity regions (modelled by indicator functions [33,37,38] or attractive-repulsive Morse potentials [24,31,26,18,16]), as well as growth and saturation terms [2,34,41,21].
The equation, in the general form (1.1), has been studied in [11,19,20,32] through its dissipative structure. To be precise, the energy functional is dissipated along solutions of Eq. (1.1); formally, we may compute Numerical schemes which preserve the energy-dissipation structure are desirable to simulate the behaviour of Eq. (1.1). Suitable methods should also ensure the conservation of mass and preserve the non-negativity of the solutions. Such schemes were proposed in [22], for non-linear diffusion equations with drift, and in [11], for equations including non-local interactions. In the latter, the authors propose a semi-discrete, upwind finite-volume scheme that preserves the 1 INTRODUCTION entropic structure of the equation at the semi-discrete level. Moreover, by construction, the method conserves the total mass and preserves the non-negativity of solutions also. A fully discrete scheme with the same properties was proposed in [3].
The generalisation of Eq. (1.1) to the case of two interacting species was studied in [15,17] by extending the aforementioned semi-discrete method. Although the system generally lacks an entropic structure (even at the continuous level), the discretisation resembles strongly that of [11]. A fully discrete, implicit finite-volume scheme for the same system was additionally developed by [15]. The authors prove the convergence of their methods (semi-discrete and fully discrete) to weak solutions of the equations through a technique of flow interchange: they construct an auxiliary functional whose dissipation contains important gradient information. The construction for a general non-linearity H had already been employed in [14] in a semi-discrete setting, and it was adapted to prove convergence of the numerical method of [15].
We present a sketch of the estimate, adapted to Eq. (1.1), at the continuous level. Given a convex internal energy density, H, we construct a function K via H (s) = sK (s), for s > 0. We readily find: having used integration by parts together with the no-flux boundary conditions, and where α ∈ (0, 1) arises from Young's inequality.
To argue the convergence of the method of [3], we will construct a family of appropriate interpolations of the discrete numerical approximations, which will be identified as elements of an L ∞ class. The adaptation of the flow interchange argument (1.4), accompanied by some additional a priori estimates, will show the strong convergence of the approximating sequence to a weak solution of Eq. (1.1). As usual, these arguments do not yield a specific rate of convergence, though thorough numerical experiments in [3] suggest first-order accuracy. Nevertheless, we establish a convergence result for their numerical method, filling this theoretical gap in the literature.
Future work might address the question of passing to the limit in the energy-dissipation inequality. While it is clear that E(ρ(t)) ≤ E(ρ(s)) holds whenever t ≥ s, the analysis of the limit of the non-linear Fisher information remains open. Furthermore, it would be interesting to ascertain if the long-time behaviour of Eq. (1.1) under suitable convex potentials is preserved by the fully-discrete scheme. This behaviour is well-understood at the continuous level [19,1], but, to the best of our knowledge, discrete results are only known for particular cases [9,8].
The rest of this work is organised as follows: in Section 2 we recapitulate the numerical method, introduce the notion of weak solution, and present the main result; Section 3 is dedicated to establishing the existence of solutions to the numerical method as well as proving certain a priori bounds; in Section 4 we employ the estimates to deduce strong compactness of the approximating sequence, and weak compactness of the discrete velocity term (which comprises the internal energy, the drift, and interaction terms); Section 5 is dedicated to proving the main result: the convergence of the approximating sequence to a weak solution of Equation (1.1); finally, in Section 6 we adapt our convergence result to the case of linear diffusion; the Boltzmann entropy is not covered by the main result because the upwind discretisation of the scheme is incompatible with the infinite speed of propagation associated with linear diffusion; nevertheless, the argument can be modified to prove the convergence of the scheme.

Numerical Scheme
In this section we introduce the fully-discrete, implicit method of [3] for the problem (1.1) in one dimension, i.e., The equation is posed on the domain Ω := (−L, L) ⊂ R, and a time interval (0, T ). The equation is supplemented with no-flux boundary conditions and a given non-negative initial datum, ρ(0, x) = ρ 0 (x) ∈ L ∞ (Ω). Throughout this work, we will denote the space-time cylinder by Q T = (0, T ) × Ω. First, we introduce the discretisation of the domain.
, is centred at x i ; the cell centre is located at x i = −L + ∆x(i − 1/2), as shown in Fig. 1. The time domain is divided into N + 1 equal intervals of length ∆t = T /(N + 1), for some integer N ∈ N. The n th interval is defined by I n := [t n , t n+1 ), with t n := n∆t.
These partitions give rise to a computational mesh, {Q n i } which divides the space-time cylinder, Q T , into a family of finite-volume cells, denoted by Q n i = I n × C i , for n ∈ N := {0, . . . , N } and i ∈ I := {1, . . . , 2M }. The coarseness of a given mesh is measured by the mesh size h = ∆x = c∆t, for a fixed c > 0.
Finally, we define the dual cells, . We now proceed to discretise the problem (1.1). First, we construct the discrete initial datum by ρ 0 := {ρ 0 i } i∈I , on the discretised spatial domain through the cell averages of the continuous datum: where − Ci f (s) ds = 1 |Ci|´C i f (s) ds denotes the average of f on the i th cell C i . To discretise Equation (1.1a), we integrate it over a test cell, Q n i , which yieldŝ is the flux. This identity can be approximated by Giving rise to the numerical scheme for i ∈ I, and n ∈ N . The discrete solution, ρ n i , approximates the continuous solution in the finite-volume sense, The numerical fluxes, F n i+1/2 , are given by the upwind discretisation where (s) + := max{s, 0} and (s) − := min{s, 0}, for any s ∈ R, denote respectively the positive and negative parts of s. Choosing F n+1 1/2 = F n+1 2M +1/2 = 0 incorporates the no-flux boundary conditions into the scheme.
In order to define the velocities, u n+1 i+1/2 , it is convenient to introduce the discrete entropy variables, and to set The contributions from the external and interaction potentials are discretised as At this stage, we highlight once more that the evolution of the density, ρ, is governed by an entropic part, which works to minimising an internal energy, as well as a potential part, that generates a drift. For convenience of the analysis, we will hereafter split these contributions at the discrete level: (2.2c) Remark 2.2 (Choice of ρ * * ). The density present in the interaction term, ρ * * , may be chosen as with i ∈ I and n ∈ N . As discussed in [3,Theorem 3.9], this choice depends on certain properties of the potential and serves to establish a discrete analogue of the energy dissipation (1.3). In fact, it can be shown that the discrete energy is unconditionally dissipated if (i) ρ * * = (ρ n + ρ n+1 )/2; (ii) ρ * * = ρ n and the potential W is negative definite; (iii) ρ * * = ρ n+1 and the potential W is positive definite.
The interaction potential, W , is called negative definite (resp. positive definite) if for any two vectors η and ζ. For the readers' convenience, we reiterate how the discrete energydissipation inequality is obtained in Section 3.1.
In the sequel, we will only consider the most general case: ρ * * = (ρ n + ρ n+1 )/2. Nonetheless, the analysis is the same for the other choices of ρ * * , since the quantity only plays a small role in the fixed-point argument of Theorem 3.1 and appears inside the (conserved) L 1 -norm elsewhere.
Before proving the analytical properties of scheme (2.1), we introduce the necessary notation for the subsequent sections.
for i ∈ I and n ∈ N . Moreover, for any function η : R → R, we use the notation η n i := η(ρ n i ), as well as Given a quantity defined at the cell interfaces, x i+1/2 , for i = 1, . . . , 2M − 1 (such as the velocity u n i+1/2 and the flux F n i+1/2 ), we define its associated piecewise constant interpolations as We now establish a bridge between the discrete approximations from scheme (2.1) and the L p -functions through the piecewise constant interpolations, using the standard norms: for t ∈ I n , as well as |η n i | p ∆x∆t 1/p ; these are analogously described for quantities defined on the dual mesh.
Last, but not least, we define the discrete gradients, which occur naturally at the cell interfaces in our analysis. Definition 2.4 (Discrete gradients). We define the discrete gradient, d x η h , for a function, η h , as for (t, x) ∈ Q n i+1/2 and 1 ≤ i ≤ 2M − 1. In the same vein, the discrete gradient of a quantity defined on the dual mesh is given by for (t, x) ∈ Q n i and i ∈ I. Note that the values of ζ n 1/2 and ζ n 2M +1/2 should either be clear from the boundary conditions or otherwise prescribed.
We conclude this section by introducing our notion of weak solution and stating the main convergence result.
Theorem 2.6 (Convergence of the scheme). Suppose that V, W ∈ C 2 ([−L, L]), and also that H ∈ C 2 ([0, ∞)), such that H (s) > 0 for all positive values of s, as well as . Moreover, suppose that there exists a function, K, bounded from below, such that sK (s) = H (s). Then, for any non-negative initial datum ρ 0 ∈ L ∞ (Ω), we find that: (i) scheme (2.1) admits a non-negative solution that preserves the initial mass regardless of the mesh size; (ii) under the condition C (iii) the limit is a weak solution to Equation (1.1) in the sense of Definition 2.5.
Remark 2.7 (Assumptions throughout this work). In the sequel, we assume that the assumptions of Theorem 2.6 are satisfied, unless otherwise stated. ). This property is employed in the proof of Proposition 3.5, in order to ensure that the function K is defined at zero. In terms of common choices of H, the assumption rules out linear diffusion, corresponding to H(ρ) = ρ log(ρ) − ρ, and porous medium diffusion, The result of Proposition 3.5 can nevertheless be obtained for the entire porous medium range, H(ρ) = ρ m /(m − 1) for m > 1, by relaxing the assumption to lim →0 H ( ) = 0. The proof, given in Lemma 3.6, handles the singularity by defining a regularised form of K .
The theorem also assumes H ∈ C 2 ([0, ∞)), because our proof of compactness (Lemmas 4.2 and 4.3 in particular) requires H to be locally Lipschitz. Independently of the previous discussion, this again excludes linear diffusion, as well as porous medium diffusion for m < 2, all of which fail at the origin. Consequently, with the additional work in Lemma 3.6, the result in Theorem 2.6 also applies to the power law case m = 2, relevant in mathematical biology applications [17,21].
The assumption on H can be relaxed further if, in turn, the initial datum is strictly positive. Given ρ 0 > 0, Proposition 3.8 guarantees the positivity of the entire solution. The proof is generalised to include linear diffusion and all the porous medium cases, bypassing the issues at the origin altogether. This extension is immediate for porous media, and the case of linear diffusion in detailed in Section 6.
The compactness proof in Section 4 relies on the local Lipschitz continuity of H as well as the bounds in Proposition 3.8. These bounds are general for the scheme with periodic boundary conditions, but require an additional assumption to handle the no-flux boundary conditions; namely, a sign on the normal spatial derivative of (V + W * ρ) at the boundary of Ω. A nonnegative (respectively non-positive) sign guarantees the upper (resp. lower) bound. Remark 2.9 (Periodic boundary conditions). The analysis and numerical solution of Equation (1.1) are often performed in the setting of no-flux boundary conditions, as is done in this work. Nevertheless, all of our results can be immediately generalised to the case of periodic boundary conditions. This is because our analysis only exploits the lack of boundary terms when performing integration/summation by parts.

A Priori Estimates
The purpose of this section is threefold. We begin by proving existence of non-negative solutions to scheme (2.1), showing that the numerical method preserves the mass, and stating the discrete energy-dissipation property. Next, by mimicking estimate (1.4) at the discrete level, we obtain discrete gradient information of H , which, in return, will yield strong compactness of the density, ρ h , itself. We conclude the section by establishing uniform L ∞ -bounds and a control from below of the solution.

Existence of Solutions and Energy Dissipation Property
We proceed to show existence through a fixed-point argument.
Theorem 3.1 (Existence of solutions). Let ρ 0 (x) ∈ L 1 (Ω) be a non-negative initial datum for Equation (1.1). Then, there exist a unique, non-negative solution {ρ n i }, (i, n) ∈ I ×N , to scheme (2.1). Furthermore, the mass of the solution is conserved, i.e., Proof. The existence proof is based on an application of Brouwer's fixed-point theorem. To this end, we consider the compact, convex set define a function G : X → R 2M , and prove that this function is, indeed, a fixed-point operator.
Suppose that the solution of scheme (2.1), {ρ n i } i , is known at time t n , for some n ∈ N . The function, G, is defined through the implicit relationθ = G(θ), wherẽ for i ∈ I. The fluxes (which would coincide with the numerical fluxes at the fixed point) read where i = 1, . . . , 2M − 1, andF 1/2 =F 2M +1/2 = 0, to account for the boundary conditions. Note that the velocity terms are evaluated at the argument of the operator, θ, rather thanθ: It is easy to see that, as a composition of continuous functions, the operator G itself is continuous (recall that we operate under the assumptions of Theorem 2.6). Thus, it remains to show that it maps into the set X , i.e., that the image is non-negative and satisfies the mass constraint. We proceed by showing thatθ, the image of the operator, is non-negative, arguing by contradiction.
Suppose there are certain values i ∈ I such thatθ i < 0. Without loss of generality we assume that the corresponding cells C i form a single contiguous cluster. Therefore, we may assume that the negative values lie on a cluster of the form j Summing Eq. (3.1) over the pathological range and evaluating the flux terms as defined in The left hand side is strictly negative, whereas all the right hand side terms are non-negative, yielding the desired contradiction. We stress that this argument applies also if several clusters of this type were to exist and if the cluster is degenerate, i.e., j = k.
Having shown the non-negativity of the image, we note using Eq. (3.1) and no-flux condition,F i+1/2 =F i−1/2 = 0. As a consequence, G maps the set X into itself. Invoking this property in conjunction with the continuity of G, we may apply Brouwer's fixed-point theorem to infer the existence of a vector ρ n+1 satisfying ρ n+1 = G(ρ n+1 ).
It is readily verified that the fixed point satisfies scheme (2.1), and that the conservation of mass follows from the same argument as in Eq. (3.3).
Before proceeding with the analysis, we recall the aforementioned energy-dissipation property. As discussed in [3], the variational structure of Eq. (1.1) persists unconditionally in the discretised problem, thanks to the scheme, for a discrete analogue of the energy (1.2): The variation of this free energy over a single time step can be computed directly: having used the definition of ρ * * = (ρ n + ρ n+1 )/2 and the convexity of H. Substituting the scheme, we obtain having used the definition of the discrete entropy variables, ξ n+1 i , the velocities u n+1 i+1/2 , and a discrete integration by parts (see Lemma A.2 in the Appendix). A straightforward simplification then yields which is a discrete version of Eq. (1.3).

Discrete Gradient Estimate
We now turn to the obtention of a discrete gradient estimate in the spirit of Eq. (1.4). A crucial step in this endeavour is the construction of an auxiliary functional whose dissipation along solutions of Eq. (1.1) can be related to an L 2 (Q T )-bound on ∂ x H (ρ).
We remark that the convexity of K follows from that of H. Furthermore, a choice which is bounded below can always be found by adding a linear function to a second primitive of ρ −1 H , since convex functions lie above their tangents.
This associated functional can be used to define an average which generalises that of [15]. It has also been used, in a more specific form, in the case of the Boltzmann entropy [36,23]. x ≤ µ H (x, y) ≤ y.
Proof. Without loss of generality, let 0 < x < y. Using the definition of the auxiliary functional we find Upon rearranging, we deduce x < µ H (x, y). Similarly, which yields y > µ H (x, y). The case x = y is trivially true, and the case of x = 0 nevertheless holds because of the assumption that s −1 H (s) ∈ L 1 loc ([0, ∞)).
Before the gradient estimate, we show that the potential terms in the velocity can be controlled.
where C Proof. The confinement potential term is readily controlled by the derivative of V , since The interaction potential may be bounded in a similar fashion: which concludes the proof.
We are now ready to reproduce the estimate (1.4) in order to obtain an L 2 -bound on the discrete gradient of H .
for any α ∈ (0, 1), where C V is the constant from Lemma 3.4.
Proof. The convexity of K, together with scheme (2.1), yields The last equality follows from a discrete integration by parts (Lemma A.2) and the no-flux boundary conditions.
Just as in the continuous estimate, cf. Eq. (1.4), we expand the flux term: using the convexity of the auxiliary functional, K. In particular, we may choose the generalised entropic average,ρ i+1/2 = µ H (ρ n+1 i , ρ n+1 i+1 ), from Lemma 3.3, and conclude using the definition of µ H . Here, the remainder is given by The remainder term can be estimated using the weighted Young's inequality: using the bounds from Lemma 3.4. Applying this to Eq. (3.6), we finally obtain which proves the statement.
Lemma 3.6 (Relaxed assumption on H). The bound of Proposition 3.5, Proof. In order to avoid the possible singularity at zero, we define the regularised auxiliary quantity, K , by This allows us to define a regularised entropic mean, which is shown to satisfy in the same vein as the proof of Lemma 3.3. To reproduce the estimate of Proposition 3.5, we note that K remains bounded. Studying the dissipation of K along the discrete solution, we find cf. Eq. (3.5). The first term of the parenthesis is estimated as in Proposition 3.5. The second, which arises from the regularisation, vanishes in the limit, as we proceed to show. To this end, let us note that where While this constant depends rather poorly on ∆x, we remark that, for our purpose, it is sufficient to bound the quantity in Eq. (3.7) for fixed ∆x, ∆t, and pass to the limit → 0. The first term on the right-hand side of the inequality is bounded, and thus vanishes in the limit. To control the second term, we recall that where the first term on the right-hand side vanishes in the limit. Choosing δ = s = /2, we find Under the original assumption, s −1 H (s) ∈ L 1 loc ([0, ∞)), the limit is zero. In the relaxed setting, where the integrand might not integrable, we nevertheless write using L'Hôpital's rule and the relaxed assumption, which completes the proof.
. Let {ρ n i } be a solution to scheme (2.1). Then there exists a constant C ∂xH > 0, independent of the mesh size, such that Proof. Performing discrete time integration of the result of Proposition 3.5 yields Using C K , the lower bound on K (see Definition 3.2), we obtain Choosing α ∈ (0, 1) concludes the proof.

Upper and Lower Bounds
We now turn our attention to deriving bounds to control the approximate solution ρ h above and below in terms of the initial datum.
Proof. We will consider the solution to the scheme at two subsequent time instances, t n and t n+1 , and compare them. Let i 0 ∈ argmax i∈I {ρ n+1 i }, and note that ρ n+1 i ≤ ρ n+1 i0 for all i ∈ I.

The scheme yields
Rewriting the numerical flux at the right cell interface, we observe We recall the definition of scheme (2.1) and split the velocity term u n+1 i+1/2 into the entropic part and the drift part as defined in (2.2). This leads to having, once again, used ρ n+1 i0+1 ≤ ρ n+1 i0 , in conjunction with the monotonicity of H . An analogous computation readily shows that through the same technique employed in the proof of Lemma 3.4; we obtain assuming ∆tC which concludes the proof of the first bound. The second bound is proven in a parallel way. Let i 0 ∈ argmin i∈I {ρ n+1 i }. We observe and, through the splitting of the velocity, find having used the convexity of H. Estimating F n+1 i0−1/2 in a similar fashion, we obtain with C (2) V > 0 as above. Rearranging and iterating the inequality yields the statement and concludes the proof.

Compactness
Equipped with the a priori estimates obtained in the preceding section, we are now ready to establish the weak compactness of the velocity terms as well as the strong compactness of the density, ρ h . To this end, we first prove the strong compactness of the entropic term, H (ρ h ).
The splitting induces the following space discretisation of the time translates: where as well as see Fig. 2 for detail. The first double sum is either 0 (when τ < ∆t, which gives K = 0) or, using the Lipschitz continuity of H with constant C H , having used scheme (2.1) in the last line. Note that C H is not the same as C ∂xH , the L 2 bound on the discrete gradient of H (ρ h ) from Corollary 3.7. Upon summation by parts and rearranging the sums we obtain which can be controlled by substituting the flux and using the L ∞ -bound from Proposition 3.8: Thus we may write where we introduced the notation We begin by bounding the term, I ∆h (k), by splitting the velocity as in (2.2): using Young's inequality on each term. Extending the sum over n to the set N permits the combination of the terms involving H , where the last line follows using the bounds of Lemma 3.4.
Treating the second term of Eq. (4.2), I ∆h (k), in an identical fashion, we may use its bound in conjunction with that of Eq. (4.3) in Eq. (4.2) to obtain a bound on the total contribution: ∆h (k) ∆h ≤ CK∆h, Using a similar argument, we find Substituting Eq. (4.4) and Eq. (4.5) into Eq. (4.1), we obtain where the last equality holds because ∆t = ∆l + ∆h and K∆t + ∆l = τ .
Proof. Fix z > 0, and suppose L − z ∈ C 2M −K . Repeating the argument of the proof of Lemma 4.2, we see K = z ∆x and K∆x ≤ z ≤ (K + 1)∆x. Once again we recover a partition, this time of the spatial interval, ∆x = ∆l +∆h; we have ∆l = z −K∆x and ∆h = (K +1)∆x−z. This splitting yields the discretisation where and To bound the term I ∆l in Eq. (4.6), we use the Lipschitz continuity of H : This expression can be controlled in terms of the estimates from Corollary 3.7 and Proposition 3.8: (4.7) The second term in Eq. (4.6), I ∆h , is controlled similarly. In fact, It is easy to see that the compactness of ρ h follows directly from the strong compactness of (H (ρ h )) h . Theorem 4.4 (L 2 (Q T )-compactness of ρ h ). Let (ρ h ) h be the piecewise constant interpolations associated to a sequence of solutions to scheme (2.1). Then, the family converges strongly (up to a subsequence) in the L 2 (Q T )-sense to a function ρ ∈ L 2 (Q T ).

Convergence of the Scheme: Proof of Theorem 2.6
This section is devoted to proving the main result, Theorem 2.6. Having established all necessary estimates, we are prepared to prove the convergence of the piecewise constant interpolations, (ρ h ), associated to scheme (2.1) to weak solutions of Equation (1.1) in the sense of Definition 2.5.
To begin, let ϕ be a smooth test function such that ϕ(T ) = 0. We define the error term In the spirit of [15], it remains to show that (h) → 0, by comparing Eq. (5.1) with the scheme, thereby proving that ρ is, indeed, a weak solution.
In order to prove this claim, we define the following cell averages of the test function, ϕ, Multiplying scheme (2.1) by ϕ n+1 i and integrating, we obtain where we have used summation by parts, ϕ N +1 i = 0, and the no-flux boundary conditions, i.e., F n+1 1/2 = F n+1 2M +1/2 = 0. As before, we manipulate the flux term and observe that Using this, we note that Eq. (5.2) can be written as In a similar fashion, we may rewrite Eq. (5.1) as with terms given by Already noticing the resemblance, we will show that as h → 0.
The Time Term. The first term is exactly zero since The Entropic Term. We observe with remainder given by This can be readily controlled as by using the Cauchy-Schwarz inequality on both terms, as well as the fact that Having dealt with the remainder, we find The first term in the sum can be shown to vanish in the limit by noting that for some constants θ i+1/2 ∈ C i+1/2 ,θ i ∈ C i ,θ i+1/2 ∈ C i ∩ C i+1 . Thus we control the term by which vanishes as the mesh size goes to zero.
Handling the second sum in a similar fashion, we see , through the use of Cauchy-Schwarz. As discussed previously, this integral term tends to zero with the mesh size, due to the convergence of ρ. Hence, we conclude |H(h) −Ĥ(h)| → 0, as claimed.
The Potential Terms. The prior computation applied to the external potential term yields with remainder term which is controlled by the analogue of the bound for the previous term. This is shown by employing the L ∞ -bound on V in place of the L 2 bound on the entropic term. Following the same strategy, we find Each of the sums is shown to tend to zero as the mesh size goes to zero, just as previously. Thus we see |V(h) −V(h)| → 0.
The discussion of the interaction term is identical: Each of the sums (as well as the remainder R W (h)) are shown to vanish in the limit, proving The Error Term. Finally, we show a control onÊ(h) using some of the previous estimates: splitting once again the velocity term as was done in the proof of Proposition 3.8, , by the repeated application of the Cauchy-Schwarz inequality, where a constant independent of the mesh (see Lemma 3.4 and Corollary 3.7). We thus find , a term which vanishes with the mesh size h due to the convergence of ρ h .
Convergence. Combining the results of the preceding steps proves that (5.3) holds indeed. Therefore, we conclude that (h) → 0, which proves the convergence of the scheme with limit a weak solution to Eq. (1.1) in the sense of Definition 2.5, as claimed in Theorem 2.6.

The Case of Linear Diffusion
The discussion of the convergence of scheme (2.1) cannot be complete without addressing the case of linear diffusion. The choice of Boltzmann's entropy, H(ρ) = ρ log(ρ)−ρ, renders Eq. (1.1) into the linear drift-diffusion equation: Equation (6.1) is within the purview of [32,35,19,20], and the variational structure (1.3) persists. However, the choice of entropy H(ρ) = ρ log(ρ) − ρ cannot be addressed by the main result, Theorem 2.6. The derivative H (ρ) = log(ρ) lacks Lipschitz continuity at zero, which causes the estimates of Lemmas 4.2 and 4.3 to break down. In fact, the problem is even more blatant: no matter the approach, one must control the quantity ρ∂ x log(ρ) in order to reproduce the estimates of Section 5. At the continuous level, this is equivalent to controlling ∂ x ρ; however, in the discrete case, the quantity appears as which may be unbounded. The only viable approach is to apply the scheme to positive initial data, ρ 0 (x) ≥ ε for all x ∈ Ω, for some ε > 0. Using the propagation of lower bounds from Proposition 3.8, we will infer the boundedness of the solution away from zero for all times, ρ n i ≥ κ −1 > 0, which permits the control of terms like (6.2) through Rather than obtaining the compactness of the family H(ρ h ) through the auxiliary functional K, we will obtain the compactness of (ρ h ) directly through the classical functional ρ 2 L 2 . To be precise, we find where the last line is obtained from the bounds on the potentials V, W , and the application of Young's inequality with parameter α ∈ (0, 1). This readily yields a bound d x ρ h 2 L 2 (Q T ) ≤ C ∂xρ , (6.4) in the style of Corollary 3.7. Armed with gradient information, we find bounds for the shifts in time and space of ρ h which are equivalent to those of Lemmas 4.2 and 4.3. The proof of the space shift is identical, and the time shift argument is parallel until the term (6.2) appears, which is readily controlled as in (6.3). For completeness, we sketch the idea while neglecting the potential terms, which can be incorporated using Young's inequality. After introducing the discretisation from Eq. (4.1), we find Using the control from (6.3), we observe ρ n+k i (h n+k i+1/2 ) + + ρ n+k i+1 (h n+k i+1/2 ) − ≤ 2C ∞ κ|d x ρ n+k i+1/2 |, which simply yields (|d x ρ n+K i+1/2 | + |d x ρ n i+1/2 |)|d x ρ n+k i+1/2 |∆t∆x∆h.
Expanding the product term by term, using Young's inequality, and extending the sum over n, we finally obtain which is analogous to the bound found in Lemma 4.2, and permits the completion of the proof in the same fashion. The compactness of (ρ h ) on Q T now follows from the Frechet-Kolmogorov-Riesz theorem, and the weak compactness of the derivatives, (d x ρ h ), once again follows from Proposition 1 in [15]. This is enough to pass to the limit by reproducing the estimates of Section 5, using the newly found estimates (6.3) and (6.4), proving that (5.3) holds indeed. Therefore, we conclude the convergence of the scheme: Theorem 2.6 revisited (Convergence of the scheme). Suppose that H(ρ) = ρ log(ρ) − ρ and that V, W ∈ C 2 ([−L, L]). Then, for any strictly positive initial datum ρ 0 ∈ L ∞ (Ω), we find that: (i) scheme (2.1) admits a strictly positive solution that preserves the initial mass regardless of the mesh size; (ii) under the condition C V ∆t < 1 (viz. Proposition 3.8), the associated piecewise constant interpolation, ρ h , converges strongly in any L p (Q T ), for 1 ≤ p < ∞, up to a subsequence; (iii) the limit is a weak solution to Equation (1.1) in the sense of Definition 2.5. Remark 6.1 (Non-strictly positive data). If the initial datum ρ 0 is positive but not strictly so, we may approximate it by ρ ε 0 (x) := ρ 0 (x) + ε, for some ε > 0. Then, choosing ε and the mesh size h sufficiently small, and using the continuity of the aggregation-diffusion equation (6.1) with respect to the datum, we nevertheless obtain the convergence of the numerical scheme.