1. Introduction
The development of multicellular organisms, cell migration and information processing by nerve cells in the brain depends on various interactions between the cells and other biological factors. In these phenomena, the functions of these interactions determine the state of the subsequent time evolution and exhibit the occurrence and various behaviour of the patterns. When modelling such phenomena, there are cases in which long-range interactions that affect distant objects globally in space naturally appear. These interactions are called nonlocal interactions. They have attracted considerable attention in various fields and have been studied extensively. The existence of nonlocal interactions has been experimentally suggested, for example, in phenomena such as neural firing in the brain [Reference Kuffler13], pigmentation patterns in animal skin [Reference Nakamasu, Takahashi, Kanbe and Kondo18, Reference Watanabe and Kondo22, Reference Yamanaka and Kondo23], development of multicellular organisms, cell migration and adhesion [Reference Katsunuma, Honda and Shinoda11].
In the experiment by Kuffler [Reference Kuffler13], the electrode was set at a ganglion cell in the receptive field of the retina of a cat. The firing rate against the light stimulus was measured by illuminating two points at different distances from the electrodes. Ganglion cells respond to light stimuli locally in space and, conversely, inhibit them laterally in space. From these observations, by considering the local excitation and lateral inhibition as positive and negative values, respectively, this interaction can be modelled by a sign-change function with radial symmetry. This function is called the local activation and lateral inhibition (LALI) interaction or Mexican hat.
Experimental results on the interactions between yellow and black pigment cells in the skin of zebrafish were reported by Nakamasu et al. [Reference Nakamasu, Takahashi, Kanbe and Kondo18]. During these experiments, a square section of black pigmented cells within a zebrafish’s skin stripe was eliminated using laser ablation. The regions of yellow pigmented cells surrounding the squares were removed by laser ablation. Several zebrafish were prepared with different patterns of yellow pigmented cells, which were removed using a laser. For each pattern of yellow pigment cells, the number of black pigment cells proliferating in the square located at the centre of the pattern was quantified for two weeks. The comparison of the number of proliferating black pigment cells revealed the existence of both long- and short-range interactions between these cells. Moreover, the derived interactions are summarised as a network. A theoretical method to reduce a given reaction-diffusion network with spatial interactions, such as metabolites or signals with arbitrary factors, into the shape of an essential integral kernel has been proposed by Ei et al. [Reference Ei, Ishii, Kondo, Miura and Tanaka9]. A Mexican hat function was theoretically derived by applying this reduction method to the network given by Nakamasu et al. [Reference Nakamasu, Takahashi, Kanbe and Kondo18]. Additionally, cells such as the pigment cells in the zebrafish extend their cellular projections to exchange the biological signals with each other. Here, we refer to the papers by Katsunuma et al. [Reference Katsunuma, Honda and Shinoda11] and Kondo [Reference Kondo12]. Hamada et al. [Reference Hamada, Watanabe and Lau14] and Watanabe and Kondo [Reference Watanabe and Kondo22] reported that pigment cells send different biological signals depending on the length of their cellular projections.
Katsunuma et al. [Reference Katsunuma, Honda and Shinoda11] investigated the behaviour of cell adhesion by using two types of cellular adhesion molecules in HEK 293 cells. These cells also have cellular projections that are ten times longer than their body size. Cells can sense the cell density around them using their total body from the tip of the leading edge, and they can decide the directions of cell migration and cell adhesion.
Based on this biological background, numerous mathematical models have been proposed and analysed. The nonlocal interactions are often modelled by convolution with a suitable integral kernel. Here, we introduce two types of the model of the nonlocal interaction. We call the convolution itself without a derivative the normal type, and the gradient of convolution in the advection term the advective type. We introduce mathematical models with a normal type of nonlocal interaction. The peaks of biological signals are located distantly from the centre of the cell body during signal transduction through cellular projections. Here, we refer to the observational results of Hamada et al. [Reference Hamada, Watanabe and Lau14] and Watanabe and Kondo [Reference Watanabe and Kondo22]. For these interactions, many models that impose a convolution with an integral kernel with peaks distant from the origin have been proposed. For example, as models of pigmentation patterns in animal skin [Reference Ninomiya, Tanaka and Yamamoto19, Reference Ninomiya, Tanaka and Yamamoto20], population dispersal for biological individuals [Reference Ei, Guo, Ishii and Wu8, Reference Hutson, Martinez, Mischaikow and Vickers15], and vegetation patterns [Reference Alfaro, Izuhara and Mimura1], the following nonlocal evolution equation is proposed:

where
$d \ge 0$
is a constant,
$u=u(x,t)$
is the density,
$f$
is a suitable reaction or growth term,
$K$
is an integral kernel and
$*$
denotes the convolution of two functions in the space variable:

Analytical results were reported by Bates et al. [Reference Bates, Fife, Ren and Wang4], Coville et al. [Reference Coville, Dávila and Martíanez7] and Ei et al. [Reference Ei, Guo, Ishii and Wu8]. It is rigorously shown by Ninomiya et al. [Reference Ninomiya, Tanaka and Yamamoto19] that the nonlocal interaction plays a role to induce the Turing instability. By imposing a convolution using the Heaviside function, a nonlocal evolution equation to investigate the dynamics of the membrane potential of neurones in the brain was proposed by Amari [Reference Amari2]. Additionally, motivated by the pattern formations observed in animal skins, a nonlocal model by applying the cut function to the convolution term was proposed by Kondo [Reference Kondo12]. This model can reproduce various patterns by changing only the kernel shape, even though it comprises only one component. The above nonlocal interactions can be derived from the continuation of spatially discretised models with intercellular interactions by Ei et al. [Reference Ei, Ishii, Sato, Tanaka, Wang and Yasugi10].
Next, we introduce mathematical models of the advective type of nonlocal interaction. As a first example, the aggregation-diffusion equation was proposed and analysed for cell migration and collective motion by Bailo et al. [Reference Bailo, Carrillo, Murakawa and Schmidtchen3] and Carrillo et al. [Reference Carrillo, Craig and Yao5]:

where
$\rho =\rho (x,t)$
denotes the cell density at position
$x$
at time
$t\gt 0$
and
$m\ge 1$
is a constant. If the potential
$W$
is positive, the convolution term
$\nabla ( W* \rho )$
determines the velocity of the advection by integrating the gradient of
$\rho$
. The density
$\rho$
at each point is advected towards the gradient of
$\rho$
. Thus, the second term provides the aggregation effect. If
$W$
is positive with a compact support, then the compact support corresponds to the total cell body. Subsequently, the term
$\nabla ( W* \rho )$
provides the effect that determines the velocity of the advection by sensing the cell density gradient in the total cell body. When
$m=1$
, this model can be classified as a nonlocal Fokker–Planck equation, whereas when
$m\gt 1$
, it can be classified as a nonlocal porous medium equation.
Another example is the cell adhesion model. It has been proposed to describe and analyse the cell adhesion phenomena by Carrillo et al. [Reference Carrillo, Murakawa, Sato, Togashi and Trush6]:

In contrast to the aggregation-diffusion equation, the velocity in the advection term is saturated by the cell density
$(1-\rho )$
in this model. Carrillo et al. [Reference Carrillo, Murakawa, Sato, Togashi and Trush6] reported that this model can replicate cell adhesion and cell sorting phenomena both qualitatively and quantitatively. For cell migration and cell adhesion processes, the integral kernel is called the potential, and the Mexican hat (or LALI) function is crucial in the local attraction and lateral repulsion in these two models.
Cell pattern formation, migration, and adhesion play pivotal roles in the biological development of various organs and tissues. Therefore, revealing the mechanisms underlying these phenomena is an important problem. However, the nonlocal term with the normal or advection type in nonlocal equations occasionally makes analysis difficult. To overcome this difficulty, the approximation of a nonlocal term by another type of term can be a solution. In light of this, we aim to reveal whether advective nonlocal interactions can be approximated by local dynamics. As a first step, we propose an approximation method for advective nonlocal terms in the nonlocal Fokker–Planck equation using a Keller–Segel system with multiple auxiliary chemotactic factors. Although the above models of (1.1) and (1.2) are considered in higher-dimensional spaces, we consider the problem on one-dimensional bounded domain in the simple case to make the regularity of the solution and fundamental solution to the elliptic equation in quasi-steady state easier to deal with. The nonlocal Fokker–Planck equation and Keller–Segel systems are basic models with advective nonlocal interactions and typical local dynamics, respectively. We show that any smooth kernel can be expanded by combining the fundamental solutions for an elliptic equation in the Keller–Segel system. Furthermore, we report that the solution to the nonlocal Fokker–Planck equation with an even smooth kernel can be approximated by that of the Keller–Segel system with specified parameters depending on the integral kernel shape. The equations used to approximate nonlocal interactions are reaction-diffusion equations, and this idea is the same as that of Ninomiya et al. [Reference Ninomiya, Tanaka and Yamamoto19, Reference Ninomiya, Tanaka and Yamamoto20]. In these papers, the results of approximation by the singular limit of the reaction-diffusion system for nonlocal interactions of normal type in one-dimensional space are reported. This is based on a theory to approximate the integral kernel by a linear sum of fundamental solutions to an elliptic equation of second order. In contrast, for the approximation of advective nonlocal interactions, it is necessary to approximate the derivative of the integral kernel by a linear sum of the aforementioned fundamental solutions, as in Corollary 2.6 below. The
$C^1$
estimate is the essential difference and difficulty between this approximation problem and the previous studies. To prove this convergence, we provide an estimation of the approximation using the Lagrange interpolation polynomial with the Chebyshev nodes. This enables us to evaluate the relationship between the error accuracy and the number of equations in the approximation.
The remainder of this paper is organised as follows: Section 2 outlines the mathematical framework and summarises the main results. In Section 3, the existence theorem is established, followed by a singular limit analysis detailed in Section 4. Section 5 starts by presenting a precise formula for the Lagrange interpolation polynomial’s coefficient. This is followed by a method to ascertain the coefficient of the linear sum in the fundamental solution for an elliptic equation, characterised by the shape of its integral kernel. The section also includes a proof of the fundamental solution’s series expansion. A linear stability analysis is then conducted in Section 6. The paper concludes with Section 7, summarising the study’s findings and implications.
2. Mathematical settings and main results
In this section, we describe the mathematical settings and results. We denote the theoretical concentration or cell density at position
$x \in \Omega \,:\!=\, [-L , L]$
at time
$t\gt 0$
by
$\rho = \rho (x,t)$
. We investigate the solution to the following nonlocal Fokker–Planck equation:

where the periodic boundary condition

is imposed and the initial datum is given by
$\rho (x, 0) \,:\!=\, \rho _0(x)$
. Here,
$W*u(x)$
is defined by

for
$W \in L^1_{\mathrm{per}}(\Omega ) \,:\!=\, \{u|_\Omega \in L^1(\Omega ) \ | \ u(x)=u(x+2L), \ x \in {\mathbb{R}} \}$
. Setting
$d_j\gt 0, (j=1,\ldots , M)$
which are the diffusion coefficients in (
$\mbox{KS}^{M,\varepsilon }$
) introduced below, we define the following function

This is actually a fundamental solution to the elliptic equation explained in Lemma 4.1 below. The typical examples of
$W$
are as follows:




with any
$a_1, a_2\gt 0$
, where
$R_1\gt R_0\gt 0$
are constants called the sensing radius,
$B(R_0)$
is a ball with radius
$R_0$
and origin centre, and

The profiles of (2.6) and (2.7) are presented in Figures 1 (a) and 2 (a), respectively. The nonlocal Fokker–Planck equation (P) with the integral kernels (2.5) and (2.6) corresponds to the parabolic-elliptic Keller–Segel systems. A linear stability analysis is presented in Section 6. Integral kernels (2.7) and (2.8) were introduced by Carrillo et al. [Reference Carrillo, Murakawa, Sato, Togashi and Trush6] and Murakawa and Togashi [Reference Murakawa and Togashi17]. They have a compact support corresponding to the cell body. If these integral kernels are used for the potential
$W$
in (P), this describes the situation in which
$\rho$
at each point detects the surrounding cell density in the own cell body and the velocity of the aggregation is determined. This corresponds to the Haptotaxis phenomenon.
First, we have the following existence result. To construct a mild solution to (P), we define a function space with a norm as

for any time
$\tau \gt 0$
. Introducing the following function:

where
$i$
is the imaginary unit and

we define the map
$\Gamma \,{:}\, E_{\tau } \to E_{\tau }$
as

We say that a function
$u \in C([0,T]; H^1(\Omega ))$
for any
$T\gt 0$
is a mild solution to (P) with
$\rho (x,0)=\rho _0(x)$
, provided
$ u=\Gamma [u]$
. The following proposition is proven using the standard argument of the fixed point theorem.
Proposition 2.1.
Let
$R\gt 0$
be an arbitrary real number, and assume that
$W \in W_{\mathrm{per}}^{1,1}(\Omega )\,:\!=\,\{u|_\Omega \in W^{1,1}(\Omega ) \ | \ u(x)=u(x+2L), \ x \in {\mathbb{R}} \}$
and

Then, for any
$T\gt 0$
, there exists a unique mild solution
$\rho$
to (
P
) in
$ C( [0, T]; H^1(\Omega ) ) \cap L^2( 0, T; H^2(\Omega ) )$
satisfying

where
$C_0 =C_0(L,R,T, \|W_x\|_{L^1(\Omega )} )$
. Moreover, this mild solution satisfies (P) in
$L^2(0,T;L^2(\Omega ))$
.
Next, we will approximate the solution to (P) with any integral kernel using that to a Keller–Segel system which is a local dynamics. Introducing the auxiliary factors
$v_j^{M,\varepsilon } = v_j^{M,\varepsilon }(x,t), \ (j=1, \ldots , M)$
, we consider the following Keller–Segel system in which the linear sum of
$v_j^{M,\varepsilon }$
is imposed in nonlocal term in (P):

Here,
$0 \lt \varepsilon \ll 1$
is a sufficiently small parameter,
$d_j \gt 0$
is the diffusion coefficient, and each
$a_j \in {\mathbb{R}}$
is a constant that determines whether
$v_j^{M,\varepsilon }$
acts as an attractive or repulsive substance in the aggregation process of
$\rho ^{M,\varepsilon }$
. Because the solutions to (
$\mbox{KS}^{M,\varepsilon }$
) depends on
$M$
and
$\varepsilon$
, we denoted them by
$(\rho ^{M,\varepsilon }, v_j^{M,\varepsilon })$
, respectively. The same periodic boundary condition as that in (P) is imposed in the equations of (
$\mbox{KS}^{M,\varepsilon }$
) as follows:

for
$j=1, \ldots , M$
. Furthermore, we impose the following initial conditions as

(
$\mbox{KS}^{M,\varepsilon }$
) is a Keller–Segel system with multiple components with the linear sensitivity function. The role of
$v_j^{M,\varepsilon }$
can be distinguished by the sign of the coefficient
$a_j$
. If
$a_j\gt 0$
, then
$v_j^{M,\varepsilon }$
is an attractive substance for
$\rho ^{M,\varepsilon }$
and
$\rho ^{M,\varepsilon }$
aggregates toward to the region in which the gradient of
$v_j^{M,\varepsilon }$
is high independently on the value of its concentration. In contrast, if
$a_j\lt 0$
,
$v_j^{M,\varepsilon }$
acts as a repulsive substance for
$\rho ^{M,\varepsilon }$
, and
$\rho ^{M,\varepsilon }$
migrates away from the region in which the gradient of
$v_j^{M,\varepsilon }$
is high.
Introducing the following function:

we define the maps
$\Psi _j: E_{\tau } \to E_{\tau }$
and
$\Phi : E_{\tau } \to E_{\tau }$
as


respectively. We now define the scalar integral equation of (
$\mbox{KS}^{M,\varepsilon }$
) as
$\rho ^{M,\varepsilon } = \Phi [\rho ^{M,\varepsilon }]$
. We say that a function
$(u,\{v_j\}_j) \in C([0,T]; H^1(\Omega ))\times C( [0, T]; C^2(\Omega ) )$
for any
$T\gt 0$
is a mild solution to (
$\mbox{KS}^{M,\varepsilon }$
) with (2.11), provided
$ u=\Phi [u]$
and
$v_j=\Psi _j[u]$
for
$j=1,\ldots ,M$
.
With the above settings, we obtain a unique solution to (
$\mbox{KS}^{M,\varepsilon }$
) using the fixed point argument.
Theorem 2.2.
Let
$R\gt 0$
be an arbitrary real number and assume


Then, there exists a real positive number
$\tau _0=\tau _0(M, \{a_j,d_j\}_{j=1}^M, L, R)$
such that for any
$\varepsilon \gt 0$
there exists a unique mild solution
$(\rho ^{M,\varepsilon }, \{ v_j^{M,\varepsilon } \}_{j=1}^M)$
to (
$\mbox{KS}^{M,\varepsilon }$
) in
$ C( [0, \tau _0]; H^1(\Omega ) ) \times C( [0, \tau _0]; C^2(\Omega ) )$
satisfying

We can extend the existence time of the solution to an arbitrary time
$T$
as follows:
Corollary 2.3.
Suppose the same assumption of Theorem
2.2
. For any
$T\gt 0$
, there exists a positive constant
$\tilde {C}_0=\tilde {C}_0(\tau _0,T)$
such that there exists a unique mild solution
$(\rho ^{M,\varepsilon }, \{v_j^{M,\varepsilon }\}_{j=1}^M)$
to (
$\mbox{KS}^{M,\varepsilon }$
) in
$ C( [0, T]; H^1(\Omega ) ) \cap L^2( 0, T; H^2(\Omega ) ) \cap H^1(0,T;L^2(\Omega )) \times C( [0, T]; C^2(\Omega ) )\cap L^2( 0, T; H^3(\Omega ) )\cap H^1( 0, T; L^2(\Omega ) )$
satisfying

Moreover, this mild solution
$(\rho ^{M,\varepsilon }, \{v_j^{M,\varepsilon }\}_{j=1}^M)$
satisfies the system (
$\mbox{KS}^{M,\varepsilon }$
) in
$L^2(0,T;L^2(\Omega ))\times L^2(0,T;C(\Omega ))$
.
In order to show the relationship between the solutions of nonlocal Fokker–Planck equation (P) with any potential
$W$
and the Keller–Segel system (
$\mbox{KS}^{M,\varepsilon }$
), we first investigate the relationship between the solution to (P) with the potential provided by
$W= \sum _{j=1}^M a_j k_j$
and the solution to (
$\mbox{KS}^{M,\varepsilon }$
).
Theorem 2.4.
Let
$M$
be an arbitrary fixed natural number, and
$\rho$
be a solution to (P) equipped with
$W = \sum _{j=1}^M a_j k_j$
and the initial value
$\rho _0 \in C^2(\Omega )$
. Let
$\rho ^{M,\varepsilon }$
be a solution to (
$\mbox{KS}^{M,\varepsilon }$
) equipped with


Then, for any
$T\gt 0$
, there exist positive constants
$C_1$
and
$C_2$
that depend on
$M$
,
$\{a_j,d_j\}_{j=1}^M$
,
$L$
,
$R$
and
$T$
such that for any
$\varepsilon \gt 0$


Remark 2.5.
We note that
$C_1$
and
$ C_2$
are independent of
$\varepsilon$
. Moreover, using the Sobolev embedding theorem for the first terms on the left-hand sides of (2.19) and (2.20) implies that

where
$\tilde {C}_1$
and
$\tilde {C}_2$
are obtained by multiplying the constants of the Sobolev embedding theorem by
$C_1$
and
$C_2$
, respectively.
The first convergence in Theorem 2.4 shows not only that the solution
$\rho ^{M,\varepsilon }$
to (
$\mbox{KS}^{M,\varepsilon }$
) is sufficiently close to that to (P) with
$W= \sum _{j=1}^M a_j k_j$
when
$\varepsilon$
is very small, but also that the convergence rate is of the order of
$\varepsilon$
. The second convergence shows that the solution of the auxiliary substances
$v_j^{M,\varepsilon }$
is also extremely close to
$k_j*\rho$
as
$\varepsilon$
tends to
$0$
. The proof of this theorem is presented in Section 4.
Using the convergence of Theorem 2.4, we can approximate the solution to (P) with any even smooth kernel
$W$
as that to (
$\mbox{KS}^{M,\varepsilon }$
) with the specified parameters
$\{d_j, a_j \}_{j=1}^M$
. Indeed, the parameters
$\{ a_j\}_{j=1}^M$
are determined by the shape of
$W$
using Theorem 5.3 in Section 5. Using the interpolation polynomial with the Chebyshev nodes, we can demonstrate the convergence in Theorem 5.3. The explicit formula for the coefficients in the Lagrange interpolation polynomial for an arbitrarily given function is constructed in Proposition 5.2 for the proof of Theorem 5.3. Although Theorem 5.3 and Proposition 5.2 are some of our main results, they are presented in Section 5 for convenience. Setting the diffusion coefficients of
$v_j^{M,\varepsilon }$
as

and
$d_1$
may be limit of infinity, we obtain

The choice of
$d_j$
is same as that in [Reference Ninomiya, Tanaka and Yamamoto19]. Because the profile of the fundamental solution
$k_j$
is unimodal even if the value of
$j$
changes, it seems to be difficult to approximate any potential
$W$
by the linear sum of
$k_j$
. However, we can obtain the following corollary. Let
$f$
be

and we set the assumption that for
$n \in \mathbb{N}$

and there exists a positive constant independent of
$n$
such that

for any
$n \in \mathbb{N}$
. We note that the function
$ L- \cosh ^{-1}(x)$
is the inverse function of
$\cosh (L-x)$
on
$[0,L]$
. Since
$(L- \cosh ^{-1}(x))'=-1/\sqrt {x^2-1}\,=\!:\,\mathscr{F}$
and it satisfies the recurrence formula
$(x^2-1)\mathscr{F}^{(n+2)} = -x(2n+1)\mathscr{F}^{(n+1)} -n^2\mathscr{F}^{(n)}$
, we see that
$f \in C^{n+1}((1, \cosh L])$
. Thus, from the assumption (2.23), we obtain that
$f \in C^{n+1}([1, \cosh L])$
.
Corollary 2.6.
Assume that
$W \in C^\infty (\Omega )$
,
$W$
is even, (2.23) and (2.24) for any
$n\in \mathbb{N}$
. Then, for any
$M \in \mathbb{N}$
, there exist constants
$\{a_j\}_{j=1}^M$
and a positive constant
$C_{W}$
that is independent of
$M$
such that

Remark 2.7. Although the convergence rate becomes worse, we can relax the condition (2.24) as

The evaluation of uniform convergence up to derivatives of the Lagrange interpolating polynomials with Chebyshev nodes for
$L$
and
$n$
, such as using the Lebesgue constant, is an open problem. Examples of the potential
$W$
are
$\cosh j(L-|x|)$
and
$-(\!\cosh ^2(L-|x|)-1)\cos (\!\cosh (L-x))$
. In the former case,
$f(x) = T_j(x)$
, and in the latter case,
$f(x) = -(x^2-1) \cos x$
for
$x \in [1,\cosh L]$
.
There are degrees of freedom in the choices of
$d_j$
and
$k_j$
, that is, the second equation of the approximation equation (
$\mbox{KS}^{M,\varepsilon }$
). It is not known whether the above choices are optimal in terms of convergence of approximation errors and dependence on the number of approximation equations, and is an open question. On the other hand, as for the second equation of (
$\mbox{KS}^{M,\varepsilon }$
), the superposition of fundamental solutions as in Corollary 2.6 can approximate even functions in
$C^1(\Omega )$
, and in this sense, (
$\mbox{KS}^{M,\varepsilon }$
) is suitable for approximating nonlocal Fokker–Planck equation.
Because we estimate the error of the solutions of the two nonlocal Fokker–Planck equations with
$W=\sum _{j=1}^M a_j k_j$
and any given potential
$W(x)$
, we prepare the following lemma.
Lemma 2.8.
Suppose that
$w_1, w_2 \in W_{\mathrm{per}}^{1,1}(\Omega )$
and let
$\rho _j, \ (j=1,2)$
denote the solution to

with the periodic boundary condition

respectively. Then for any
$T\gt 0$
, there exists a positive constant
${\tilde C_T} = {\tilde C_T} (\rho _0, L, T, \|w_{1,x} \|_{L^1(\Omega )}, \|w_{2,x} \|_{L^1(\Omega )})$
such that

This lemma shows that the difference between the solutions to the two Fokker–Planck equations is bounded by the difference between the two potentials. The proof is presented in Subsection 4.3
Referring to
$\{ \alpha _j^{M-1}\}$
in Theorem 5.3, we put

The approximation of
$W$
in Theorem 5.3 requires a constant term with
$j=0$
in
$\sum _{j=0}^n\alpha _j^n \cosh j(L-|x|)$
. From this requirement, it is necessary to define
$d_1$
as a sufficiently large value. We estimate the differences between the two solutions to (P) with an arbitrary even
$W$
in
$C^\infty (\Omega )$
and
$W = \sum _{j=1}^M a_j k_j$
by using (2.25) in Corollary 2.6 and (2.26) in Lemma 2.8. Moreover, because we can estimate the solutions to (
$\mbox{KS}^{M,\varepsilon }$
) and (P) with
$W = \sum _{j=1}^M a_j k_j$
from Theorem 2.4, we obtain the following main result.
Theorem 2.9.
For any even 2L-periodic function
$W$
in
$C^\infty (\Omega )$
with (2.23) and (2.24) for any
$n \in \mathbb{N}$
, any time
$T\gt 0$
, any
$M \in \mathbb{N}$
and any
$\varepsilon \gt 0$
, there exists a Keller–Segel system (
$\mbox{KS}^{M,\varepsilon }$
) with
$M+1$
component, a positive constant
$C_T^{(1)}$
that is independent of
$M$
and
$\varepsilon$
, and a positive constant
$C_T^{(2)}=C_T^{(2)}(M)$
that is independent of
$\varepsilon$
such that

where
$\rho$
is the solution to (P) equipped with
$\rho _0 \in C^2(\Omega )$
and
$\rho ^{M,\varepsilon }$
is the first component of the solution to (
$\mbox{KS}^{M,\varepsilon }$
) equipped with (2.17) and (2.18).
Remark 2.10.
The Sobolev embedding theorem yields that the convergence (2.28) holds in
$C([0,T];C(\Omega ))$
. We note that the limit of
$\varepsilon \to 0$
can be taken in (2.28). It indicates that
$({\rm{KS}}^{M,0})$
can also approximate (P). Moreover, the limit of
$M \to \infty$
can be taken in
$({\rm{KS}}^{M,0})$
.
This theorem shows that the solution to the nonlocal Fokker–Planck equation with any potential can be approximated by that to the multiple components of the Keller–Segel system with specified parameters. This convergence result shows a relationship between any advective nonlocal interactions and local dynamics. Furthermore, the convergence rate is specified with respect to
$M$
and
$\varepsilon$
. Since the right-hand side of (2.28) includes the power and factorial terms, this indicates that the error in
$({\rm{KS}}^{M,0})$
can converge to
$0$
rapidly.
3. Mild solution
To show Theorem 2.2, we present some lemmas. Throughout this section, we set
$d_j\gt 0$
, and
$a_j$
are constants for
$j=1,\ldots ,M$
.
3.1. Fundamental solution and boundedness
First, we provide the following lemma:
Lemma 3.1.
Functions
$G$
of (2.9) and
$G_j^{\varepsilon }$
of (2.12) are the fundamental solution for
$u_t = u_{xx}$
and

with periodic boundary condition of (2.3) type, respectively.
The proof of this lemma is obtained by substituting and calculating.
Next, we have the following lemma.
Lemma 3.2.
Assume
$\phi \in E_{\tau }$
. Then (2.13) satisfies

The proof of this lemma is provided by substituting (2.13) into the equation and calculating. Next, we estimate the boundedness of
$\Psi _j[\phi ]$
in the following lemma.
Lemma 3.3.
Assume that
$ \phi \in E_{\tau }$
and (2.16), and let
$C_3$
be a positive constant given by
$C_3 \,:\!=\,\sqrt { L} /(d_j \sqrt {6})$
. Then, (2.13) satisfies the following estimates


Proof of Lemma 3.3. From the triangular inequality, we have

By the maximum principle for the heat equation (3.29), we compute the first term.

The last equality is based on the Fourier series expansion. Next, we denote the Fourier coefficient of
$\phi$
by

Before estimating the second term of (3.32), we can compute that

Then, we see that


Thus, we obtain the first estimation. Replacing the functions
$(v_{j})_{0,x}$
and
$ \phi$
with
$(v_{j})_{0,xx}$
and
$ \phi _{x}$
in (3.32), respectively, we have also the second assertion of the boundedness by the same calculation as above.
Because we will use the following boundedness several times, we give the following lemma:
Lemma 3.4.
Assume
$f=f(x) \in L^2(\Omega )$
and
$ g = g(x,t)\in C( [0,T];L^2(\Omega ))$
for any
$T\gt 0$
, and let
$C_4$
be a positive constant given by
$C_4\,:\!=\,\pi /(d_jL)$
. Then, we obtain that, for all
$t\in (0,T]$
,





We note that the right-hand sides of (3.36) and (3.37) do not depend on
$\varepsilon$
. The proof is based on the estimation of the Fourier coefficients of
$f$
and
$g$
. As this is a straightforward calculation, we present the proof in A.
3.2. Contraction map
We show that the map
$\Phi$
becomes a contraction map from

to
$B_R$
taking the sufficiently small time
$\tau \gt 0$
.
Lemma 3.5.
Assume that
$\phi \in B_R$
and (2.16). Then, (2.13) satisfies the following estimate

where

Since this proof is straightforward, we put it in B. Next, we have the following lemma.
Lemma 3.6.
Assume that
$\phi \in B_R$
and (2.15). Then, (2.14) satisfies

Proof of Lemma 3.6. From the Minkowski inequality, we see that

Then, by using (3.33) in Lemma 3.4, we can compute the first term of (3.38) as follows:

Next, we estimate the second term of (3.38). Utilising (3.34) in Lemmas 3.4 and 3.5, we compute that

Similarly, (3.35) in Lemma 3.4 yields that

Consequently, we observe that
$\Phi : E_{\tau } \to E_{\tau }$
, and that there exists
$\tau _1 \gt 0$
independent of
$\varepsilon$
such that for
$\tau \lt \tau _1$
$\Phi$
is a map from
$B_R$
to
$B_R$
.
Lemma 3.7.
Assume that
$\phi , \psi \in E_{\tau }$
. Then, (2.13) satisfies that


where
$C_4$
denotes the constant defined in Lemma 3.4
.
Since the proof is simply based on using Lemma 3.4, it is put in B. Finally, we show that map
$\Phi$
becomes a contraction map by taking a sufficiently small
$\tau _2\gt 0$
and setting
$\tau = \tau _2$
.
Lemma 3.8.
Assume that
$\phi , \psi \in E_{\tau }$
. Then there exists a positive constant
$C_5$
independent of
$\varepsilon$
such that

Proof of Lemma 3.8. Since

we estimate each term on right-hand side. We can compute that

Using (3.35) in Lemma 3.4, we estimate that

where the boundary integration from the integration by parts in the first equality becomes
$0$
because of the periodicity, we used
$\phi \in C(\Omega \times [0,\tau ])$
from the Sobolev embedding theorem, the Minkowski inequality, and (3.39) in Lemma 3.7.
Similarly, to this estimation, (3.34) in Lemma 3.4 yields that

where we utilised the Minkowski inequality, the boundedness (3.30) and (3.31) in Lemma 3.3, and

We note that
$C_6$
does not depend on
$\varepsilon$
.
Next, we estimate the second term of (3.41) on right-hand side. First, we write as

Similarly, to the previous estimations, using (3.35) in Lemma 3.4, we obtain that

where we used the boundedness (3.30) in Lemma 3.3, the boundedness that
$\phi \in C([0,\tau ]; C(\Omega ))$
, (3.39) and (3.40) in Lemma 3.7, and we put

It should also be noted that
$C_7$
does not depend on
$\varepsilon$
.
Finally, using (3.35) in Lemma 3.4, we obtain that

where
$C_6$
is as defined in (3.42). Putting

we complete the proof.
Consequently, taking a sufficiently small value
$\tau =\tau _2\gt 0$
which is independent of
$\varepsilon$
, we obtain

Thus, the map
$\Phi \,:\, B_R \to B_R$
is a contraction map.
Proof of Theorem
2.2. By setting
$\tau _0 \,:\!=\, \min \{ \tau _1, \tau _2\}$
, we see that the map
$\Phi : B_R \to B_R$
is a contraction map. From the Banach fixed-point theorem, the equation
$\rho ^{M,\varepsilon } = \Phi [\rho ^{M,\varepsilon }]$
has a unique solution in
$C([0,\tau _0]; H^1(\Omega ))$
. Since
$ \Phi [\rho ^{M,\varepsilon }](-L,t)=\Phi [\rho ^{M,\varepsilon }](L,t)$
and
$\Phi _x[\rho ^{M,\varepsilon }](-L,t)=\Phi _x[\rho ^{M,\varepsilon }](L,t)$
for
$t\gt 0$
, this mild solution satisfies the periodic boundary condition.
Proof of Corollary 2.3. Repeating to use Theorem 2.2, we can connect the mild solutions for
$t \in [0,T]$
at any time
$T\gt 0$
. Using the term-by-term of the weak derivative for the integral equations with respect to
$x$
and
$t$
, we observe that the solutions
$\rho ^{M,\varepsilon }$
and
$v_j^{M,\varepsilon }$
satisfy (
$\mbox{KS}^{M,\varepsilon }$
) in
$L^2(0,T;L^2(\Omega ))$
and
$L^2(0,T;C(\Omega ))$
, respectively.
Differentiating the right-hand sides of
$\rho ^{M,\varepsilon } = \Phi [\rho ^{M,\varepsilon } ]$
and
$v_j^{M,\varepsilon } = \Psi _j[ \rho ^{M,\varepsilon } ]$
with respect to
$x$
by applying the term-by-term differentiation, respectively, we can show that
$\rho ^{M,\varepsilon } \in L^2( 0, T; H^2(\Omega ) )$
and
$v_j^{M,\varepsilon } \in L^2( 0, T; H^3(\Omega ) )$
. Similarly, differentiating the right-hand side of
$\rho ^{M,\varepsilon } = \Phi [\rho ^{M,\varepsilon }]$
with respect to
$t$
by applying the term-by-term differentiation, we see that
$\rho ^{M,\varepsilon } \in H^1(0,T; L^2(\Omega ))$
.
4. Singular limit analysis
To demonstrate Theorem 2.4, we prepare the lemmas in the following subsections. Throughout this section, we set
$d_j\gt 0$
, and
$a_j$
are constants for
$j=1,\ldots ,M$
.
4.1. Fundamental solution
First, we have the following lemma.
Lemma 4.1.
$k_j$
defined in (2.4) is the fundamental solution to

where
$\delta$
denotes the Dirac delta function. Moreover,

holds in the weak sense.
Proof. For the proof of the second assertion, we can directly obtain the weak derivatives by multiplying
$k_j$
by the test function
$\varphi _x \in C_0^\infty (\Omega )$
, and the integration by parts. For any
$\varphi \in C_0^\infty (\Omega )$
, we can compute that

This implies that

Lemma 4.2.
Let
$C_8$
be a positive constant given by

Then,
$k_j$
satisfies that
$\| k_j \|_{ L^1 (\Omega )} = 1,\ \| (k_j)_x \|_{ L^1 (\Omega )} = C_8$
, and

As the proof is elementary, we omit it.
4.2. Boundedness of auxiliary factors
Next, we estimate the boundedness of the solutions
$(\rho ^{M,\varepsilon },\{v_j^{M,\varepsilon }\}_{j=1}^M)$
to (
$\mbox{KS}^{M,\varepsilon }$
) in this subsection. First, we obtain the following lemma.
Lemma 4.3.
Let
$\rho ^{M,\varepsilon }$
be the solution to the first equation of (
$\mbox{KS}^{M,\varepsilon }$
) with
$\rho _0 \in C^2(\Omega )$
. Then,
$\rho ^{M,\varepsilon } \in C^1([0,T];L^2(\Omega ))$
and there exist positive constants
$C_9$
and
$C_{10}$
that depend on
$(\rho _0)_{xx}, M, \{a_j,d_j\}_{j=1}^M, L, R$
and
$T$
but are independent of
$\varepsilon$
such that

Although we should explicitly write the dependence of
$t$
on the norm of the spatial direction for functions depending on the position
$x$
and time
$t$
, for example,
$\|\cdot \|_{L^2(\Omega )}(t)$
, we abbreviate the symbol of
$(t)$
for the simple descriptions from here.
Proof of Lemma
4.3. First, we denote the Fourier coefficient of
$ ( \rho ^{M,\varepsilon } (\sum _{j=1}^M a_j v_j^{M,\varepsilon })_x )_x$
by

Since
$\rho ^{M,\varepsilon }$
is the mild solution, we have

We will estimate each term. Since
$G_t*\rho _0 = G_{xx} * \rho _0 = G* (\rho _0)_{xx}$
, the maximum principle yeilds that

From the Fourier series expansion, we obtain that

where the last term is bounded by Lemma 3.5. Finally, we compute that

Putting

we obtain the first and the second assertions. The Young inequality and Lemma 4.2 implies the third and last assertions, where we put

Now, we estimate the difference between the solutions in the following auxiliary equations


for
$j= 1, \ldots , M$
, where the third and fourth conditions of (2.10) are imposed, respectively, and
$\rho ^{M,\varepsilon }$
denotes the solution to (
$\mbox{KS}^{M,\varepsilon }$
). We note that the solution to (4.45) is given by
$v_j = k_j * \rho ^{M,\varepsilon }$
. We set the difference as

Lemma 4.4.
Let
$(\rho ^{M,\varepsilon },v_j^{M,\varepsilon })$
and
$v_j$
be the solutions to (
$\mbox{KS}^{M,\varepsilon }$
) with
$\rho _0\in C^2(\Omega )$
and (2.18) and (4.45), respectively, and let
$C_9$
and
$C_{10}$
be positive constants in Lemma 4.3
. Then, for any time
$T\gt 0$
the following estimates hold:

When excluding the terms multiplied by
$d_j$
on the left-hand sides, the above inequality holds without
$T/2$
on the right-hand sides.
Proof of Lemma 4.4. Taking the difference between equations of (4.44) and (4.45), we see that

Multiplying this equation by
$V_j^\varepsilon$
, integrating over
$\Omega$
and using Lemma 4.3, we have

Applying to the classical Gronwall lemma to

we have that

from the initial conditions given in (2.18). Furthermore, integrating (4.47) over
$(0,T)$
, we see that

Applying
$\sup _{t \in [0,T]}$
to (4.48) and adding it to (4.49), we obtain the first assertion.
Similarly, multiplying (4.46) by
$-\big(V_{j}^\varepsilon \big)_{xx}$
and integrating over
$\Omega$
, we have

where we used Lemma 4.3. From the Gronwall inequality, we have

by the initial condition (2.18). Integrating (4.50) over
$(0,T)$
, we see that

Applying
$\sup _{t \in [0,T]}$
to (4.51) and adding it to (4.52), we have the second assertion.
Because of
$v_j^{M,\varepsilon } \in L^2(0,T,H^3(\Omega ))$
, the similar calculation can be applied. Differentiating (4.46) with respect to
$x$
and multiplying it by
$-(V_j^\varepsilon )_{xxx}$
, we see that

Thus, the Gronwall lemma yields that

from the initial condition given in (2.18). Integrating (4.53) over
$(0,T)$
with respect to
$t$
and (4.54) implies the final assertion of this lemma.
4.3. Order estimation
Under the above preparation, we estimate the difference of solutions. Set the difference between the solutions to the first component of (
$\mbox{KS}^{M,\varepsilon }$
) and (P) with
$W= \sum _{j=1}^M a_j k_j$
as

We will show the following convergence.
Lemma 4.5.
Suppose that
$M$
is an arbitrarily fixed natural number. Let
$\rho$
be the solution to (P) equipped with
$W= \sum _{j=1}^M a_j k_j$
and the initial value
$\rho _0 \in C^2(\Omega )$
, and let
$\rho ^{M,\varepsilon }$
be the solution to the first component of (
$\mbox{KS}^{M,\varepsilon }$
) with (2.17) and (2.18). Then, for any
$T\gt 0$
, there exists a positive constant
$C_{11} = C_{11}((\rho _0)_{xx}, M, \{a_j, d_j\}_{j=1}^M, L, R, T )$
such that for any
$\varepsilon \gt 0$

We note that
$C_{11}$
is independent of
$\varepsilon$
.
Proof of Lemma
4.5. Taking the difference between the first equation of (
$\mbox{KS}^{M,\varepsilon }$
) and the equation of (P), we have

Subsequently, multiplying (4.56) by
$U^\varepsilon$
and integrating over
$\Omega$
, we obtain

where each term of the integral is set as

respectively. First, we compute
$I_1$
. Using the Cauchy–Schwartz inequality, Sobolev embedding theorem and Lemma 4.4, we have

where

and
$C_{\mbox s}$
is the constant from the Sobolev embedding theorem. Next, we compute that

where we used the estimate in Lemma 4.2 and

Finally, the Sobolev embedding theorem, the Young inequality and the boundedness in Lemma 4.2 yield that

where the constant is defined as

Summarising these estimations, we have

where we put
$C_{15}\,:\!=\, 2(C_{13}^2 + 2C_{14}^2 )$
. Applying the classical Gronwall inequality with the initial condition (2.17) to

we have

Furthermore, integrating (4.57) over
$(0,T)$
, we also obtain that

Defining
$C_{11} \,:\!=\, C_{12} (e^{ C_{15} T } - 1)/ C_{15} + 4 C_{12} e^{C_{15} T} T$
and adding (4.58) and (4.59) imply the assertion of this lemma.
Similarly to this lemma, we can obtain the following estimate.
Lemma 4.6.
Suppose the same assumptions as Lemma 4.5
. Then, for any
$T\gt 0$
, there exists a positive constant
$C_{16} = C_{16}((\rho _0)_{xx}, M, \{a_j, d_j\}_{j=1}^M, L, R, T)$
such that for any
$\varepsilon \gt 0$

We note that
$C_{16}$
is independent of
$\varepsilon$
.
Proof of Lemma
4.6. Similarly to Lemma 4.5, multiplying (4.56) by
$-U^\varepsilon _{xx}$
and integrating over
$\Omega$
, we have

where we defined each term of energy term by the integral as

First, we estimate
$\mathscr{I}_1$
. From Lemma 4.4 and the Sobolev embedding theorem, there exists a positive constant
$C_{\mbox s}$
such that
$\left \| V_{j,x} \right \|_{ C (\Omega )} \le C_{\mbox s} \left \| V_{j,x} \right \|_{ H^1 (\Omega )}$
. Then, we see that

where

Next, we compute
$ \mathscr{I}_2$
as

where

From
$\rho ^{M,\varepsilon } \in C(\Omega )$
and Lemma 4.2, we see that

where we put

Similarly to that of
$\mathscr{I}_3$
, we obtain that

where we used the estimate (4.55) in Lemma 4.5 and we put

Hereafter, we often use the estimate (4.55) in Lemma 4.5. We compute
$\mathscr{I}_5$
as

where we used (4.43) in Lemma 4.2 and (4.55) in Lemma 4.5, and we put

Similarly, we see that

where

Combining these estimation and setting a positive constant as

we have

Applying to Gronwall inequality to (4.60) without
$-\left \| U^\varepsilon _{xx} \right \|_{ L^2 (\Omega )}^2/4$
, we obtain that

by the initial condition in (2.17).
Finally, integrating (4.60) over
$(0,T)$
, we obtain that

Adding (4.61) and (4.62) yields the assertion of this lemma, where we put

Then, we obtain the following proof.
Proof of Theorem 2.4. Putting

where we used Lemma 4.5, Lemma 4.6 for
$C_1$
and Lemma 4.4 for
$C_2$
. We obtain the convergence of the assertion in this theorem.
Next, we prove Lemma 2.8.
Proof of Lemma 2.8. Set the differences between the solutions and between the kernels as

respectively. The method for this proof is similar to that of Lemmas 4.5 and 4.6. Taking the difference between the equations
$(\mbox{P}_1)$
and
$(\mbox{P}_2)$
, we have

Multiplying it by
$U$
and integrating it over
$\Omega$
, we have

Since

we can compute that

where we put

and
$C_{\mbox s}$
is the coefficient from the Sobolev embedding theorem, and
$C_0^{(1)}$
and
$C_0^{(2)}$
correspond to the constant
$C_0$
in Proposition 2.1 for
$\rho _1$
and
$\rho _2$
, respectively. Applying the Gronwall inequality to this, we have

Integrating (4.64) over
$[0,T]$
and adding it and (4.65), we have

Next, multiplying (4.63) by
$-U_{xx}$
and integrating it over
$\Omega$
, we consider the following equation of energy

Since

we can obtain that

with suitable positive constants from
$C_{26}$
to
$C_{29}$
, where we put

Therefore, the Gronwall inequality yields that

Finally integrating (4.64) and (4.66) over
$(0,T)$
and adding them and (4.65) and (4.67), we obtain the assertion of Lemma 2.8 by putting

5. Coefficients of linear sum
We now explain the method used for determining the coefficient
$\{ a_j \}_{j=1}^M$
of the linear sum of the fundamental solution for a given even potential function
$W$
. Furthermore, we will perform the numerical simulations of the approximation of
$W$
by sum of
$\cosh j(L-|x|)$
, and numerical simulations of (P) and (
$\mbox{KS}^{M,\varepsilon }$
) with this series expansion. Since
$W$
is even, we only consider
$[0,L]$
. Throughout this section, we set
$d_1$
is sufficiently large and
$d_j=1/(j-1)^2$
for
$j=2,\ldots ,M$
.
First, we provide the following lemma with respect to the
$n$
degree Chebyshev polynomial
$T_n(x) \,:\!=\, \cos (n\arccos x)$
. We set coefficients as

for
$n \in \mathbb{N}$
, where
$[\!\cdot\!]$
denotes the Gauss symbol. By this constant, the Chebyshev polynomial of the first kind of
$n$
degree
$T_n$
can be expressed as
$T_n(x) = \sum _{k=0}^{[n/2]}C^n_k x^{n-2k}$
for
$x\in [-1,1]$
. The properties of the Chebyshev polynomial can be referred from [Reference Mason and Handscomb16] by Mason and Handscomb. Utilising this equation, we have the following Lemma regarding the change of the variable for
$T_n$
.
Lemma 5.1. Setting

for
$n \in \mathbb{N}$
, we define the coefficient as

Then,

holds.
Proof of Lemma 5.1. We compute that

where we used the binomial expansion in the third equality.
Next, we explicitly provide the coefficient of the linear sum of the
$n$
degree Lagrange interpolation polynomial with the Chebyshev nodes for the arbitrary function
$F=F(x)$
for
$x\in [a,b]$
. We will replace the arbitrary function
$F$
with the function
$f$
defined in (2.22) to prove Theorem 5.3. The root of the
$n$
degree Chebyshev polynomial, called Chebyshev nodes, in an arbitrary interval
$[a,b]$
is given by

We have that

Moreover, setting the coefficient as

for
$n \in \mathbb{N}$
, we see that the
$n$
degree Chebyshev polynomial for the function
$F$
is given by

Note that
$L_n( r^{n+1}_j ) = F ( r^{n+1}_j )$
. Then, we obtain the following proposition.
Proposition 5.2. Set

for
$n \in \mathbb{N}$
. Then, the
$n$
degree Lagrange interpolation polynomial
$L_n$
for an arbitrary function
$F$
on
$[a,b]$
can be described as

Proof of Proposition 5.2. Using Lemma 5.1 and (5.68), we can compute that

As
$x - r^{n+1}_j$
for
$j = 0,\ldots ,n$
are the factors of
$T_{n+1}( (2x - (b+a) )/(b-a) )$
, respectively,
$\sum _{k=0}^{n+1} \xi ^{n+1}_k x^k/(x - r^{n+1}_j )$
must be divisible from the factor theorem. Thus, we obtain that

Before the proof of Theorem 5.3, we introduce the following constant for
$n \in \mathbb{N}$
and
$j \equiv n\pmod 2$
:

Using this constant, we can describe the formula as
$x^n = \sum _{j=0, j \equiv n\pmod 2}^n \delta ^n_j T_j(x),$
[Reference Mason and Handscomb16]. Differentiating it, we then obtain that
$x^n = \sum _{j=0, j \equiv n\pmod 2}^n (j+1)\delta ^{n+1}_{j+1} U_j(x)/(n+1)$
, where
$U_n(x)\,:\!=\, (\sin (n+1)\theta )/\sin \theta , \ (x=\cos \theta )$
is the Chebyshev polynomial of the second kind. In addition, from the induction, we note that
$T_n(\!\cosh (L-|x|))= \cosh n(L-|x|)$
for
$n \in \mathbb{N}$
holds. Using
$f$
in (2.22) instead of
$F$
in (5.69), we reconsider the coefficients
$\zeta ^n_j$
in (5.69) and
$b^n_l$
in (5.70) in Theorem 5.3.
Theorem 5.3.
Assume that
$W \in C^m([0,L])$
for given
$2 \le m \in \mathbb{N}$
and (
2.23
) for
$n \le m-1$
. Let
$\zeta ^n_j$
and
$b^n_l$
be (
5.69
) with
$f$
and (
5.70
) with a natural number
$n \le m-1$
for
$a=1$
and
$b=\cosh L$
, respectively. Set the coefficient
$\alpha ^{n}_j$
as

Then, there exists a constant
$C_L$
that is independent of
$n$
such that for any
$n \le m-1$

holds.
We note that this theorem can be applicable to the smooth even potential
$W$
with (2.23) on
$\Omega$
since
$\cosh j(L-|x|)$
is even.
Proof of Theorem 5.3. Recall that
$f(y) = W(L- \cosh ^{-1} y)$
. From the property of the Lagrange polynomial for
$f(y)$
defined in (2.22) for any
$y \in [1, \cosh L]$
, there exists a constant
$\min _i\{ r^{n+1}_i \} \lt c \lt \max _i\{ r^{n+1}_i\}$
such that

with
$b=\cosh L$
and
$a=1$
and thus,

Here, we used
$\| T_n \|_{C[0,1]}=1$
because of the definition. Using Proposition 5.2 and changing the variable to
$y = \cosh (L-x)$
for
$x \in [0,L]$
, we can compute the right-hand side of the above inequality as

Next, differentiating (5.71), we see that

because of
$T'_{n+1}(x) = (n+1) U_n(x),$
[Reference Mason and Handscomb16]. Setting the function as
$\mathscr{G} = \sqrt {y^2-1}$
for
$y \in [1,\cosh L]$
, we see that
$\mathscr{G}(\!\cosh (L-x)) = \sinh (L-x)$
for
$x \in [0,L]$
. Then, we can compute that

where we used
$x^n = \sum _{j=0, j \equiv n\pmod 2}^n (j+1)\delta ^{n+1}_{j+1} U_j(x)/(n+1)$
in the second equality. Thus, we obtain that

Here, we used
$\| U_n \|_{C([0,1])} = U_n(0)=n+1$
because the value of
$(\sin (n+1)\theta )/\sin \theta$
on extreme points in
$\theta \in [0,\pi /2]$
is decreasing. Putting
$C_L\,:\!=\,\max \{ (\!\cosh L -1)/4, \sinh L \}$
implies the assertion of this theorem.
As any continuous function can be approximated by the sum of
$(\!\cosh (L-x))^j$
by Theorem 5 in [Reference Ninomiya, Tanaka and Yamamoto19], if the coefficients of
$(\!\cosh (L-x))^j$
, that is
$b^n_j$
, are obtained, it was conceived that the error between
$W$
and
$\sum _{j=1}^n a_j k_j$
could be estimated. Indeed, the smooth function can be approximated by the Lagrange interpolation polynomial. Thus, it was considered that it should be shown that
$L_n$
can be expressed by the form of
$\sum _{j=0}^n b_j y^j$
. This is generally unknown. However, when the Chebyshev nodes are given to the Lagrange interpolation polynomial, computing
$b^n_j$
resulted in the proof since
$T_n(x) = \sum _{k=0}^{[n/2]}C^n_k x^{n-2k}$
holds and the numerator of the Lagrange polynomial must be divisible by
$x-r^{n+1}_j$
. In this sense, Proposition 5.2 is essential.
Next, we prove Corollary 2.6
Proof of Corollary 2.6. Using Theorem 5.3 and
$a_j$
defined in (2.27), we can obtain the same estimate as the proof of Theorem 5.3.
Then, we explain the proof of Theorem 2.9.
Proof of Theorem 2.9. According to Theorem 5.3, for any even function
$W$
in
$C^\infty ([0,L])$
with (2.23) and (2.24) for any
$n\in \mathbb{N}$
, and arbitrary
$ M\in \mathbb{N}$
, there exist constants
$\{ \alpha _j^{M-1}\}_{j=0}^{M-1}$
such that

Here, we set
$n=M-1$
in Theorem 5.3. Putting the parameters as (2.21) and (2.27), then we have

Let
$\bar { \rho }(x,t)$
be the solution to (P) with the integral kernel
$\sum _{j=1}^M a_j k_j$
. Since (5.72) holds and
$C(M)$
converges to zero as
$M \to \infty$
, there exists a positive constant
$\delta$
that is independent of
$M$
such that
$\|\sum _{j=1}^M a_j k_j\|_{W^{1,1} (\Omega )} \le \delta + \|W \|_{W^{1,1}(\Omega )}$
. Thus, we see that
$\tilde C_T$
in Lemma 2.8 does not depend on
$M$
. Then, for this fixed
$M$
, Theorem 2.4 and Lemma 2.8 yield that for any
$\varepsilon \gt 0$

Thus, we put
$C_T^{(1)}\,:\!=\, 2L C_L \sqrt { {\tilde C_T} }$
and
$C_T^{(2)}(M)\,:\!=\, C_1$
, respectively.
We performed a numerical simulation of the approximation for the potential
$W$
by the linear combination of
$\cosh j(L-|x|)$
. The results are shown in Figure 3. The linear combination of
$\cosh j(L-|x|)$
covers a potential
$W$
. The longer the length of the interval
$L$
becomes, the worse the rate of convergence becomes from the numerical simulations. However, as the rate of convergence can be exponential as given by Theorem 5.3, the method for determining the coefficient of
$a_j$
is compatible with and useful for numerical simulations.

Figure 3. Results of a numerical simulation of the approximation for
$W$
by the linear combination of
$\cosh j(L-|x|)$
. We set
$W(x)=e^{-5x^2}(\cos (3\pi ) x- 1/2 \cos (2\pi x))$
, and
$L=2$
. (a) Profiles of
$W$
and the linear sum of
$\cosh j(L-|x|)$
. (b) Profiles of
$f$
and the Lagrange interpolation polynomial on
$[1, \cosh L]$
. (c) Distribution of
$\{ \alpha ^9_j \}_{j=0}^9$
.
Figures 4 and 8 shows the numerical results of (P) with the potential
$W(x)=e^{-5x^2}$
and (
$\mbox{KS}^{M,\varepsilon }$
) with parameters
$\varepsilon =0.001$
and
$\{\alpha _j^6\}^6_{j=0}$
specified by Theorem 5.3. We can observe that the solution
$\rho$
in (P) is approximated by that
$\rho ^{7,\varepsilon }$
in (
$\mbox{KS}^{M,\varepsilon }$
), even though there are seven auxiliary factors
$v_j^{7,\varepsilon }$
. In Figure 4, (c) shows the profiles of both
$W(x)=e^{-5x^2}$
and
$\sum _{j=0}^6 \alpha _j^6 \cosh (j(L-|x|))$
. Since
$\sum _{j=0}^6 \alpha _j^6 \cosh (j(L-|x|))$
has good accuracy for the approximation of
$e^{-5x^2}$
, both curves are seen to overlap.

Figure 4. Results of numerical simulations for (P) with a potential
$W(x) = e^{-5x^2}$
and
$\mu$
defined in (
$\mbox{P}_\mu$
), and (
$\mbox{KS}^{M,\varepsilon }$
) with
$M=7$
. The parameters are given by
$L=1$
,
$\varepsilon =0.001$
,
$d_1=1000000$
,
$\mu =5$
and
$d_j$
and
$a_j$
are provided by (2.21) and (2.27), respectively. (a) Profiles of the numerical result of (P) at
$t=200.0$
. The horizontal and vertical axes correspond to the position
$x$
and
$\rho$
, respectively. The red curve is the numerical result of
$\rho$
. (b) Profiles of the numerical result of (
$\mbox{KS}^{M,\varepsilon }$
) at
$t=200.0$
. We impose the same initial data for
$\rho ^{7,\varepsilon }$
as that of
$\rho$
and
$(v_j)_0= k_j*\rho _0, \ (j=1,\ldots ,M)$
. The axes are set same as that of (a). The red and the other colour curves correspond to
$(\rho ^{7,\varepsilon },\{v_j^{7,\varepsilon }\}_{j=1}^7)$
, respectively. (c) Profiles of
$W$
and
$\sum _{j=0}^6 \alpha _j^6 \cosh (j(L-|x|))$
. The orange dashed and blue curves corresponding to
$W$
and
$\sum _{j=0}^6 \alpha _j^6 \cosh (j(L-|x|))$
, respectively, are drawn in a same plane. (d) The distribution of
$\{\alpha _j^6\}_{j=0}^6$
.
6. Linear stability analysis
In this section, we perform a linear stability analysis around the equilibrium point for (P) and (
$\mbox{KS}^{M,\varepsilon }$
) with three components to specify the role of advective nonlocal interactions in pattern formations. We demonstrate that the eigenvalue of the linearised operator of (
$\mbox{KS}^{M,\varepsilon }$
) converges to that of (P) as
$\varepsilon \to 0$
when the integral kernel is given by
$k_j$
of (2.4). We analyse the following equation with the parameter
$\mu \gt 0$
:

We explain the instability of the solution near the equilibrium point. Let
$\rho ^*\gt 0$
and
$ \xi =\xi (x,t)$
be an arbitrary constant and a small perturbation, respectively.
$\rho = \rho ^*$
becomes a constant stationary solution of (P). Putting
$\rho (x,t) = \rho ^* + \xi (x,t)$
and substituting it for (
$\mbox{P}_\mu$
), we have

Focusing on the linear part of above, we denote the linear operator
$\mathscr L$
by

Because this linearised operator has
$\mu \rho ^*$
, the effects of the strength of aggregation and the mass volume on the pattern formation around the constant stationary solution are equivalent. Therefore, we replace
$ \mu \rho ^*$
with
$\mu$
. Defining the Fourier coefficient of
$W$
as

we have the following lemma with respect to the eigenvalues and eigenfunctions:
Lemma 6.1. Setting the eigenvalues

then we have

Proof. The proof follows from a direct calculation

Using this lemma, we find the solution to
$\xi _t = \mathscr L [\xi ]$
around
$\rho ^*$
in the form of
$ \sum _{n\in \mathbb{Z}} \hat {\xi }_n e^{\lambda _n t} e^{i \sigma _n x }$
, where
$\{ \hat {\xi }_n \}$
is the Fourier coefficient.

Figure 5. Results of a numerical simulation for (
$\mbox{P}_\mu$
) with (2.6). The parameters are
$\mu =5.0$
,
$d_1=0.1$
and
$d_2=3.0$
, and the initial data are given by
$1.0$
with small perturbations. The horizontal and vertical axes correspond to the position
$x$
and value of solution
$\rho$
, respectively. The red curve corresponds to the solution
$\rho$
. The left, middle left, middle right and right pictures exhibit the profiles of solutions of (
$\mbox{P}_\mu$
) with (2.6) in the interval
$[0, 10]$
at
$t = 0, 0.5, 1.0$
and
$3.0$
, respectively.

Figure 6. Results of a numerical simulation for (
$\mbox{P}_\mu$
) with (2.7). The parameters are
$\mu =4.0$
and
$R=1.0$
, and the initial data are given by
$1.0$
with small perturbations. The horizontal and vertical axes correspond to the position
$x$
and value of solution
$\rho$
, respectively. The red curve corresponds to the solution
$\rho$
. The left, middle left, middle right and right pictures exhibit the profiles of solutions of (
$\mbox{P}_\mu$
) with (2.7) in the interval
$[0, 10]$
at
$t = 0, 0.8, 2.0$
and
$12.0$
, respectively.
Here, we recall the concept of the diffusion-driven instability in pattern formations proposed by Turing [Reference Turing21]. Diffusion-driven instability is a paradox where diffusion, typically leading to concentration homogenisation, destabilises the uniform stationary solution and induces nonuniformity due to the difference in the diffusion coefficients. By using the eigenvalue
$\lambda =\lambda (n)$
for the linear operator of the reaction-diffusion system, the diffusion-driven instability can be defined as the eigenvalue
$\lambda$
satisfies
$\lambda (0)\lt 0$
and there exists
$n\in \mathbb{Z}$
such that
$\lambda (n)\gt 0$
. For the model (
$\mbox{P}_\mu$
), the similar situation occurs. If
$W$
satisfies that
$\lim _{n \to \pm \infty } \omega _n=0$
and that there exists
$0 \neq n_1 \in \mathbb{N}$
such that
${\displaystyle } Re \omega _{n_1} \gt 0$
, the maximum eigenvalue may be attained around
$n=n_1$
. In that case, the unstable mode around a stable equilibrium point is given by
$e^{i \sigma _{n_1} x }$
.
We performed numerical simulations of (
$\mbox{P}_\mu$
) with the integral kernel (2.6) and (2.7) with the finite volume method. Figures 5 and 6 present the results. The Fourier coefficients for the integral kernels (2.6) and (2.7) are given by

respectively. Figures 1 (b) and 2 (b) show the distributions of the eigenvalues with
$\omega _{n,1}$
and
$\omega _{n,2}$
. The number of peaks of the solution at the beginning of the pattern formation corresponds to the maximum wave number in Figures 5 and 6, respectively.
For (2.6), by introducing
$v^\varepsilon _1=k_1*\rho$
and
$v^\varepsilon _2=k_2*\rho$
into
$W$
, the solution of (
$\mbox{P}_\mu$
) can be approximated by that of the 3-component Keller–Segel system from Theorem 2.4:

with
$0 \lt \varepsilon \ll 1$
. In (6.73),
$v_1^\varepsilon$
and
$v_2^\varepsilon$
represent the attractive and repulsive substances in chemotactic process, respectively. By determining the solution as a Fourier series expansion, the linearised problem is provided by the following system:

where
$(\hat {\varphi _3})_n$
is the Fourier coefficient for the perturbations of
$v_2^\varepsilon$
. The characteristic polynomial is given by

Then, we can see that

Only one eigenvalue converges to a bounded value when
$\varepsilon \to 0+0$
from the implicit function theorem. Setting the eigenvalue as
$\lambda _\varepsilon$
, we can solve that

and thus

This implies that the solution to (6.73) and (
$\mbox{P}_\mu$
) with (2.6) are sufficiently close, but also that the Fourier mode of (
$\mbox{P}_\mu$
) when the pattern forms around an equilibrium point is also extremely close to that of the 3-component attraction-repulsion Keller–Segel system (6.73). Because the constant stationary solution is destabilised by the auxiliary factors
$v_1^\varepsilon$
and
$v_2^\varepsilon$
, the mechanism of the pattern formation is almost the same as that of diffusion-driven instability. In other words, if the integral kernel is provided by (2.6), the solution to (P) is sufficiently close to that of the Keller–Segel system, which can cause the diffusion-driven instability, and thereby suggesting that the kernel
$W$
is crucial in generating the diffusion-driven instability in the nonlocal Fokker–Planck equation (P) in linear sense.
We performed a numerical simulation of (6.73) with
$\varepsilon =0.001$
. The profile of the solution
$\rho ^\varepsilon$
at each time point in Figure 7 is similar to that of
$\rho$
in Figure 5. These figures also indicate that stationary solutions to (P) may be approximated by those to (
$\mbox{KS}^{M,\varepsilon }$
). As explained above, by approximating the dynamics of nonlocal evolution equations using Keller–Segel systems, we can describe the nonlocal dynamics within the framework of local dynamics and identify both mechanisms.

Figure 7. The results of a numerical simulation for (6.73). The parameters
$\mu , d_1, d_2$
and the initial data are same as that in Figure 5, and
$\varepsilon =0.001$
and
$((v_1)_0,(v_2)_0)=( k_1*\rho _0, k_2*\rho _0 )$
. The horizontal and vertical axes correspond to the position
$x$
and value of solutions
$\rho ^\varepsilon$
,
$v_1^\varepsilon$
and
$v_2^\varepsilon$
, respectively. The red, green and blue curves correspond to the solution
$\rho ^\varepsilon$
,
$v_1^\varepsilon$
and
$v_2^\varepsilon$
, respectively. The left, middle left, middle right and right pictures exhibit the profiles of solutions of (6.73) in the interval
$[0, 10]$
at
$t = 0, 0.5, 1.0$
and
$3.0$
, respectively.

Figure 8. Comparison of the time evolutions of numerical results given in Figure 4 (a) and (b) until
$t=1.0$
.
7. Concluding remarks
We approximated the solutions of the nonlocal Fokker–Planck equation with any even advective nonlocal interactions (P) by those of multiple components of the Keller–Segel system (
$\mbox{KS}^{M,\varepsilon }$
). This indicates that the mechanism of the weight function for determining the velocity by sensing the density globally in space can be realised by combining multiple chemotactic factors. Additionally, our results show that this diffusion–aggregation process can be described as a chemotactic process. We propose a method in which the parameters
$\{d_j, a_j\}$
can be determined based on the profile of the potential
$W$
. Using the Keller–Segel type approximation, we rigorously demonstrate that the destabilisation of the solution near equilibrium points in the nonlocal Fokker–Planck equation closely resembles diffusion-driven instability. This type of analysis can be applied to other nonlocal evolution equations with advective nonlocal interactions, such as cell adhesion models.
The Keller–Segel approximation also benefits the numerical algorithm in (P). By approximating the potential
$W$
by
$\sum _{j=1}^Ma_j k_j$
using Theorem 5.3 and solving (
$\mbox{KS}^{M,\varepsilon }$
) or
$({\rm{KS}}^{M,0})$
numerically, we can remove the nonlocality from (P). By calculating these local systems instead of (P) using a simple integral scheme, a numerical simulation can be performed more rapidly. Indeed, the calculation cost as performed in Figures 4 (b) and 7 is
$O(MN)$
in the time loop iteration, where
$O$
is the Landau symbol and
$N$
is the number of the spatial mesh. Here, we used the finite volume method and LU decomposition with a tridiagonal matrix.
Theorem 2.9 indicates that local dynamics, such as the Keller–Segel system, and nonlocal dynamics, such as the nonlocal Fokker–Planck equation, can be bridged in the sense of continuous evolutionary behaviour. Thus, we can treat the problem (P) within the framework of (
$\mbox{KS}^{M,\varepsilon }$
) if (
$\mbox{KS}^{M,\varepsilon }$
) is easier. As demonstrated by the linear stability analysis of (
$\mbox{KS}^{M,\varepsilon }$
), we can characterise the solutions to (P) for local dynamics.
According to Ninomiya et al. [Reference Ninomiya, Tanaka and Yamamoto19], the existence of parameters
$\{ a_j\}$
was shown for the continuous integral kernel
$W$
although the explicit formula of the coefficients
$\{ a_j\}$
was not obtained. This suggests that the condition of Theorem 5.3 for determining
$\{ a_j\}$
for potential
$W$
may be relaxed. We aim to intensify this investigation in the future.
Financial support
The authors were partially supported by JSPS KAKENHI Grant Number 22K03444 and 24H00188. HM was partially supported by JSPS KAKENHI Grant Number 21KK0044 and by the Joint Research Center for Science and Technology of Ryukoku University. YT was partially supported by JSPS KAKENHI Grant Number 20K14364 and 24K06848 and by Special Research Expenses of Future University Hakodate. A visualisation software, GLSC3D, was used to visualise numerical solutions.
Competing interests
The authors declare that they have no conflict of interest.
Data availability statement
The source codes used to produce the numerical solutions in this manuscript are available on Zenodo at link: https://zenodo.org/records/15583545.
Appendix A. Proof of Lemma 3.4
Proof. First, we denote the Fourier coefficients of
$f$
and
$g$
by

for
$n \in \mathbb{Z}$
, respectively.
Then, using the orthogonality and the Parseval identity, we compute that

Straightforwardly, we can calculate that

where we used the Maclaurin series expansion, that is, for
$\sigma _n^2 t$
, there exists
$\theta _n \in (0,1)$
such that

Similarly, the Maclaurin series expansion (A.1) yields that

Next, we estimate the boundedness with
$G_j^\varepsilon$
. Utilising the orthogonality and the Parseval identity, we obtain that

and

Appendix B. Proofs of Lemma 3.5 and 3.7
Proof of Lemma 3.5. Using the Minkowski inequality and Lemma 3.3, we see that

Proof of Lemma 3.7. Using (3.36) in Lemma 3.4, We compute that

Similarly, (3.37) in Lemma 3.4 shows that
