{'fraction_non_alphanumeric': 0.05460237338316305, 'fraction_numerical': 0.019147232266362512, 'mean_word_length': 3.9762106025842288, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An entropy structure preserving space-time formulation for cross-diffusion systems: Analysis and Galerkin discretization

Cross-diffusion systems are systems of nonlinear parabolic partial differential equations that are used to describe dynamical processes in several application, including chemical concentrations and cell biology. We present a space-time approach to the proof of existence of bounded weak solutions of cross-diffusion systems, making use of the system entropy to examine long-term behavior and to show that the solution is nonnegative, even when a maximum principle is not available. This approach naturally gives rise to a novel space-time Galerkin method for the numerical approximation of cross-diffusion systems that conserves their entropy structure. We prove existence and convergence of the discrete solutions, and present numerical results for the porous medium, the Fisher-KPP, and the Maxwell-Stefan problem.

1. Introduction. In this paper we develop a new space-time approach to the celebrated boundedness by entropy method by Ansgar Jüngel [26]. For a textbook version see [27]; see also [11,31].
Cross-diffusion systems are systems of nonlinear parabolic partial differential equations that are commonly used to describe dynamical processes appearing in modeling, for example, population dynamics, ion transport through nanopores, tumor growth models, and multicomponent gas mixtures. The challenge in the analysis of these systems is that the diffusion matrix is not necessarily symmetric nor positive semi-definite, and thus no maximum principle is available. Following [26], the remedy is to make use of the entropy structure of the system. Introducing the entropy function, a transformation of the solution, allows us to examine long-term behavior and show that the solution is nonnegative and bounded. Here, we present a space-time approach to the proof of existence of bounded weak solutions of cross-diffusion systems. The main tool will be the method of compensated compactness, which is a special technique of applying the classical div-curl lemma [48]. The key difference to the existing literature is that we do not make use of time-stepping, but instead consider time and space altogether. This naturally leads to a novel space-time Galerkin method for the numerical approximation of cross-diffusion systems. The space-time approach entails test and trail spaces, as well as the mesh, where time is included as additional dimension. This provides an easy way to increase the approximation degree simultaneously in space and time, and makes space-time hp-refinement possible. In a schematic way, our overall approach consists of the following four steps: 1. space-time variational formulation, 2. transformation to entropy variables, 3. regularization with a space-time H 1 inner product, 4. Galerkin discretization.
Existing numerical schemes for cross-diffusion systems rely on time-stepping methods. An entropy/energy conserving time-stepping algorithm for thermomechanical problems was developed in [41] being of second order in time. In [32], assuming existence of sufficient regular strong solutions on some time interval [0, T ] of a scalar diffusion equation, Runge-Kutta methods were studied using maximal regularity. Although maximal regularity also applies to a certain type of cross-diffusion systems [42], Runge-Kutta methods were only applied to very restrictive classes; an example (semi-discrete Runge-Kutta scheme) can be found in [29]. In [25], an entropy diminishing/mass conserving fully discrete variational formulation for a cross-diffusion system was presented. An alternative discretization for cross-diffusion systems based on the change to entropy-variables has been proposed in [16], where a dissipation-preserving approximation by Galerkin methods in space and discontinuous Galerkin methods in time has beed studied.
Maxwell-Stefan systems, see [37,46], describe multicomponent diffusive fluxes in non-dilute solutions or gas mixtures, and are a prime example for the cross-diffusion systems considered here. The first result on global solutions for the Maxwell-Stefan equations close to the equilibrium is given in [23]. The global existence of solutions close to equilibrium and the large time convergence to this equilibrium can be found in [21,Chapter 9], [22,24], and [42,Chapter 12]. The proof of existence of local classical solutions to the Maxwell-Stefan equations can be found in [6]. For a textbook on this topic, see [42]. The fact that the Maxwell-Stefan equations satisfy the assumptions made in this paper, see (H1)-(H3) below, is due to [30], where the entropy structure of the Maxwell-Stefan system was used to prove the existence of globally bounded weak solutions. An entropy structure was also identified for a generalized Maxwell-Stefan system coupled to the Poisson equation in [28], where the existence of global weak solutions was proven as well. The unconditional convergence to the unique equilibrium for given mass was shown in [24,36] without reaction terms. Those results were extended to also include reaction terms using mass-action kinetics in [13], whenever a detailed balance equilibrium exists. The heat equation can be recovered from the Maxwell-Stefan equation as a relaxation limit [43].
As of yet, numerical schemes for the Maxwell-Stefan equations commonly employ timestepping. A finite differences approximation can be found in [34,35]. Fast solvers for explicit finite-difference schemes were studied in [20]. A posteriori estimates for finite elements in the stationary case are given in [10]. In [40], a mass conserving finite volume scheme was presented. Existence of solutions for a mixed finite element scheme under some restrictions on the coefficients was proven in [38]. The scheme of [14] was proven to also conserve the L ∞ bound by making use of a maximum principle. A scheme using finite elements in space and implicit Euler in time was used to approximate a Poisson-Maxwell-Stefan system in [28]. That scheme, which is based on a formulation in entropy variables, admits solutions that conserve the mass as well as the entropy structure. As a by product, the solution satisfy an L ∞ bound. Another scheme that is mass conserving and conserves the L ∞ bounds of the solutions was presented in [8].
On simultaneous space-time finite element approaches for parabolic problems, there is a rich literature on the linear case, focusing on the heat equation; see, e.g. the recent overview [33,Ch. 7]. We point out that due to the different orders of derivatives present, conforming discretizations are typically based on a Petrov-Galerkin approach; see e.g. [2,4,45,47,50]. For nonlinear parabolic equations, the existing literature on space-time methods is much sparser. The adaptive finite element scheme introduced in [17] for linear parabolic problems was extended in [18] to the scalar version of the nonlinear reaction-diffusion equation treated in this paper. A space-time discontinuous Galerkin method for scalar nonlinear convection and diffusion was introduced in [49]. A space-time method for nonlinear PDEs using adaptive wavelets was introduced in [1].
The structure of this paper is as follows. In section 2, we state the problem and make the necessary assumptions for the existence of an entropy function. In section 3, we present the space-time Galerkin method on a regularized formulation of the problem in the entropy variable unknown, and state our two main results in Proposition 3.3 and Proposition 3.4, namely existence and convergence of discrete solutions, respectively. Existence of discrete solutions is proven in subsection 3.1. The proof of convergence will be split into two parts, first showing convergence with respect to mesh size in subsection 3.2, then proving convergence as the regularization parameter goes to zero in subsection 3.3. In subsection 3.4, we are then able to prove existence of a weak solution of the continuous problem. Numerical tests for the porous medium, the Fisher-KPP, and the Maxwell-Stefan problem are presented in section 4. All numerical results 1 were obtained using the finite element software NGSolve, see [44]. Additionally, in section 5, we reformulate the Maxwell-Stefan system with implicitly given currents in terms of the concentrations, and test it numerically.
Here, A(ρ) ∈ R N ×N is the diffusion matrix, f (ρ) : R N → R N represents the reactions, and ν is the outward pointing unit normal vector at ∂Ω; moreover, for 1 ≤ i ≤ N , We make the following hypotheses, which are slightly stronger assumptions compared to those made by A. Jüngel in [26].
, such that the following two conditions are satisfied: (H2a) There exists a constant γ > 0 such that Additionally, we make the following assumption on ρ 0 : (H3) The initial datum ρ 0 : Ω → D is measurable. A discussion on when it is possible to find a convex function s such that (H2) is satisfied for cross-diffusion equations can be found in [12] (see [12,Lemma 22]).
Remark 2.1. In [26], a more degenerate version of (H2a) is permitted where, in a nutshell, the coercivity inequality is replaced by z · s (ρ)A(ρ)z ≥ γ i ρ 2(m−1) i z 2 i for some m. In that context, an L 2 (0, T ; H 1 (Ω)) estimate for ρ i might be out of reach. Instead, the system is rewritten there by using ρ m i as an L 2 (0, T ; H 1 (Ω)) function. Moreover, also entropy densities s, which are not bounded, such as s(u) = u − log u, are allowed in [26]. As a consequence, a different version of the hypothesis (H2b) is considered there. We believe that our ansatz can be extended to that scenario applying the ideas from [26]. This will mainly affect the entropy inequalities and the proof of Proposition 3.8. However, we chose to use our simplified assumptions that already cover a large class of parabolic systems. By this, we try to keep the idea of our proof as fundamental as possible.
3. Space-time Galerkin method. Let the time T ∈ (0, ∞) be fixed. We denote by Q T = (0, T ) × Ω the space-time cylinder for a domain Ω ⊂ R d , d ≥ 1. We derive our method in four steps.
Step 1 (space-time variational formulation). The first step is to perform integration by parts in the time variable in (2.2), and to use the embedding which can be proved exactly as in [19,Chapter 5.9,Theorem 3]. Then, we define the following variational formulation of (2.1).
The following lemma, which will be proven in section 3.4 below (see Remark 3.11), establishes that the variational formulation in Definition 3.1 is actually equivalent to the one in (2.2).
Step 2 (transformation to entropy variables). In the second step, we express ρ in formulation (3.2) as a function of the entropy variable w, namely ρ = u(w). The resulting spacetime variational formulation is the following: Find w ∈ H 1 (Q T ) N such that Here, we use the notation φ(t) := tr(φ)(t, ·), where tr denotes the trace operator tr : Step 3 (regularization). Then, we introduce the following regularized problem with regularization parameter ε: The regularization term is given by the scaled H 1 (Q T ) N inner product defined as for f, g ∈ H 1 (Q T ) N . Its associated norm will be denoted by · H 1 ε (Q T ) .
Step 4 (Galerkin discretization). Finally, we discretize equation (3.4). Let {V h } h>0 be a family of finite dimensional spaces, parametrized by h > 0, such that, for every h, V h ⊂ H 1 (Q T ) N . We make the following approximability assumption on the family of spaces Therefore, we consider the following space-time Galerkin scheme in the entropy variable unknown: Find w ε h ∈ V h such that The first term in (3.5) can be interpreted as a stabilization term for the Galerkin scheme, with parameter ε > 0. This is used to obtain a control of the entropy variable. Note that, due to the nonlinearity of u, we expect that The following two propositions constitute the main result of this paper. Proposition 3.3 establishes that the Galerkin problem (3.5) admits solutions w ε h and these solutions satisfy an entropy estimate. Then, this is exploited in Proposition 3.4 to obtain existence of weak solutions ρ to the continuous problem (2.1), together with related entropy estimates. Propositions 3.3 and 3.4 will be proven in section 3.1 and section 3.4, respectively. Here and in the following, |Ω| denotes the volume of Ω, and γ and C f are as in Assumption (H2). Proposition 3.3 (Existence of discrete solutions). There exists a solution w ε h ∈ V h of method (3.5). Moreover, every solution w ε h ∈ V h of (3.5), for ε, h > 0, satisfies the entropy estimate Proposition 3.4 (Convergence). Let w ε h ∈ V h be a solution of (3.5) for ε, h > 0. Then there exist a weak solution for all r ∈ [1, ∞). Moreover, ρ satisfies the entropy estimate Remark 3.5 (Non-closed systems). The general setting given by (2.1) describes closed systems, i.e., without influx or outflux at ∂Ω. These systems are of particular interest, as they obey the second law of thermodynamics proving the decay of the entropy. However, in some cases (e.g. the lung model presented in section 5.3), one is interested in a non-closed subsystem involving, for instance, inhomogeneous Dirichlet boundary conditions on some part Γ D ⊂ ∂Ω of the boundary.
If ρ = g is prescribed on Γ D , for a given g ∈ H 1 (Ω) N taking values in D, we assume that the approximation spaces V h are affine subspaces of s (g) . Unfortunately, we can neither guarantee the discrete entropy estimate (3.6) nor its continuous version (3.7). Instead, we have to work with the relative entropy which is still convex in the variable ρ. We conjecture that, under the assumption s (g) ∈ L ∞ (Ω), one can prove corresponding versions of Proposition 3.3 and Proposition 3.4 by employing an estimate of the relative entropy of the form for all τ ∈ (0, T ) and some C A,g,f only depending on γ, g H 1 (Ω) and s (g) L ∞ (Ω) , as well as on A L ∞ (D) and f L ∞ (D) . This estimate can be derived from (2.2), with test function φ = s (ρ) − s (g). However, the proof of convergence of the discrete scheme (3.5) remains an open problem and is under ongoing investigation.

Existence of a solution of the numerical scheme.
Proof of Proposition 3.3. The idea is to use the Leray-Schauder fixed-point theorem for the mapping Φ : Since A, f, u are continuous, so is Φ. Since V h has finite dimension, Φ is also compact. Then by the Leray-Schauder fixed-point theorem, we obtain that Φ admits a fixed-point if we can show that the set ). Thus, by the fundamental theorem of calculus, Note that, by definition, s (u(w ε h )) = w ε h . The convexity of s then implies that The next step is to use (H2a) in combination with w ε h = s (u(w ε h )), which yields that Thus, the Leray-Schauder theorem is applicable and yields that Φ has a fixed point, and therefore the scheme (3.5) admits a solution. Using these calculations for σ = 1, it follows that every solution has to satisfy the entropy inequality (3.6).

3.2.
Convergence of the numerical scheme as h → 0. We will show that, for a fixed ε > 0, the numerical scheme (3.5) converges as h → 0.
. Moreover, w ε solves (3.4) and satisfies the entropy estimate Proof. The first part of the assertion follows from the fact that w ε h is uniformly bounded in due to the Banach-Alaoglu theorem, and w ε h → w ε in L 2 (Q T ) N , due to Rellich's theorem. In particular, we can choose this subsequence in such a way that w ε h converges a.e. to w ε . As u is bounded (see Assumption (H2)), the dominated convergence theorem entails the strong convergence of Combining this with the entropy estimate (3.6), there exists another subsequence (which we do not relabel) Using φ h as a test function in (3.5), we obtain (3.4) in the limit h → 0, as each integral in (3.5) converges separately. The entropy inequality (3.8) is a consequence of Fatou's lemma and the weak lower semicontinuity of the norm.
The following corollary will be used in the analysis of the limit for ε → 0 (see proof of Proposition 3.8 below).
Corollary 3.7. Let τ, δ ≥ 0 be such that τ + δ ≤ T . It holds true that Similarly as in the proof of Proposition 3.3, we use ρ ε := u(w ε ) and and, since ψ(T ) = 0 and ψ(0) = 1, Thus, using the definition of ψ, and treating the last term of the previous equation as in the proof of Proposition 3.3, we get From (3.4) tested with φ = ψw ε and the previous inequality, we get which, due to the properties of ψ and the assumption (H2), entails Finally, we can estimate the first term as Using the Cauchy-Schwarz inequality and the definition of the H 1 ε norm yields and therefore Note that we cannot estimate the first term on the right-hand side by the first term on the left-hand side, because the domain of the norms are disjoint. Fortunately, we have the entropy estimate (3.8), which we add √ ε/δ times to this inequality to get which, since s(ρ ε (T )) ≥ 0, implies the assertion.
3.3. Limit of ε → 0. We consider the limiting problem As above, we use the notation φ(t) := tr(φ)(t, ·), where tr denotes the trace operator tr : Moreover, ρ satisfies the entropy inequality In the proof of Proposition 3.8, the key ingredient to prove strong convergence of (at least a subsequence of) ρ ε will be the idea of compensated compactness, which is a special technique applying the classical div-curl lemma; see, e.g. [48,Lemma 7.2].
. Proof of Proposition 3.8. Let w ε , ρ ε := u(w ε ) denote the solution of (3.4) satisfying the entropy inequality (3.8). For any fixed i, i = 1, . . . , N , we define the vector-valued functions with (1 + d) components Note that, by assumption, D is bounded and so is ρ ε = u(w ε ). Thus, thanks to the entropy estimate (3.8), α ε , β ε are bounded uniformly in L 2 (Q T ) 1+d w.r.t. ε ∈ (0, 1). By the Banach-Alaoglu theorem, there exist α, β ∈ L 2 (Q T ) 1+d and a subsequence ε → 0 such that Clearly, β has the form (ρ i , 0) for some ρ i ∈ L 2 (Q T ). Due to the entropy estimate (3.8), . In order to apply the div-curl lemma, it remains to prove that the space-time divergence of α ε is bounded. For this, we require the equation for ρ ε i in the interior of Q T , i.e., from equation (3.4), . We can rewrite this by using the weak space-time divergence of α ε as We observe that the right-hand side defines a bounded operator on L 2 (Q T ) due to the entropy estimate (3.8) and the fact that f i is uniformly bounded as a continuous function defined on a compact set (see (H2)). This yields that div (t,x) α ε is uniformly bounded in L 2 (Q T ). Therefore, we can apply the div-curl lemma and obtain that In particular, there exists a subsequence not being relabeled such that ρ ε i → ρ i a.e. in Q T . For almost every (t, x) ∈ Q T , we know that ρ ε (t, x) ∈ D and that D is bounded. Thus, we can apply the dominated convergence theorem, which yields that and that ρ(t, x) ∈ D for almost every (t, x) ∈ Q T . Moreover, the entropy inequality (3.8) also states that ∇ρ ε i is bounded in L 2 (Q T ) d independently of ε. Since |ρ ε | = |u(w ε )| = |(s ) −1 (w ε )| ≤ sup v∈D |v| 2 , according to (H2), then, using 2 The fact that the L 2 norm of ∇ρ ε i is uniformly bounded is ultimately a consequence of our hypothesis (H2a), i.e., the matrix s (ρ)A(ρ) being coercive. Using instead the original assumptions made by A. Jüngel in [26], one can only assure that ∇ρ m i is bounded in L 2 for some m. However, one can circumvent this issue by defining β ε := (ρ m i , 0) T in this case.

Existence of a weak solution.
In this section, we prove that problem (2.1) possesses a weak solution ρ in the sense of Definition 3.1. Moreover, we prove the equivalence stated in Lemma 3.2 between the weak formulation (3.2) in Definition 3.1 and the weak formulation (2.2).
Corollary 3.12. Let ρ be given by Proposition 3.8. Then ρ is a solution of (2.2). Proof. Thanks to Proposition 3.10, we know that ρ possesses enough regularity such that we can integrate in (3.10) w.r.t. t, which yields (2.2) for all φ ∈ H 1 (Q T ) N with φ(T ) = 0. Using a density argument yields the assertion.

The proof of Proposition 3.4 is now straightforward.
Proof of Proposition 3.4. We only have to collect the previous results to obtain the proposition using a diagonal sequence argument.  4). For the Maxwell-Stefan system, the discussion of the general setting and of an alternative space-time Galerkin method for the case of N > 2 is postponed to section 5. We remark that we apply this nonlinear setting to the linear heat equation for validation purposes and, in particular, in order to stress its unconditional stability on a simple test problem.
In all cases, we consider the entropy density s : D → [0, +∞) defined by Then s ∈ C 2 (D, [0, ∞)) ∩ C 0 (D) and is convex. Moreover, u : R N → D defined as a choice first used in [9] to investigate the case N = 2, is in C 1 (R N , D), and is the inverse of s . Thus, the preamble of assumption (H2) is satisfied.
In the numerical experiments below, we use continuous space-time finite element discretization spaces. On the space-time cylinder Q T = Ω × (0, T ), with Ω bounded interval (d = 1) or Lipschitz polytope (d > 1), we consider families of shape-regular simplicial or Cartesian meshes {T h } h>0 .
The parameter h denotes the mesh granularity, namely T h = {K i , i = 1, . . . , N h }, h K := diam(K), and h := max K∈T h h K .
As discretization spaces, we choose where P p (K) denotes the space of polynomial functions on K of degree at most p, if K is a simplex, or of degree at most p in each variable, if K is a cuboid. Therefore, the approximability assumption (H4) in the first part of section 3 is satisfied. Defining the space-time Galerkin method (3.5) can be rewritten more explicitly in terms of the entropy variable unknown as follows: Throughout this section, we measure the absolute numerical error defined by ρ − u(w ε h ) L 2 (Q T ) . Remark 4.1 (On the choice of ε). For the analysis the regularization term with the factor ε > 0 is crucial, as it delivers essential bounds on the entropy variable w. In the numerical examples we will specify the choice of ε for every example, and show on several occasions that it is possible to choose ε = 0. However, it is crucial to note that in these examples the entropy variable representing the exact solution has 'nice' bounds. Conversely, when the solution approaches the singularities of the entropy, it is required to choose ε > 0 large enough for the solver to converge. In general, it is best to choose ε as small as possible, as we observe in the second example in subsection 4.2.
For the numerical tests, we take Ω = (0, 1) 2 and ρ 0 (x) = 0.5 cos(πx 1 ) cos(πx 2 ) + 0.5, so that the problem has the analytical solution given by ρ(t, x) = 0.5 exp(−8π 2 t/τ ) cos(2πx 1 ) cos(2πx 2 ) + 0.5, where we use τ = 7 to rescale the time. The solution is shifted and scaled in order to avoid the singularities of s at 0 and 1. Without this rescaling, the system matrix is highly ill-conditioned,  which prohibits optimal convergence rates. We solve (4.3), setting ε = 0 and solving the nonlinearity by Newton's method. We use unstructured space-time simplicial meshes.The Newton method converges in 6 steps, for all considered values of h and p. We measure the L 2 error on the whole space-time domain. In Figure 1, the convergence rates of the h-and the p-version of the method are shown. We observe optimal rates, exponential in p and of order p + 1 in h. In the case of p = 4, we observe a preasymptotic region for very large mesh sizes; the exact rates are shown in Table 1.  In Figure 2 we highlight another feature of the space-time approach, namely the ability to use adapted mesh in space and time. As a comparison, we use a mesh made of time slabs, a mesh structure similar to what would result from classic time-stepping methods. The timeslab height is given by h t ≈ h 2 x , with h x being the mesh size of the spatial mesh. For the spacetime adapted mesh, we start with a unstructured simplicial mesh of size h = 0.2 and then apply adaptive refinement. For this simple example, we use a flux based error estimator and Dörfler marking with θ = 0.5. We observe the same rate of convergence on both meshes, however using the space-time adapted mesh allows us to obtain a given accuracy with fewer degrees of freedom.

The porous medium equation.
Let m > 1. The porous medium equation is given by in Ω.
In Figure 3, we show the convergence rates of the scheme. Regardless of the nonlinearity, we match the convergence rates of the heat equation, i.e. exponential in p and of order p + 1 in h. The convergence rates in terms of h are also considered for different values of ε. We observe that ε introduces a lower bound on the error. Therefore, choosing it as small as possible, such that the solver still converges, gives the best results. On the other hand, in the next example, we can see that for certain solutions, that produce a very ill-conditioned system, we must choose ε fairly large. In contrast to the heat equation, the power law in the porous medium equation introduces a finite propagation speed of the solution. This is best observed by the interesting behavior of certain initial conditions that induce a waiting time. That is, the solution keeps a fixed support until the waiting time is reached. On Ω = (0, π), the initial condition given by produces this behavior. It is shown in [39] that the corresponding solution has a waiting time of t * = m−1 2m(m+1) . As we choose m = 2, here t * = 0.083. We modify the initial condition to ρ 0 (x) = 10 −16 for x / ∈ [0, π] to avoid ill-conditioning. Furthermore, to ensure convergence of the Newton method used as a nonlinear solver, we had to choose ε = 10 −6 , making use of the regularization term. We solve on a Cartesian space-time mesh until final time T = 0.2, with spatial mesh size h s = 0.05, and temporal mesh size h t = h s /2, and fix p = 5. The results are shown in Figure 4. Looking at snapshots of the numerical solution we can observe that it keeps a compact support set. In Figure 4, on the right, we plot the value of the solution on the left interface against time, marking the expected waiting time t * with the vertical line.
Choosing an entropy density such that assumption (H2b) is satisfied with C f = 0 allows for the right-hand side of the entropy estimate (3.6) to be independent of time. Motivated by this, we now investigate the rescaled entropy density s : D → (0, +∞) given by with n to be chosen. Note that f (ρ) > 0 for ρ ∈ (0, 1), and n/ρ − 1 > 1 if and only if ρ < n/2. Thus, for all for ρ ∈ (0, 1) if and only if n ≥ 2. We choose n = 2 so that the hypothesis (H2b) is fulfilled with C f = 0.
We start again by investigate convergence towards a smooth solution. We choose Ω = (0, 1), and initial conditions and Neumann boundary conditions such that is the exact solution for A = 1. We set ε = 0 and solve on unstructured simplicial space-time meshes. The results are presented in Figure 5. We observe again optimal convergence rates in both p and h, namely exponential in p and of order p + 1 in h. Next, we aim to reproduce the experiments presented in [5], considering an initial condition with a jump, given by ρ 0 (x) = 1 if 0 < x < 1/2 and 0 elsewhere, with diffusion coefficient A = 10 −4 . We solve using p = 3 on a Cartesian mesh with h s = 0.025, h t = 0.4 up to T = 8. We choose ε = 10 −8 to avoid ill-conditioning in the solver. Snapshots of the numerical solution are taken every 1.3 seconds, the results are shown in Figure 6 on the left. In Figure 6 on the right, we consider different choices for the entropy up to T = 15. Note that at the point in time the solution has already converged to ρ ≡ 1. The choice for the entropy density in [5] was ρ log(ρ) − ρ + 1. This choice is not covered by our assumptions, however, using it produces the correct results, as conjectured in Remark 2.1. We compare it to the entropy in (4.4) for different values of n in Figure 6. For the choice of n = 2, we recover a similar behavior of the entropy, namely, a region with slow decay followed by an exponential decay. As the solution converges to 1 it can easily be seen that for n > 2 the entropy does not converge to zero exponentially, as exemplified by the choice of n = 2.1 in the figure, since it is not the correct relative entropy with respect to the equilibrium.

4.4.
The three-component Maxwell-Stefan system. The Maxwell-Stefan system for N = 2 can be written as The unknowns ρ 1 and ρ 2 represent the concentrations of the first two gases (ρ 3 = 1 − (ρ 1 + ρ 2 )); the parameters d 1 , d 2 , and d 3 are the related to the binary diffusion coefficients of the three gases. In section 5 below, we derive this form of the Maxwell-Stefan system, prove that it fits our framework, and discuss the case N > 2.   In [7, Sec. 2] numerical results were presented for the three component gas diffusion experiment originally performed by Duncan and Toor in [15]. The setting is the following. Consider two spherical bulbs of volume 77.99 cm 3 (radius 26.49 mm) and 78.63 cm 3 (radius 26.58 mm), respectively, which are connected by a capillary tube of length 85.9 mm and diameter 2.08 mm, with a valve in the middle. We consider the Maxwell-Stefan equations with N = 2, corresponding to the gas mixture composed of hydrogen (ρ 1 ), nitrogen (ρ 2 ), and carbon dioxide (ρ 3 ). We consider the following initial gas mixture in the left-and right-hand side of the device: Left: (ρ 0 ) 1 = 0.000, (ρ 0 ) 2 = 0.501, (ρ 0 ) 3 = 0.499, For these gases, the coefficients d 1 , d 2 , and d 3 are given in terms of the binary diffusion coefficients (see section 5.1 below) as follows: As in [7], we can reduce the domain to two dimensions, as the device and initial conditions are axially symmetric and the flux vector has no angular component. In Figure 7, the computational domain is shown. We choose the spatial mesh size h s = 2.08 mm, equal to the diameter of the tube. The size of the Cartesian product mesh in time is chosen as 20.8 s. We solve iteratively on these slabs, restarting the computations with the previous solution as initial condition. We fix p = 2 and ε = 10 −10 .
The results are shown in Figure 8. We recover the same behavior shown in [7]. Both hydrogen and carbon dioxide converge monotonically to the expected equilibrium. Nitrogen shows the peculiar behavior known from the experiment. In Figure 9, we show the relative entropy and its dissipation, i.e. the time derivative of the entropy, both converge exponentially. This can be expected, as this behavior has been proven for similar types of equations, see for example [3, eq. (1.16)] where it is formally stated that the entropy dissipation converges exponentially for the Fokker-Planck equation by the use of the Bakry-Emery method.
5. The Maxwell-Stefan system revisited. In this section, we derive the formulation of the Maxwell-Stefan system as that used in section 4.4, and show that it fits into the general framework of section 2 (section 5.1). For the case N > 2, in which an explicit representation of the currents may not be easily derived, we introduce and analyze an alternative space-time Galerkin method, which is based on a formulation that is implicit for the currents (section 5.2). Let The Maxwell-Stefan equations are given by the continuity equations for i = 1, . . . , N + 1, where the currents J i are implicitly given by for some D ij = D ji > 0.

Explicit formula for the currents.
In this section, we establish an explicit representation of the currents, which allows us to derive the formulation of the Maxwell-Stefan system in the concentration variable unknowns. We follow [6] (see also [30]).
Using ρ i ≥ 0 and D ij = D ji > 0, it is easy to see that M (ρ) is quasi-positive (M ij (ρ) ≥ 0 for i = j). Moreover, provided that ρ i > 0 for all , is symmetric, thus all the eigenvalues of M (ρ) are real. By the Perron-Frobenius theory for quasi-positive, irreducible matrices, one deduces that the eigenvalue zero has multiplicity one (we refer to [6] or [30] for details). We deduce (5.3) Ker M (ρ) = span{ρ} and Im M (ρ) = v : As M (ρ) is not invertible, we have to restrict ourselves to a subspace of all possible currents J in order to obtain an explicit formula for J. For this, we make the assumption that the total current vanishes. Then by summing in (5.1) over all i = 1, . . . , N + 1, we see that Using this, we can rewrite the implicit formulation of the currents as As before, we can define a matrix , i, j = 1, . . . , N.
From (5.3), the matrix M(ρ) has full rank, and hence it is invertible. We have Proof. Let A(ρ) = −M −1 (ρ). The fact that M is smooth directly implies that A is smooth. Similarly as in the proof of [30,Lemma 3.2], one can show that for some γ > 0 and all smooth w.
In order to prove (H2a), we have to show that Let ρ ∈ D, x 0 ∈ Ω, and z ∈ R N . We define the following vector-valued function of x: whereê 1 denotes the unit vector (1, 0, . . . , 0) ∈ R d . We have and, for u = (s ) −1 , This, together with (5.6), implies that which proves the assertion.
For N = 1, the matrix M(ρ) is actually a scalar, which is given by Hence, J 1 = D 12 ∇ρ 1 . Therefore, in this case the Maxwell-Stefan system reduces to the heat equation.
For three species/gases (N = 2), we have and recall that D 21 = D 12 . One can verify that Let A(ρ) denote the inverse of −M(ρ). We can rewrite the Maxwell-Stefan equations as the system in section 4.4.

Implicit formulation for the currents.
In subsection 5.1, we have seen that the Maxwell-Stefan system (5.1)-(5.2), can be written in the form (2.1), with f = 0 and A(ρ) being given by the inverse of −M(ρ) for Moreover, we have computed A(ρ) explicitly for N = 1 and N = 2. However, for large N , it is more complicated to find the explicit formulation for A(ρ). In any case we do not expect a simple formulation in these cases. Therefore, this section provides a space-time Galerkin scheme, which avoids the explicit computation of the inverse of M. Let q, p ∈ N. We consider the following problem: Proof of Proposition 5.3. The idea of the proof is to proceed similarly to the proof of Proposition 3.3. We define the mapping where w is (uniquely) defined via the equation and J µ denotes the unique solution (see below for a justification) of Note that the mapping Φ is well-defined, as (5.8) admits a unique solution for given v ∈ V p h according to the Lemma of Lax-Milgram: we see that for y := A(u(v)) −1 z = −M(u(v))z. Moreover, the mapping Φ is continuous since A and u are continuous. Then by the Leray-Schauder fixed-point theorem, we obtain that Φ admits a fixedpoint if we can show that the set Using that σ ∈ (0, 1] is bounded from above yields a uniform bound on w in V q h and on M(u(w))J µ in L 2 (Q T ) N . As V q h is finite dimensional, we directly obtain that w L ∞ (Q T ) N is uniformly bounded. Thus, is also uniformly bounded. As all norms are equivalent on V q h , this directly implies that J µ is uniformly bounded in V q h . Thus, the Leray-Schauder theorem is applicable and yields that Φ has a fixed-point, and therefore the scheme (5.7) admits a solution.
Proof. The proof is analogue to the proof of Proposition 3.4. We only need to replace Proposition 3.6 by Lemma 5.6 below.
Proof. The fact that w h is uniformly bounded in H 1 ε (Q T ) N yields that there exists w ∈ H 1 (Q T ) N and subsequence h → 0 such that w h w in H 1 ε (Q T ) N , due to the Banach-Alaoglu theorem, and w h → w in L 2 (Q T ) N due to Rellich's theorem. As u is bounded, the dominated convergence theorem entails the convergence for ρ h ≡ u(w h ) to ρ := u(w) along another subsequence (which we do not relabel).
For the second part, we note that, due to the Banach-Alaoglu theorem and the boundedness of M(u(w h ))J µ h in L 2 (Q T ) N , we know that there exist ξ µ ∈ L 2 (Q T ) N such that, for a subsequence (not being relabeled), In particular, Using φ µ h as a test function in (5.7), in the limit h → 0, we obtain as each integral in (5.7) converges separately. In particular, by the fundamental lemma of calculus of variations, we see that ∂ xµ w = s (u(w))M(u(w))J µ and equivalently which implies that ρ solves (3.4). Finally, the entropy inequality is a consequence of Fatou's lemma.

Numerical Tests.
We again turn to [7,Sec. 3] for numerical results we can compare our method to. This time, we consider a model for the lung. The computational domain resembles a branch of the tree structure found in the bottom of the lung. The domain, depicted in Figure 10, consists of the inflow, Γ 1 , on top, the outflow, Γ 2 , located on the bottom of the two branches, and the alveoli, Γ 3 , located in the middle of each of the branches. The remaining boundary Γ 4 is a wall where nothing goes in or out. Opposed to the domain presented in the reference, we consider the branches of the lung to be symmetrical and perpendicular to each other. The paper does not mention the angle between the branches used there. Also the size of the alveoli is left  unspecified in the paper. Here, we split the boundary of the branches into three equal parts, with the alveoli (Γ 3 ) in the middle. On Γ 1 , Γ 2 , Γ 3 we impose Dirichlet boundary conditions to model the gas exchange with the other parts of the lung. On the wall, Γ 4 , we take homogeneous Neumann boundary conditions.
We make use of the implicit formulation (5.7) to find the numerical solution. To incorporate the Dirichlet boundary condition, we use Nitsche's method and add to (5.7) the following terms: for a parameter η > 0, h s being the spatial mesh size, on the Dirichlet boundary Γ D .
In the examples below, we use η = 1. The first term comes from the integration by parts. The second and third terms are productive zeros that weakly enforce the Dirichlet boundary condition, and are chosen such that they agree with Nitsche's method for the heat equation in the degenerative case.
5.3.1. Diffusion of air. In the following example, compare [7,Sec. 3.4], we choose alveolar air as initial condition and as the Dirichlet data on the outflow and alveoli. On the inflow boundary we choose humidified air as Dirichlet data. See Table 2 for the gas components of the different types of air, and Table 3 for the diffusion coefficients.  Since there is no helium present we can reduce the number of species involved, setting N = 3. For the numerical calculations we choose spatial mesh size h s = 0.3 and measure the value of the gas every 0.001 seconds. The discrete system is not ill-conditioned and we are able to choose ε = 0. In Figure 11 we show the numerical results for Oxygen and Carbon dioxide as the other gases stay (almost) constant. Both converge to their equilibrium value. Comparing the results to [7], we can see that the equilibrium value slightly differs, which is likely due to the symmetry of the domain and size of the alveoli.

5.3.2.
Diffusion of air/heliox. Next, we try to reproduce the results form [7,Sec. 3.5]. We consider alveolar heliox as initial condition. As the Dirichlet data on the outflow and alveoli, we also choose alveolar heliox, whereas we put humidified air on the inflow. The discrete system is very ill-conditioned due to the gas components taking zero values. In order for the solver to converge, we had to choose ε = 10 −3 . Furthermore, to avoid the singularity of the entropy density, we adjust the helium content in air and the nitrogen content in heliox to be 10 −7 , subtracting the same amount of water, in order to keep them summing to one. Note that this is not unreasonable, for example, the correct amount of helium in air is about 5.3 · 10 −7 . With these adjustments, the solver converges. The numerical results are shown in Figure 12. Both oxygen and carbon dioxide levels rise above the values in provided gas mixtures, before they start to decrease towards the equilibrium value. This is the expected behavior. However, the maximum values reached here are slightly lower than the ones found in [7]. This can be attributed to the perturbations of the zero concentrations and, as already seen, to the approximation of the geometry. 6. Conclusions. We have presented and analyzed a continuous space-time Galerkin method for cross-diffusion systems in entropy variable formulation, proving existence and convergence of discrete solutions, as well as existence of a weak solution of the continuous problem using the space-time approach. As opposed to time-stepping schemes, this approach provides an easy way to increase the approximation order simultaneously in space and time, makes hp-refinement in space-time possible, without the need for a globally fixed time-step size, and delivers numerical solutions, which can be evaluated at arbitrary points in time. Furthermore, at the same time, positivity and boundedness of the solutions are preserved also at the discrete level.
In the numerical examples, we have observed optimal convergence rates, given that the solution stays away from the singularities of the entropy. Lifting this restriction could be the topic of future works. Also, a more efficient numerical treatment of the space-time system is of interest, for example using a suitable preconditioner, a fine tuned solver, and making use of the mesh structure when using a tensor-product mesh.