{'fraction_non_alphanumeric': 0.058315374897822694, 'fraction_numerical': 0.0262131232815635, 'mean_word_length': 3.735550277117973, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 47, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient Numerical Algorithms based on Difference Potentials for Chemotaxis Systems in 3D

In this work, we propose efficient and accurate numerical algorithms based on Difference Potentials Method for numerical solution of chemotaxis systems and related models in 3D. The developed algorithms handle 3D irregular geometry with the use of only Cartesian meshes and employ Fast Poisson Solvers. In addition, to further enhance computational efficiency of the methods, we design a Difference-Potentials-based domain decomposition approach which allows mesh adaptivity and easy parallelization of the algorithm in space. Extensive numerical experiments are presented to illustrate the accuracy, efficiency and robustness of the developed numerical algorithms.

constant, γ ρ and γ c are the reaction coefficients. The parameter α is equal to either 1 or 0, which corresponds to the "parabolic-parabolic" or reduced "parabolic-elliptic" coupling, respectively. In this work, we will focus on the PKS chemotaxis model (1) with α = γ ρ = γ c = 1, but the developed methods are not restricted by this assumption and can be easily extended to more general chemotaxis systems and related models.
In the past several years, chemotaxis models have been extensively studied (see, e.g., [18,19,20,21,32] and references below). A known property of many chemotaxis systems is their ability to represent a concentration phenomenon that is mathematically described by fast growth of solutions in small neighborhoods of concentration points/curves. The solutions may blow up or may exhibit a very singular behavior. This blowup represents a mathematical description of a cell concentration phenomenon that occurs in real biological systems, see, e.g., [1,4,5,6,10,11,30,33].
Capturing blowing up or spiky solutions is a challenging task numerically, but at the same time design of robust numerical algorithms is crucial for the modeling and analysis of chemotaxis mechanisms. Let us briefly review some of the recent numerical methods that have been proposed for the "parabolic-parabolic" coupling of chemotaxis models in the literature. High-order discontinuous Galerkin methods have been proposed in [14,15] and [27] for chemotaxis models in 2D rectangular domains. A flux corrected finite element method is designed in 2D in [36], and is extended to chemotaxis models on stationary surface domains and cylindrical domains in [35,37]. A finite-volume based fractional step numerical method is proposed in [39] for models in 2D domains. However, the operator splitting approach may not be applicable when the convective part of the chemotaxis system is not hyperbolic. A simpler and more efficient second order positivity-preserving finite-volume central-upwind scheme is developed for 2D rectangular domains in [8], and is extended to fourth-order accurate numerical method in [9]. Finally, an upwind Difference Potentials Method is proposed in [12] to approximate chemotaxis models in 2D irregular domains, using uniform Cartesian meshes and Fast Poisson Solvers. Note that among the methods that have been proposed, only [37] is designed to handle chemotaxis models in 3D irregular domains by the use of unstructured meshes. For a more detailed review on recent developments of numerical methods for chemotaxis problems, the reader can consult [7] and [38].
In general, the design of numerical methods on unstructured meshes is more computationally intensive, in comparison to design of methods on structured meshes. Thus, in this work, we extend the numerical algorithm designed in [12] to chemotaxis systems in 3D irregular domain. The developed numerical methods based on Difference Potentials Method can handle 3D irregular geometry with the use of only Cartesian meshes and employ Fast Poisson Solvers. In addition, to further enhance computational efficiency of the numerical algorithms, we design a Difference-Potentials-based domain decomposition approach which allows mesh adaptivity and easy parallelization of the algorithm in space. The proposed numerical algorithms are shown to be positivity-preserving, robust, accurate and computationally efficient.
The paper is organized as follows. In Section 2, we formulate a positivity-preserving upwind Difference Potentials Method for accurate and efficient approximations of solutions to chemotaxis models in spherical domains. Next, in Section 3, a domain decomposition approach based on Difference Potentials Method is introduced, which improves the efficiency of the algorithm at no loss of accuracy. In Section 4, extensive numerical experiments (convergence studies in space and in time, long-time simulations, etc.) are presented to illustrate the accuracy, efficiency and robustness of the developed numerical algorithms.

An Algorithm Based on DPM
The current work is an extension of the work [12] in 2D, to 3D chemotaxis systems and to adaptive algorithms in space that allows easy parallelization. For the time being, we will consider the PKS chemotaxis model in a spherical domain, but the proposed methods can be extended to more general domains in 3D (and the main ideas of the algorithms stay the same). We employ a finite-volume-finite-difference scheme as the underlying discretization of the model (1) in space, combined with the idea of Difference Potentials Method ( [34] and some very recent work [3,12,13,28], etc.), that provides flexibility to handle irregular domains accurately and efficiently by the use of simple Cartesian meshes.
Introduction of the Auxiliary Domain: As a first step of the proposed method, we embed the domain Ω of model (1) into a computationally simple auxiliary domain Ω 0 ⊂ R 3 that we will select to be a cube. Next we discretize the auxiliary domain Ω 0 using a Cartesian mesh, with uniform cells D j,k,l = [x j− 1 2 , x j+ 1 2 ] × [y k− 1 2 , y k+ 1 2 ] × [z l− 1 2 , z l+ 1 2 ] of volume ∆x∆y∆z centered at the point (x j , y k , z l ), ( j, k, l = 1, . . . , N), and we assume here for simplicity, h := ∆x = ∆y = ∆z. Note that, we select the same auxiliary domain and the same mesh for the approximation of ρ and c, whereas in general, the auxiliary domains and meshes for ρ and c need not be the same. After that, we define the standard 7-point stencil with center placed at (x j , y k , z l ) that we will consider as a part of the discretization of model (1): j,k,l := {(x j , y k , z l ), (x j±1 , y k , z l ), (x j , y k±1 , z l ), (x j , y k , z l±1 )} Now we are ready to define point sets that will be used in the proposed hybrid finite-volume-finitedifference approximation combined with DPM.

Definition 2.1. Introduce following point sets:
• M 0 = (x j , y k , z l ) | (x j , y k , z l ) ∈ Ω 0 denotes the set of all the cell centers (x j , y k , z l ) that belong to the interior of the auxiliary domain Ω 0 ; ∈ Ω denotes the set of all the cell centers (x j , y k , z l ) that belong to the interior of the original domain Ω (see Fig. 1a); ∈ Ω 0 \Ω} is the set of all the cell centers (x j , y k , z l ) that are inside of the auxiliary domain Ω 0 , but belong to the exterior of the original domain Ω; The point sets N ± and N 0 are the sets of cell centers covered by the stencil N 7 j,k,l for every cell center (x j , y k , z l ) in M ± and M 0 respectively (see Fig. 1a); • γ = N + ∩ N − defines a thin layer of cell centers that straddles the continuous boundary Γ and is called the discrete grid boundary (see Fig. 1b); • γ in = M + ∩ γ and γ ex = M − ∩ γ are subsets of the discrete grid boundary that lie inside and outside of the spherical domain Ω respectively (see Fig. 1b).  (1). In each cell D j,k,l of volume |D j,k,l |, we denote the approximation of the cell average of the density, the approximation of point values of the density and the point values of the chemoattractant concentration at time t as:

Construction of the System of Discrete Equations for Model
We will also denote byρ i j,k,l , ρ i j,k,l and c i j,k,l the computed cell average of the density, point value of the density and the point value of the chemoattractant concentration at cell center (x j , y k , z l ) at the discrete time level t i , respectively.
Next, we use a hybrid finite-volume-finite-difference scheme in space, similar to 2D algorithms in [9,12] and first order implicit-explicit (IMEX) scheme in time as the underlying discretization for the Patlak-Keller-Segel model (1) at time t i+1 for every (x j , y k , z l ) in M + : where ∆ h is the discrete Laplace operator obtained using standard second-order centered finite difference approximation, I denotes the identity matrix of the same size with ∆ h and ∆t := t i+1 − t i . Here in the right hand side of density equation, g i j,k,l denotes the discretization of the convection term ∇ ·( χρ∇c) in (1) on , ∇ h y c i j,k±1/2,l , ∇ h z c i j,k,l±1/2 in the discrete gradient of c are computed similarly as in (7).
In addition, ρ i j±1/2,k,l , ρ i j,k±1/2,l , ρ i j,k,l±1/2 are the approximation of density values at the center of the six faces of the same cell D j,k,l , which are evaluated in an upwind manner. For example, the component in x-direction ρ i j+1/2,k,l is computed using the following piecewise linear construction: In the second order piecewise linear construction (9), x − x j , y − y k , z − z l denotes the vector defined by two points: (x, y, z) ∈ D j,k,l and (x j , y k , z l ) ∈ D j,k,l , and ∇ hρi j,k,l is the discrete gradient ofρ i at cell center (x j , y k , z l ) and time t i . Each component of the gradient ∇ hρi j,k,l is calculated using the minmod slope limiter. For example, the x-component ∇ h xρ i j,k,l is computed as where the minmod function is defined by The density values ρ i j− 1 2 ,k,l , ρ i j,k±1/2,l , ρ i j,k,l±1/2 at the center of the other faces in cell D j,k,l are computed similarly as in (8).
The Discrete Auxiliary Problem (AP). One of the important steps of DPM-based methods is the introduction of the auxiliary problem (AP). For brevity, we will denote u i+1 :=ρ i+1 or u i+1 := c i+1 . Thus, the two difference equations in (5), which are decoupled after discretization of model (1) with IMEX, can be cast into a compact form, i.e.
where L h,∆t := (I − ∆t∆ h ) and f i j,k,l is the right hand side function which is either f i j,k,l =ρ i j,k,l − ∆tg i j,k,l (for the density equation) or f i j,k,l = (1 − ∆t)c i j,k,l + ∆tρ i j,k,l (for the chemoattractant concentration equation), as in (5).
Next, we define the discrete Auxiliary Problem, which will play a key role in construction of the Particular Solution and the Difference Potentials as a part of DPM-based algorithm proposed in this work.

Definition 2.2.
At time t i+1 , given the grid function q i on M 0 , the following difference equations (12)- (13) are defined as the discrete Auxiliary Problem (AP): Here L h,∆t is the linear operator similar to the one in (11), but is defined now on a larger point set M 0 . (13) in the AP is chosen merely for efficiency of our algorithm, i.e. we employ Fast Poisson Solvers to solve the AP. In general, other boundary conditions can be selected for the AP as long as the defined AP is well-posed and can be solved computationally efficiently.

Construction of the Particular Solution.
Let us denote by G h,∆t f i j,k,l , (x j , y k , z l ) ∈ N + the Particular Solution which is defined on N + of the fully discrete problem (11) at time level t i+1 . The Particular Solution is obtained by solving the AP (12)-(13) with the following right hand side: and by restricting the computed solution from N 0 to N + .

Remark 2.
Note that, if the center (x j , y k , z l ) of cell D j,k,l belongs to γ in , some points (x j±1 , y k , z l ), (x j , y k±1 , z l ) or (x j , y k , z l±1 ) from stencil N 7 j,k,l with the center point (x j , y k , z l ) will lie in γ ex , and hence, outside of the domain Ω. To construct the Particular Solution for the density approximation in (5)-(6) and (12)- (14), one needs to compute the discretized convective term g i j,k,l for every point (x j , y k , z l ) ∈ M + . Hence, we need to approximate values of ρ and c at the points that belong to set γ ex . We will take the following strategy to define values of ρ and c on γ ex and γ in : • Initially, we approximate values ρ 0 j,k,l and c 0 j,k,l for (x j , y k , z l ) ∈ γ using 2-term extension operator (21) and zero Neumann boundary conditions, i.e. ρ 0 j,k,l ≈ ρ 0 (x j , y k , z l ) and c 0 j,k,l ≈ c 0 (x j , y k , z l ).
Here, (x j , y k , z l ) is the orthogonal projection on the continuous boundary Γ corresponding to a cell center (x j , y k , z l ) ∈ γ, and ρ 0 and c 0 are the initial conditions.
• At later time level t i+1 , g i j,k,l is computed using the solutionρ i j,k,l and c i j,k,l obtained from the discrete generalized Green's formula (26) at previous time level t i .

Construction of the Difference Potentials.
To construct the Difference Potentials, let us first define a linear space V γ of all grid functions v i+1 γ (x j , y k , z l ) at t i+1 on γ, and are extended by zero to other points in M 0 set. These grid functions v i+1 γ are called densities on the discrete grid boundary γ at the time level t i+1 .

Definition 2.3. The Difference Potential associated with a given density v i+1
γ ∈ V γ is the grid function P N + γ v i+1 γ defined on N + at time level t i+1 and is obtained by solving the AP (12)-(13) with the following right hand side: and by restricting the solution from N 0 to N + . In Definition 2.3, P N + γ denotes the operator which constructs Difference Potential P N + γ v i+1 γ from the density v i+1 γ at time t i+1 . Note that Difference Potential is a linear operator of the density function, is the index of the grid point in the set N + , n is the index of the grid point in the set γ, P N + γ v i+1 γ | m is the value of the Difference Potential at the grid point with index m and A nm are the coefficients of the Difference Potential.
Next we will introduce the trace operator. Given a grid function v i+1 defined on the point set N + , we denote by Tr γ v i+1 the trace or restriction of v i+1 from N + to the discrete grid boundary γ. Similarly, we define Tr γ i n v i+1 as the trace or restriction of v i+1 from N + to γ in ⊂ γ. We are ready to define an operator The operator P γ is a projection operator. Now we will state the key theorem in Difference Potentials Method, which allows us to reformulate the difference equation (11) defined in M + into an equivalent Boundary Equation with Projections (BEP) defined on the discrete grid boundary γ only (see [34] for more details).
where G h,∆t f i γ := Tr γ G h,∆t f i j,k,l is the trace of the Particular Solution restricted to the discrete grid boundary γ.

Remark 3.
Note, using that Difference Potential is a linear operator, we can recast (17) as where m is the index of the grid point in the set γ and G h,∆t f i m is the value of the Particular Solution at the grid point with index m in the set γ. Next we introduce the reduced BEP (19) defined only on γ in and show that it is equivalent to the BEP (17) defined on γ.
Moreover, the reduced BEP (19) contains only linearly independent equations.
Similarly to (17)- (18), the reduced BEP (19) can be recast as Remark 4. The BEP (17) or (19) reduces degrees of freedom from O(h −3 ) in the difference equation (11) to O(h −2 ). In addition, the reduced BEP (19) defined on γ in reduces the number of equations in BEP (17) by approximately one half, since |γ in | ≈ |γ|/2. Thus, using the reduced BEP (19) will further reduce the computational cost in our numerical algorithm and we will use the reduced BEP as a part of the proposed numerical algorithm.
Additionally, let us note the BEP (17) or BEP (19) will admit multiple solutions u i+1 γ since the system of equations (17) is equivalent to the system of difference equations (5) without imposing boundary conditions yet. Therefore, to construct a unique solution to BEP (19), we need to supply the BEP (19) with zero Neumann boundary conditions for the density and concentration. To impose boundary conditions efficiently into BEP, we will introduce the extension operator (21) below, similarly to [3,12,13,28], etc.
where n is the unit outward normal vector on Γ, d is the signed distance between a point (x j , y k , z l ) ∈ γ and the point of its orthogonal projection on the continuous boundary Γ in the direction of n. The parameter β controls the number of terms that will be used in the extension operator. If β = 0, we call (21) the 2-term extension operator and if β = 1, we call (21) the 3-term extension operator. We select β based on the regularity of the solution and to achieve the overall second-order accuracy of the numerical approximation in space.
Basically, the extension operator (21) defines values of the density u i+1 γ at the point of the discrete grid boundary (x j , y k , z l ) ∈ γ through the values of the continuous solution and its gradients at time t i+1 at the continuous boundary Γ of the domain with the desired accuracy.
Spectral Approach. With the extension operator (21) defined, the homogeneous Neumann boundary condition is readily incorporated, i.e. the second term in the right hand side of (21) will vanish. Hence, the extension operator (21) for our problem reduces to: Therefore, we incorporate the extension operator (22) into BEP (19) to solve for the unique density , (x j , y k , z l ) ∈ γ at time t i+1 as described below. However, to determine the unique density u i+1 γ at time t i+1 using the BEP (19) together with (22), we need to solve for the unknown solution u i+1 | Γ and its second order normal derivative ∂ 2 u i+1 ∂n 2 Γ at the continuous boundary Γ at t i+1 . To do this efficiently and accurately, we will employ spectral approximations for the two unknown terms on the boundary of the domain: where we take the basis functions φ ν (θ, ϕ) to be spherical harmonics and (θ, ϕ) are the polar and azimuthal angles respectively on the continuous boundary Γ. The spectral coefficients {C 0,i+1 . . , M) are the unknown coefficients that will be computed using BEP (19) at every time level t i+1 .
After we incorporate the spectral approximation (23) into the extension operator (22), we will have: Here, (θ, ϕ) are the polar and azimuthal angles for every point on the continuous boundary Γ which is the orthogonal projection of each point in the discrete grid boundary γ. Therefore, at every time level t i+1 , the BEP (19) becomes an over-determined linear system of dimension where Least Squares Method can be used to solve for the coefficients {C 0,i+1

Remark 5. (i)
In numerical experiments, we consider initial conditions for ρ and c in the form of f (−a(x 2 + y 2 +b(z−c) 2 )). Hence, we will only need the zonal modes in the spherical harmonics, i.e. φ ν (θ, ϕ) = P 0 ν (cos θ) where P 0 ν is the associated Legendre polynomial of degree ν and order 0. If the highest degree of the zonal spherical harmonics is L, the total number of harmonics for each term in (23) is only L + 1, which is significantly reduced comparing to (L + 1) 2 when full spectrum of spherical harmonics up to degree L is used.
(iii) In addition, we used the same number of harmonics M + 1 for the terms u i+1 | Γ and ∂ 2 u i+1 ∂n 2 Γ for simplicity in our implementation. In general, the number of harmonics can be chosen independently based on the regularity of each term. We should note that one can select a different basis functions φ ν (θ, ϕ) too, for example, spherical radial basis functions (see [22]).
Discrete Generalized Green's Formula. The final step of DPM is to use the computed density u i+1 γ to construct the approximation to continuous solutions of the chemotaxis model (1) in domain Ω subject to zero Neumann boundary conditions on ρ and c.

Proposition 2.4 (Discrete Generalized
Green's formula.). The discrete solution u i+1 :=ρ i+1 or u i+1 := c i+1 on N + constructed using Discrete Generalized Green's formula: is the approximation to the exact solution u := ρ or u := c, respectively, at (x j , y k , z l ) ∈ Ω at time t i+1 of the continuous chemotaxis model (1) subject to homogeneous Neumann boundary condition. We also conjecture that we have the following accuracy of the proposed numerical scheme: Remark 6. Indeed, in our numerical results (Section 4), we observe second order convergence in space and first order order convergence in time in the approximation of the solutions ρ and c. Moreover, we use ∆t = O(h 2 ) before blow-up for second order accuracy of the proposed algorithms. See also [3,28] for a more detailed discussion.
Positivity in a Spherical Domain. Now, we establish a positivity-preserving property of the proposed hybrid finite-volume-finite-difference DPM scheme for the approximation of solutions to (1).
Theorem 2.5 (Positivity in a Spherical Domain.). The discrete solution u i+1 =ρ i+1 or u i+1 = c i+1 obtained using the Discrete Generalized Green's formula (26) is non-negative on N + , provided that the following conditions are satisfied: where (x j , y k , z l ) ∈ M + and χ is the chemotactic sensitivity constant; 3. The density u i+1 γ is non-negative on the discrete grid boundary γ at time level t i+1 .
Proof. First, the discrete solution u i+1 =ρ i+1 or u i+1 = c i+1 is constructed by using the Discrete Generalized Green's formula (26). Also, we assume that the associated density u i+1 γ that was obtained by solving the BEP as discussed above is non-negative. From the construction of our algorithm above, it can be seen that u i+1 satisfies the discrete system on N 0 : where L h,∆t ≡ I − ∆t∆ h on M 0 . Hence, u i+1 satisfies the following discrete system defined on M + : where the right hand side f i is defined as the right hand side in (5) ( f i is either the right hand side for ρ or the right hand side for c). We will use discrete system (31) with the known density u i+1 γ to prove the non-negativity of the discrete solution u i+1 .
To this end, we will first establish the CFL-type condition (28) to guarantee the non-negativity of right hand side f i ≥ 0 in (31). When u i = c i , the right hand side f i in (31) is automatically non-negative provided Condition 1 in Theorem 3 is satisfied, since we assume that ∆t is sufficiently small, namely, ∆t < 1. When u i =ρ i , the right hand side f i in (31) is non-negative provided the CFL-type condition (28) is satisfied. To show this, we will proceed similarly to [9] and [12]. Let us first define the values at the center of six faces in a cell D j,k,l , (x j , y k , z l ) ∈ M + as: where the functionρ i is the piecewise linear reconstruction defined in (9). Note that ρ W j,k,l , ρ E j,k,l , ρ S j,k,l , ρ N j,k,l , ρ U j,k,l , ρ D j,k,l are non-negative due to Condition 1 of Theorem 2.5 and the positivity preserving property of the piecewise linear reconstruction (9). Further, the valueρ i j,k,l at the center of the cell D j,k,l can be expressed as: due to conservation property of the cell densities. Next, we rewrite right hand side f i j,k,l in (31) at point (x j , y k , z l ) ∈ M + as: where and F S j,k,l , F N j,k,l , F U j,k,l , F D j,k,l in the other directions are defined similarly. We will take the first term F E j,k,l as an example to establish the CFL-type condition (28): In this case, we require the coefficient of ρ E j,k,l to be non-negative, and it leads us to the constraint on ∆t: which will also ensure F E j,k,l ≥ 0. The details for non-negativity of F W j,k,l , F S j,k,l , F N j,k,l , F U j,k,l , F D j,k,l are similar to F E j,k,l . Thus, we obtain the CFL-type condition (28) to ensure non-negative right hand sides f i for the cell density ρ.
Thus far, we have shown that the right hand side f i j,k,l ≥ 0 in (31) is non-negative on M + under Conditions 1 and 2. What remains to be shown is that the solution u i+1 to the discrete system (31) is non-negative on M + , provided that the discrete boundary condition u i+1 γ e x is non-negative on the point set γ ex . Without loss of generality, let us assume that the center (x j , y k , z l ) of a cell D j,k,l belongs to M + , and only the point (x j−1 , y k , z l ) in the 7-point stencil of (x j , y k , z l ) is in γ ex . Since the boundary condition u i+1 γ e x to the discrete system (31) is given and non-negative, u i+1 j−1,k,l is non-negative for (x j−1 , y k , z l ) ∈ γ ex . Note that from the equation (31) for cell D j,k,l , we have Re-arranging the non-negative known value u i+1 j−1,k,l to the right hand side would give us: Now consider every point (x j , y k , z l ) in M + and we would have a modified version of (31): where∆ h is the modified Laplace operator as we discussed in (38). Note that the non-negative boundary condition u i+1 γ e x is already incorporated into the modified right hand sidef i . Also, the discrete equation (39) admits a unique solution, since (31) supplemented with condition u i+1 ≡ u i+1 γ e x on γ ex has a unique solution. Note that if the stencil N 7 j,k,l at point (x j , y k , z l ) ∈ M + has more than one point in γ ex , we simply move more boundary terms to the right hand side, like in (38). In other words, we always add non-negative boundary terms to the right hand side, so the modified right hand sidef i will always be non-negative if Conditions 1-3 of Theorem 2.5 are satisfied. Also, note that we only remove off-diagonal entries in ∆ h to obtain the modified Laplace operator∆ h , hence I − ∆t∆ h is an M-matrix with a non-negative inverse, similar to I − ∆t∆ h as discussed in [8,12,23]. Thus we conclude that the discrete system (31) admits a unique, non-negative solution.
Finally, since u i+1 (constructed using the Discrete Generalized Green's formula (26)) is a solution to the discrete system (31), we conclude thatρ i+1 and c i+1 are non-negative on N + , if Conditions 1-3 of Theorem 2.5 are satisfied.

Remark 7. Theorem 2.5 is not restricted to a spherical domain and can be applied to a general bounded
domain Ω ⊂ R 3 . In our algorithm, indeed when we require the density u i+1 γ to be non-negative, we observe no negative values in the discrete solutions, as can be seen in Fig. 9 and Fig. 11, for example. • Step 3: Inside the time loop at each time level t i+1 : update the time step ∆t using the minimum between CFL-type condition (28) and 0.5h 2 .

An Outline of Main
• Step 4a: At time level t i+1 : construct the Particular Solution G h,∆t f i using the updated ∆t on N + for the density ρ i+1 and concentration c i+1 respectively.
• Step 4b: At time level t i+1 : if the time step ∆t becomes smaller than the time step at previous time level t i , recompute the BEP matrix using the same time step ∆t as in Step 4a. Then, recompute the QR decomposition of the BEP matrix as in Step 2. Otherwise, skip this step.
• Step 6: At time level t i+1 : construct the Difference Potentials P N + γ u i+1 γ using the same time step ∆t as in Step 4a by solving the AP with the right hand side (16) and density u i+1 γ obtained in Step 5. • Step 7: At time level t i+1 : construct the discrete solution u i+1 (ρ i+1 or c i+1 ) via the Discrete Generalized Green's formula (26). • Step 8: Repeat Steps 3-7 till final time or the ∞-norm of the discrete solution ofρ i+1 is above certain threshold.
Step 4b is only required near blow-up, when time step decreases due to the CFL-type condition (28). Hence, the proposed algorithm based on Difference Potentials approach is very efficient for such problems.

A Domain Decomposition approach based on DPM
Note that the most computationally expensive step in our algorithm proposed in Section 2 for a single spherical domain is Step 4b, where we need to recompute the BEP (25) by solving the APs as many times as the total number of basis functions. In particular, we need to recompute BEP (25) at every time level near blow-up time in Step 4b, when time step ∆t decreases due to the CFL-type condition (28). In addition, as usual with any numerical algorithm, the computational cost of the entire algorithm increases significantly with global mesh refinement in 3D. However, the considered solutions (ρ and c) to chemotaxis model (1) have a compact support, which means global mesh refinement introduces unnecessary computation. Hence, to develop a more efficient scheme for chemotaxis models in 3D, we introduce the adaptivity in space that is compatible with Fast Poisson Solvers for the AP (12) and (13).
To this end, we consider a domain decomposition approach based on DPM. Our proposed domain decomposition algorithm follows and extends the numerical algorithms for interface problems in 2D [2,3,28]. In the domain decomposition approach, we decompose the spherical domain Ω into two non-intersecting sub-domains Ω 1 and Ω 2 . Next, similarly as in Section 2, we introduce computationally simple auxiliary domains Ω 0 (cubic) and embed each sub-domain Ω into its corresponding auxiliary domain ( = 1, 2). We assume the large values or the non-smooth parts of the solutions to chemotaxis system (1) near blow-up are located in sub-domain Ω 1 at any given time. Next, we discretize each auxiliary domain Ω 0 using uniform Cartesian meshes of dimension N × N × N and grid size h ( = 1, 2).  We denote the artificial interface between the two sub-domains as Z. In this work, we consider the following two types of interfaces: Z ∩ Γ = (Fig. 2) and Z ∩ Γ (Fig. 3), where Γ := ∂Ω is the continuous boundary of the original spherical domain Ω. Therefore, across the artificial interface Z, we impose the continuous interface conditions: where ∂ k u i+1 ∂n k denotes the k-th order normal derivative of u i+1 at time t i+1 , n denotes the unit normal vector on the interface Z, and u i+1 :=ρ i+1 or u i+1 := c i+1 in sub-domain Ω ( = 1, 2).
After we decompose the spherical domain Ω into two sub-domains, we proceed in each sub-domain as in the single domain approach (see Section 2). Again, we define the point sets N 0 , M 0 , N ± , M ± ( = 1, 2) as in Definition 2.1 for each Ω and Ω 0 ( = 1, 2). Next, we define the discrete grid boundary along the continuous boundary Γ and the discrete grid interface along the continuous interface Z in two cases respectively:

Case 1 of domain decomposition (suitable for the blow-up at the center of the domain):
The discrete grid boundary γ for sub-domain Ω 2 (Fig. 2b) is defined as the intersection of N + 2 and N − 2 within a small neighborhood of the boundary Γ. The discrete grid interfaces ζ are defined as the intersection of N + and N − ( = 1, 2) within a small neighborhood of the continuous interface Z for sub-domain Ω ( = 1, 2). See Fig. 2 for example of γ, ζ 1 , ζ 2 in cross-sectional view.

Case 2 of domain decomposition (suitable for the blow-up at the boundary of the domain): Z ∩ Γ
. In this case, we take the discrete grid boundary γ 1 for sub-domain Ω 1 to illustrate how to define the discrete grid boundary and the discrete grid interface. First, define the "discrete grid boundary" γ for the whole boundary Γ (similarly to Case 1). Next, include a point p ∈ γ into the discrete grid boundary γ 1 for sub-domain Ω 1 , if the polar angle of the orthogonal projection of p on the continuous boundary Γ belongs to [0, θ * + ]. Here, θ * is the polar angle at the intersection points of Γ and Z and the tolerance is introduced to include extra layer of points near the "wedge". The discrete grid boundary γ 2 and the discrete grid interfaces ζ 1 and ζ 2 are defined similarly. See Fig. 3 for example of γ 1 , γ 2 , ζ 1 , ζ 2 for geometry with wedges in cross-sectional view. BEPs for Case 1: Z ∩ Γ = . Note that ζ 1 is the "discrete grid boundary" for sub-domain Ω 1 and the union γ ∪ ζ 2 constitutes the "discrete grid boundary" for sub-domain Ω 2 (see Fig. 2). Next, we treat each sub-domain as a single domain and recall Theorem 2.1 to construct BEP (41) for sub-domain Ω 1 and BEP (42) for sub-domain Ω 2 : where the projection operators are defined as P ζ 1 := Tr ζ 1 P N + 1 ζ 1 and P γ∪ζ 2 := Tr γ∪ζ 2 P N + 2 (γ∪ζ 2 ) . As discussed in Section 2 and for efficiency of our algorithms, we reduce the above coupled BEPs (41) and (42) to the interior of discrete grid boundary γ in and the discrete grid interface ζ 1,in and ζ 2,in : Next, we follow the single domain approach in Section 2 and supplement the BEPs (43) and (44) with extension operators and spectral approximations, to incorporate the zero Neumann boundary condition and the continuous interface condition (40). Then we solve the coupled BEPs (43) and (44) and obtain the densities u i+1 ζ 1 , u i+1 γ and u i+1 ζ 2 . Finally, the discrete solution u i+1 in each sub-domain Ω ( = 1, 2) at time level t i+1 is constructed using the Discrete Generalized Green's formula (26) for each sub-domain.

BEPs for Case 2: Z ∩ Γ
. Similarly to Case 1, note that γ 1 and ζ 1 constitute the "discrete grid boundary" for sub-domain Ω 1 , and γ 2 and ζ 2 form the "discrete grid boundary" for sub-domain Ω 2 . Next, we treat each sub-domain as a single domain and recall Theorem 2.1 to formulate the BEP for sub-domain Ω : where u i+1 γ ∨ζ denotes the densities defined on both sets γ and ζ , and P γ ∨ζ is defined as: Remark 10. Note that any point p ∈ γ ∩ ζ contributes two values in the density u i+1 γ ∨ζ : one from density u i+1 γ and the other from density u i+1 ζ ( = 1, 2). Next, we follow [29] and take the average of the two values corresponding to point p and assign the average to be the "effective" density at point p. The "effective" density will be used in computing the Difference Potentials P N + (γ ∨ζ ) u γ ∨ζ . Also, note that we impose two equations at point p ∈ γ ∩ ζ in the BEPs (45): one corresponds to u i+1 γ and the other for u i+1 ζ ( = 1, 2). Again for efficiency of our algorithms, we reduce the coupled BEPs (45) to the interior of the discrete grid boundary γ ,in and interior of the discrete grid interface ζ ,in : on γ ,in and ζ ,in . Similarly, we follow the single domain approach in Section 2 and supplement the BEPs (47) with extension operators and spectral approximations, to incorporate the zero Neumann boundary condition and the continuous interface condition (40). After we solve the coupled BEPs (47) and obtain densities u i+1 γ 1 , u i+1 ζ 1 , u i+1 γ 2 and u i+1 ζ 2 , the discrete solution u i+1 in each sub-domain Ω ( = 1, 2) at time level t i+1 is constructed using the Discrete Generalized Green's formula (26) for each sub-domain.

Numerical Results
Throughout this section, we assume the single spherical domain Ω is centered at the origin and is of radius r = 0.5. Also, to illustrate the idea of domain decomposition approach, we assume the sub-domains Ω 1 and Ω 2 are non-intersecting parts of the spherical domain, and Ω = Ω 1 ∪ Ω 2 with r 1 = 0.25 for sub-domain Ω 1 and r 2 = 0.5 for sub-domain Ω 2 for both Cases 1 and 2 (see Section 3). Initial Conditions. We will study the following two examples of test problems to demonstrate the accuracy, efficiency and robustness of our proposed algorithms: • Initial conditions of ρ and c that will give blow-up at the center of the spherical domain: ρ 0 = 1000e −100(x 2 +y 2 +z 2 ) , c 0 = 500e −50(x 2 +y 2 +z 2 ) . (A) • Initial condition of ρ and c that will give blow-up at the boundary of the spherical domain:

Choice of Meshes.
We illustrate the choices of uniform Cartesian meshes using the single domain approach, and the domain decomposition approach follows similarly. In the single domain approach, we take the auxiliary domain to be a cubic domain: Here h := 2r/(N − 4), where r is the radius of the spherical domain and N denotes that the cubic auxiliary domain is discretized by mesh of dimension N × N × N, and N is the number of cells in x (or y, z) direction. Note that the choice of auxiliary domain here is only a few layers larger than the spherical domain Ω to have efficient computations.
Next we introduce the notations for meshes in the single domain (SD) and in the domain decomposition (DD) approaches: • "N" denotes that the mesh for Ω in the SD approach is dimension N × N × N cells, with grid size h = 2r/(N − 4); • "N 1 /N 2 " denotes that the meshes for sub-domains Ω 1 and Ω 2 in the DD approach are dimension N 1 × N 1 × N 1 cells and N 2 × N 2 × N 2 cells respectively, with h 1 = 2r 1 / (N 1 − 4) for the sub-domain Ω 1 and h 2 = 2r 2 /(N 2 − 4) for the sub-domain Ω 2 .
These choices of uniform Cartesian meshes are key to using Fast Poisson Solvers in 3D. We employed Fast Poisson Solver based on FFTW3 library [16] with openMP which enabled for better parallel efficiency. For other choices and comparison of high performance Fast Poisson Solvers, see [17].

Basis functions, extension operators and BEPs.
As mentioned in Section 2, the basis functions are selected to be the zonal spherical harmonics in the form of: We use a different number of zonal harmonics in the spectral approximations along the boundary and along the interface for Test (A) and Test (B). For Test (A), only 1 zonal harmonic is needed for each term in the 3-term extension operator along the boundary or the interface, due to: (i) the radial symmetry of the initial data and the solutions at later time levels; and (ii) the symmetry of the spherical domain. For Test (B), we use 20 harmonics for each term in the 3-term extension operator along the continuous interface Z, since we expect the solutions to be smooth along the interface Z at any time. However, we vary the number of harmonics used for each term in the 2-term extension operator along the boundary Γ in both the SD and DD approaches, and the number of harmonics used will be specified in the numerical results. The 2-term extension operators are used along the boundary for Test (B) since solution ρ loses regularity at the boundary due to blow-up at the boundary.
We should also note that all results from the SD approach are obtained using the reduced BEP (19) and results from the DD approach are obtained using the reduced BEPs (43)-(44) or (47) with = 1, 2.

Comparisons between single domain (SD) and domain decomposition (DD) approaches
In this subsection, we compare and contrast the numerical results obtained from the SD approach and the DD approach using Test (A) in the following aspects: • Convergence in space: Since there is no exact solution to the chemotaxis system (1), we use the discrete solutions on a finer mesh to construct the reference solutions. The ∞-norm error for the density ρ in space is computed as: whereρ i h * andρ i h * ( = 1, 2) are reference solutions with grid sizes h * and h * in the SD and DD approaches respectively, and the function 1 S denotes the characteristic function for a point set S: Here S = M + in the SD approach, or S = M + ( = 1, 2) in the DD approach. The ∞-norm error of the chemoattractant concentration c in space for the SD or DD approach is computed similarly.
• L 2 -norm of the relative error in time of the maximum value of the density ρ: where N T is the total number of time steps taken, and || · || ∞ is the maximum value of ρ in space at time t i .
• Convergence in time: We test convergence in time by fixing the meshes and refining the time steps: ∆t = ∆T/τ, where ∆T is a constant chosen to be 10 −7 and τ is a multiple of 2. Then, the error of the density ρ at the final time t i = 10 −6 is computed as: whereρ i ∆t * andρ i ,∆t * are the reference solutions at the same final time with the time step ∆t * in the SD and DD approaches respectively. The convergence in time for the concentration c at the final time t i = 10 −6 for the SD or the DD approach is computed similarly.
• Max density at time t i : which we will use to approximate the blow-up time.
• Second moment of ρ at time t i (see [24]): • Free energy (see [24,32]) of the chemotactic system (1) at time t i : In particular, we expect to see a decrease in free energy over time according to the second law of thermodynamics.

Comparison of the convergence at time 10 −6
In this subsection, we present the comparison of convergence studies both in space and in time for the SD and the DD approaches using Test (A). The final time for the convergence test is set to be 10 −6 so that (i) the solutions ρ and c are sufficiently smooth; and (ii) the reference solution on the finest mesh 516 × 516 × 516 in the SD approach can be obtained within a reasonable wall-clock time.

Remark 13.
Note that the single domain approach is more computationally expensive than the domain decomposition approach.  (49) and (50)) and the order of the convergence between SD and DD approaches to approximate ρ in Test (A) at time 10 −6 , with fixed time step ∆t = 10 −8 . The reference solutions in both SD and DD approaches are obtained using the discrete solution on mesh 516 × 516 × 516 in the SD approach.   (49) and (50)) and the order of the convergence between SD and DD approaches to approximate c in Test (A) at time 10 −6 , with fixed time step ∆t = 10 −8 . The reference solutions in both SD and DD approaches are obtained using the discrete solution on mesh 516 × 516 × 516 in the SD approach.

Convergence in Space.
Observe that in Table 1 and Table 2 the errors in sub-domain Ω 1 from both DD (h 1 /h 2 = 1/2) and DD (h 1 /h 2 = 1/4) are identical to those from using SD approach, if the grid spacing h and h 1 are the same. Also, the errors in sub-domain Ω 2 are less in magnitude than errors in sub-domain Ω 1 , even though the meshes for sub-domain Ω 2 are coarser than meshes for sub-domain Ω 1 . This shows that our strategy of using domain decomposition as in Fig. 2 is effective and the error in the entire domain Ω is dominated by the error from Ω 1 . In addition, the errors in sub-domain Ω 2 from using DD (h 1 /h 2 = 1/2) are smaller than using DD (h 1 /h 2 = 1/4), since DD (h 1 /h 2 = 1/2) employs finer mesh for sub-domain Ω 2 .
Overall, errors from sub-domain Ω 1 dominate the errors in Ω, and second order convergence in space are observed as expected in the SD approach, DD (h 1 /h 2 = 1/2) or DD (h 1 /h 2 = 1/4) for both ρ and c.

Convergence of max density in space.
In Table 3, we observe identical relative errors and second order convergence among SD, DD (h 1 /h 2 = 1/2) and DD (h 1 /h 2 = 1/4), which shows the robustness of the domain decomposition approach.  (52) and (53)) for Test (A). The results are obtained using fixed time step ∆t = 10 −8 to final time 10 −6 . The reference solutions for both SD and DD approaches are obtained using the discrete solution on mesh 516 × 516 × 516 in the SD approach.

Convergence in Time.
In Table 4 and Table 5, we observe first order convergence in time for SD approach and for both sub-domains in DD approach as we expected. The errors in sub-domain Ω 1 are similar to errors in the single domain Ω, regardless of the ratio h 1 /h 2 = 1/2 or h 1 /h 2 = 1/4. Although the meshes for Ω 2 are coarser than meshes for Ω 1 , the errors in sub-domain Ω 2 are much smaller than errors in sub-domain Ω 1 . We should note that errors from using SD approach and from using DD approach are similar, but are not exactly the same, because the mesh sizes are slightly different, i.e. h 1 /h 2 ≈ 1/2 as opposed to h 1 /h 2 = 1/2. Overall, the errors from sub-domain Ω 1 dominate the errors in the entire Ω as in the previous tests.  (54) and (55)) and the order of the convergence between SD and DD for ρ in Test (A) at the final time 10 −6 with time step ∆t = 10 −7 /τ. The underlying mesh for SD is N = 255 and the meshes for DD are N 1 /N 2 = 127/127 (h 1 /h 2 ≈ 1/2) and N 1 /N 2 = 127/63 (h 1 /h 2 ≈ 1/4). The reference solutions are obtained using time step ∆t = 10 −7 /512 at the same final time 10 −6 .   The reason for choosing the final time 6×10 −5 is that blow-up already occurs around T = 6×10 −5 (see Fig. 8).
We consider mesh size h for the SD approach and mesh sizes h 1 /h 2 for the DD approach in sub-domain Ω 1 and sub-domain Ω 2 respectively. Before blow-up, the time step is selected to be ∆t = 0.5h 2 for the SD approach and ∆t = min{0.5h 2 1 , 0.5h 2 2 } for the DD approach. Near blow-up, the time step is selected according to the CFL-type condition (28). Through Figs. 4-6, we obtain similar results in all aspects (|| ρ|| ∞ , free energy, etc.) between the SD and the DD approaches, as long as the finest grid sizes h and h 1 are close. DD approaches with ratios h 1 /h 2 ≈ 1/2 and h 1 /h 2 ≈ 1/4 give similar approximations, due to smooth solutions in Ω 2 . From numerical simulations on the CHPC clusters at the University of Utah, we observe that DD approach with h 1 /h 2 ≈ 1/4 gives approximately 8× speedup, while DD approach with h 1 /h 2 ≈ 1/2 gives approximately 4× speedup, in comparison to using the SD approach with similar accuracy in space. The speedup is mainly due to the reduction of degrees of freedom.

Numerical Results for Test (A) to final time 10 −4
As can be seen in the last two subsections, DD (h 1 /h 2 ≈ 1/4) gives similar accuracy to SD and DD (h 1 /h 2 ≈ 1/2), but is much more efficient. Hence we will continue with DD (h 1 /h 2 ≈ 1/4) for Test (A) to explore and simulate the chemotactic process on finer meshes to a longer time T = 10 −4 . This final time is already post blow-up, as can be seen in Fig. 8. Before blow-up, the time step is also selected to be  ∆t = 0.5h 2 for the SD approach and ∆t = min{0.5h 2 1 , 0.5h 2 2 } for the DD approach. Near blow-up, the time step is selected according to the CFL-type condition (28).
Before commenting on the longer-time simulation results, we should note that DD with mesh N 1 /N 2 = 255/127 gives similar results to SD with mesh N = 511, but SD with mesh N = 511 would take significantly longer time to finish as we observed in simulations.
We plot the evolution of || ρ|| ∞ versus time in Fig. 7a, and this test can be used to detect blow-up time, similar to [8,9,12]. In Fig. 7b and Fig. 7c, we show that the second moment and the free energy are both decreasing, which obey the second law of thermodynamics. In addition, the decrease rates of free energy are similar on different meshes before time step becomes restrictive due to the CFL-type condition (28) near blow-up time. After the time step becomes smaller, the decrease rates of the free energy are even larger, especially on finer meshes. This might be explained by the faster aggregations of cells on finer meshes.
In Fig. 8, snapshots of the cutline plots along z-axis of the solutions ρ at different times are given. We can see that the max value of ρ increases by two magnitudes from t = 3 × 10 −5 to t = 5 × 10 −5 . Density ρ already starts to develop singularity at the origin at t = 5 × 10 −5 . Eventually at t = 6 × 10 −5 , ρ develops a almost singular solution.
In Fig. 9, we present the 3D view of the isosurface plots of the solutions ρ and c at time t = 3 × 10 −5 , from which we saw that indeed the peak values of both solutions ρ and c occur at the origin only. We chose this time t = 3 × 10 −5 to give a better 3D visualization of the isosurface plots. As can be seen in Fig. 8, ρ at later times will be highly concentrated at the origin, which makes the 3D isosurface plots at later time levels uninformative. Moreover, the outlines of two cubic auxiliary domains (for Ω 1 and Ω 2 respectively) are also  given in Fig. 9, from which we can observe that the peak values indeed occur in sub-domain Ω 1 .

Selected Numerical Results for Test (B)
In this subsection, we present selected simulation results for Test (B). The blow-up is expected to occur at the north pole and our simulation results agree with our expectations. Again, the time step is selected to be ∆t = 0.5h 2 for the SD approach and ∆t = min{0.5h 2 1 , 0.5h 2 2 } for the DD approach before blow-up. Near blow-up, the time step is selected according to the CFL-type condition (28). In our simulations, the codes are terminated once the change in || ρ|| ∞ over two consecutive time steps exceeds 1000, i.e.
∆|| ρ|| ∞ ≥ 1000 (62) We use criterion (62) to determine that the solutions have reached blow-up. As can be seen in Fig. 10a, Fig. 10b and Table 6, the max density || ρ|| ∞ undergoes a drastic increase in a very small time interval. The second moment in Fig. 10c is increasing as the cells are moving away from the origin. The free energy in Fig. 10d are decreasing, which also obeys the second law of thermodynamics. We should note  that the free energies highly agree on difference meshes. This is due to the relatively small value of c (in the magnitude of 10). Overall, we observe in Fig. 10 that the SD approach, DD (h 1 /h 2 ≈ 1/2) and DD (h 1 /h 2 ≈ 1/4) approaches give similar results for Test (B) as well. Finally, in Fig. 11, we present the 3D view of isosurface plots of ρ and c for Test (B) at different times: • At time t = 0.03, we observe that both solutions ρ and c are in motion towards the north pole and continuity of solutions are also observed across the interface.
• At time t = 0.04, we observe that the peak value of c already occurs at the north pole, while the peak value of ρ is still in motion towards the north pole. This unveils the role that the chemoattractant plays and might help to explain why the blow-up solution of ρ occurs at the north pole eventually.
• At time t ≈ 0.079812, indeed we observe an almost singular solution at the north pole only. For both ρ and c, the values are stratified away from the point of aggregation, i.e. the north pole. Moreover, no negative values are observed. We should also note that the simulations can proceed to a longer final time by choosing a different stopping criterion, but the simulation will not give a significant difference (see Fig. 12 for simulation results at t ≈ 0.079846).
The two cubic domains with white outlines in Fig. 11 and Fig. 12 are auxiliary domains for sub-domain Ω 1 and sub-domain Ω 2 respectively.

Proof to Theorem 2.1
Proof. For the reader's convenience, let us recall proof of Theorem 2.1 in [12,34]. First of all, let us assume that there is a density u i+1 γ that satisfies the BEP (17). Using u i+1 γ , we can construct a superposition of the Particular Solution using (12)(13)(14) and the Difference Potential using (12)(13) and (16): j,k,l = P N + γ u i+1 γ + G h,∆t f i , (x j , y k , z l ) ∈ M 0 .
It can be seen that u i+1 defined in (63) is a solution to the following discrete system on N 0 : L h,∆t u i+1 j,k,l = f i j,k,l , (x j , y k , z l ) ∈ M + , L h,∆t [u i+1 γ ], (x j , y k , z l ) ∈ M − , u i+1 j,k,l = 0, (x j , y k , z l ) ∈ N 0 \M 0 , which implies that u i+1 j,k,l , (x i , y j , z k ) ∈ M + in (63) is some solution to (11) on M + , and note that the trace of u i+1 satisfies Tr γ u i+1 = Tr γ P N + γ u i+1 from the definition of trace operator. Moreover, from our assumption on the density u i+1 γ , we have P γ u i+1 γ + G h,∆t f i γ = u i+1 γ . Thus, Tr γ u i+1 = u i+1 γ and we can conclude that if some density u i+1 γ satisfies the BEP (17), it is a trace to some solution of system (11). Now assume the density u i+1 γ is a trace to some solution of system (11): u i+1 γ = Tr γ u i+1 where u i+1 j,k,l , (x j , y k , z l ) ∈ N + is the unique solution to the difference equation (11) subject to certain boundary condition. Then by our assumption on u i+1 , u i+1 with zero extension from N + to N 0 is also a solution to the system of difference equation (64) on N 0 . Since P N + γ u i+1 γ + G h,∆t f i is also a solution the system of difference equation (64), we conclude that:  by the uniqueness argument. Finally, restriction of the u i+1 from N + to the discrete grid boundary γ would give us and u i+1 γ = P γ u i+1 γ + G h,∆t f i γ , since u i+1 γ = Tr γ u i+1 by our assumption. In other words, if u i+1 γ is a trace to some solution of system (11), it satisfies the BEP (17).

Proof to Proposition 2.2
Proof. The proof is similar to ideas in [12,34]. Without boundary conditions, the difference equations (11) would admit infinite number of solutions and so will the BEP (17), due to Theorem 2.1. However, if we provide density values u i+1 γ e x on γ ex , the difference equation (11) will be well-posed and will admit a unique solution. Hence, the BEP (17) will have a unique solution if u i+1 γ e x is given. The solution to BEP (17) is uniquely determined by densities on γ ex , thus the solution u i+1 γ to BEP (17) has dimension |γ ex |, which is the cardinality of set γ ex . As a result, the BEP (17) has rank |γ| − |γ ex | = |γ in |.

Proof to Theorem 2.3
Proof. The proof is based on ideas in [12,34]. First, define the grid function: where P i+1 is a solution to the AP (12)-(13) on N 0 with right hand side (16) using density u i+1 γ , G i+1 is a solution to the AP (12)-(13) on N 0 with right hand side (14), and u i+1 γ is extended from γ to N 0 by zero. By the construction of v i+1 , one can see that v i+1 is a solution to the following difference equation: Thus, we conclude that v i+1 solves the following homogeneous difference equations on M − : In addition, due to the construction of v i+1 , P i+1 and G i+1 , the grid function v i+1 satisfies the following boundary condition: v i+1 = 0, on N 0 \M 0 .
Hence, it is enough to show that (72) is equivalent to (73) to prove the equivalence between the BEP (17) and the reduced BEP (19). First, note that if (72) is true, then (73) is obviously satisfied. Now, assume that (73) is true and let us show that (72) holds. Let us consider problem (70): L h,∆t v i+1 = 0 on M − , subject to boundary conditions (71) and (73), since the set γ in ∪ (N 0 \M 0 ) is the boundary set for set M − . Then we have the following discrete boundary value problem: which admits a unique zero solution: v i+1 = 0 on M − . Since γ ex ⊂ M − , we conclude that v i+1 = 0 on γ ex , as well as on γ ≡ γ ex ∪ γ in , which shows that (73) implies (72). Thus, we showed that (72) is equivalent to (73), and therefore, BEP (17) is equivalent to the reduced BEP (19). Moreover, due to Proposition 2.2, the reduced BEP (19) consists of only linearly independent equations.