Nonlocal Reaction–Diffusion Models of Heterogeneous Wealth Distribution

Abstract

Dynamics of human populations can be affected by various socio-economic factors through their influence on the natality and mortality rates, and on the migration intensity and directions. In this work we study an economic–demographic model which takes into account the dependence of the wealth production rate on the available resources. In the case of nonlocal consumption of resources, the homogeneous-in-space wealth–population distribution is replaced by a periodic-in-space distribution for which the total wealth increases. For the global consumption of resources, if the wealth redistribution is small enough, then the homogeneous distribution is replaced by a heterogeneous one with a single wealth accumulation center. Thus, economic and demographic characteristics of nonlocal and global economies can be quite different in comparison with the local economy.

1. Introduction

Distribution of wealth is a factor that defines social justice and, as such, is believed to determine the stability of society [1]. It also affects the rate of economic development; in particular, it is thought that a certain degree of inequality is necessary to stimulate the economic growth [2,3], although the issue remains controversial [4]. Patterns of wealth distribution have been attracting keen scientific attention ever since the seminal work by Pareto [5,6]. Traditional approach to this matter is based on data collection (e.g., through the tax office) and their statistical analysis, which allows one to reveal the trends such as, for instance, an increase or decrease of the inequality [7,8]. Although providing valuable information, this approach alone has a limited power to explain the reasons behind the corresponding social and financial dynamics. Mathematical modelling is increasingly regarded as a powerful research tool across life sciences and social sciences [9,10]. Over the last few decades, along with the research based on analysis of relevant financial, social and political data, there have been growing literature concerned with mathematical modelling of the wealth distribution [11,12,13,14,15,16].

Although there has been some growing interest in spatially explicit models over the last two decades [12,15], most of modelling studies available in the literature are concerned with the distribution of wealth over different social groups in a non-spatial context (e.g., see [16] and the references there), for instance in a particular area or city, or on the nationwide scale, hence largely ignoring the heterogeneity of wealth distribution across space. Meanwhile, such heterogeneity can be remarkable, even within the same country; an example is shown in Figure 1. Interestingly, while on a smaller spatial scale (e.g., down to the county level, see Figure 1, left) the distribution looks random, on a larger spatial scale of individual states level it exhibits features reminiscent of periodicity. One example of an apparently nearly periodic behavior is shown in Figure 1 (right) where small and large values of the Gini coefficient (which quantifies the inequality in the wealth distribution) alternates in the East-West direction.

Although it is not our aim in this paper to compare the modelling results with the data, Figure 1 have been an inspiration for this study, where we endeavor to look for a generic mechanism that may result in a spatially heterogeneous wealth distribution. We mention here that obviously, such heterogeneity can arise for a variety of reason. Heterogeneity of the transportation system, e.g., in the road network is one reason as it may affect the availability of jobs [18]. Arguably, heterogeneity in the distribution of natural resources may have a similar effect. Here we hypothesize that apart from the above and other similar reasons, there may exist a generic dynamical mechanism of wealth distribution that can lead to the formation of a heterogeneous spatial pattern. Mathematically, our model is described by coupled nonlinear partial differential equations of reaction–diffusion type that are extended to take into account the nonlocal nature of some of the processes. Being inspired by the success of similar approaches in describing complex population dynamics systems [19,20,21,22], in this paper we consider a spatially explicit model that links the wealth distribution to the population density.

2. Reaction–Diffusion Model of Population-Wealth Distribution

2.1. Basic Model with Constant Resources

We will describe the dynamics of human populations with a conceptual economic–demographic model suggested in [23,24]. The model quantifies the population density $p(x,t)$ and the distribution of wealth $u(x,t)$. It consists of two partial differential equations of reaction–diffusion type:

$$\frac{\partial u}{\partial t}=D_u\frac{{\partial }^{2}u}{\partial x^2}+F(u,p),$$

(1)

$$\frac{\partial p}{\partial t}=D_p\frac{{\partial }^{2}p}{\partial x^2}+G(u,p),$$

(2)

where we neglect a possible effect of cross-diffusion introduced in [23]. The first term in the right-hand side of Equation (2) describes the movement of the population in space considered to be random. The diffusion term in Equation (1) characterizes local wealth redistribution due to the economic activities such as trade and investments, and/or taxes.

The population growth in Equation (2) is described by the function

$$G(u,p)=\alpha p\left(K\right(u)-p)-\sigma \left(u\right)p$$

with the logistic growth characterized by the fertility rate $\alpha$ and the carrying capacity $K\left(u\right)$. The second term describes mortality of the population with rate $\sigma \left(u\right)$. Both the carrying capacity and the mortality rates depend on the wealth. For the mortality rate, we consider it to be a monotonously decreasing function of wealth. It takes into account the general observation that on average, the mortality rate is lower for wealthier people due to access to better health services and healthier lifestyle [25]. More specifically, we consider the dependence:

$$\sigma \left(u\right)={\sigma }_0-\frac{{\sigma }_{1}u}{1+u},$$

where ${\sigma }_0$ and ${\sigma }_1$ are positive parameters.

In order characterize the function $K\left(u\right)$, we note that it describes an equilibrium population density determined by the availability of resources [26]. When the resources decrease, the carrying capacity decreases too and may ultimately go to zero when the resources are completely depleted and hence the environment is not capable any more to support the population. Therefore, for small u, $K\left(u\right)$ is an increasing function. However, when the resource (wealth) becomes abundant, $K\left(u\right)$ can lose its monotonicity. Please note that in our formulation of the logistic growth, the linear per capita fertility rate is $\alpha K\left(u\right)$. Thus, changes in the carrying capacity reflect the changes in the fertility rate. It has been observed that that there is a certain wealth-dependent difference between low-income and high-income society groups, so that wealthier people, on average, have less children [27]. Accordingly, we consider the carrying capacity $K\left(u\right)$ to be an increasing function of wealth for small u but decreasing function for large u. Specifically, we consider the carrying capacity in the following form:

$$K\left(u\right)=\frac{a_{2}u}{u^2+c_2^2}+{\psi }_0,$$

where $a_2$ and $c_2$ are some positive parameters, and ${\psi }_0$ is a part of carrying capacity independent of wealth.

Reaction term in Equation (1) is considered in the following form:

$$F(u,p)=W(u,p)-S(u,p),$$

where W and S are the rates of the wealth production and consumption, respectively. Wealth production is described by the Cobb-Douglas dependence:

$$W=b{H}^{\nu }{Q}^{\beta }{M}^{\gamma },$$

where H is the labor, Q is the capital, and M available resources, parameter b characterizes the level of technology, $\nu ,\beta$, and $\gamma$ are positive constants [28,29]. Please note that in our model we consider “resources” in a broad sense, as any element that is essentially involved in the production of wealth and facilitates the production. Therefore, M quantifies not only the availability of natural resources, both renewable and non-renewable, but also technological resources such as qualification and skills of the labor force (and hence the availability of relevant training and education), the existence and availability of relevant software, etc. We further assume that capital Q is a function of wealth, $Q=h\left(u\right)$, and labor is a function of the population density, $H=g\left(p\right)$. Then

$$W(u,p)=h\left(u\right)g\left(p\right)M$$

(where we consider $\nu =\beta =\gamma =1$ for the sake of simplicity). We assume that $h\left(u\right)$ and $g\left(p\right)$ are increasing functions of their variables with saturation. We choose them in the generic form:

$$h\left(u\right)=\frac{a_{1}u}{u+c_1},g\left(p\right)=\frac{p}{p+c_2},$$

where $a_1,c_1$ and $c_2$ are some positive parameters.

We assume that the wealth consumption S is determined by the depreciation (in particular in case of buildings, machinery, etc.) and due to the consumption of the goods and products by the people. Depreciation is considered to be a linear process with a constant rate a. The rate of the individual (per capita) consumption, c, can be described by the Keynes linear consumption function, $c=r+sy$, where y is the per capita income, r and s are some positive coefficients. Assuming additionally that average income is proportional to the wealth, we arrive at the following expression: $S(u,p)=au+(r+su)p$. Altogether, we obtain the following expressions for the functions F and G:

$$F(u,p)=\frac{a_{1}up}{(u+c_1)(p+c_2)}M-(au+(s+ru)p),$$

$$G(u,p)=\alpha p\left(\frac{a_{2}u}{u^2+c_2^2}+{\psi }_0-p\right)-\left({\sigma }_0-\frac{{\sigma }_{1}u}{1+u}\right)p.$$

Please note that in the corresponding non-spatial system with the reaction terms defined as above, $\Omega =\{p>0,-\infty <u<\infty \}$ is an invariant domain but ${\mathbf{R}}_{+}^2$ is not. Negative values of u are interpreted as debt; for example, and further discussion see [24].

In the next section we will discuss how to determine available resources M in the wealth production function.

2.2. Variable Resources

System of Equations (1) and (2) with functions $F(u,p)$, $G(u,p)$ defined above was introduced in [24] under the assumptions that available resources M are constant. In this work we will study the case of variable resources, and we discuss in this section their dependence on wealth taking into account that consumption of resources can be local or nonlocal:

• Local dependence on wealth: $M_1\left(u\right)=(\mu +\theta u)(1-ku)$. The first factor $(\mu +\theta u)$ shows that extraction of resources depends on wealth, and the second factor $(1-ku)$ that resources are consumed proportionally to wealth. Here $\mu$, $\theta$, and k are some positive parameters.
• Nonlocal dependence on wealth: $M_2(u,J\left(u\right))=(\mu +\theta u)(1-kJ\left(u\right))$, where $J\left(u\right)={\int }_{-\infty }^{\infty }\varphi (x-y)u(y,t)dy$, where the kernel $\varphi (x-y)$ shows how the consumption of resources depends on the distance $|x-y|$. The kernel function $\varphi \left(x\right)$ is even and non-negative. It can take into account the cost of transportation and other factors limiting the usage of distant resources.
• Global dependence on wealth: $M_3(u,I\left(u\right))=(\mu +\theta u)(1-kI\left(u\right))$, where $I\left(u\right)={\int }_{-\infty }^{\infty }u(y,t)dy$. In this case, resources are available independently of their location. In particular, this may be the case of intellectual property, software, and other “immaterial” resources, or the cases where the cost of transportation and other distance-related expenses can be neglected.

Let us note that the form of the function $M_1\left(u\right)$ implies that it vanish for $u=1/k$. This means that resources are exhausted for the same level of wealth, which is for the same level of consumption, in spite of the fact that the production of resources increases as a function of wealth. This assumption can be justified if there are several different types of resources such that the production of some of them increases with wealth, while the production of some other remains constant. Hence, there are some limiting resources whose exhaustion occurs for the same level of consumption. For example, agricultural production requires earth and technical equipment. Additional equipment does not increase the level of production if earth is limited, assuming that this equipment does not increase the production per unit area.

In a more general case, we set $M_1\left(u\right)=(\mu +\theta u)(1+\kappa u-ku)$, where the term $1+\kappa u$ is the exhaustion limit, which now depends on wealth. If $\kappa <k$, then this case is reduced to the previous one. The opposite inequality leads the unrealistic situation of the unlimited production and consumption of resources.

In the case of nonlocal consumption of resources, we now have $M_2\left(u\right)=(\mu +\theta u)(1+\kappa u-kJ\left(u\right))$. Contrary to the local consumption of resources, the case $\kappa >0$ cannot be reduced to the case $\kappa =0$, unless it is neglected if $\kappa$ is small enough. A similar remark concerns the global consumption of resources.

3. Wealth Distribution in Excess of Human Resources

We begin the analysis of system (1) and (2) with a special case where the population grows on a much longer time scale, and hence in the dynamics of wealth the available human resources (quantified by variable p) can be considered to be constant. Correspondingly, we consider the equation

$$\frac{\partial u}{\partial t}=D_u\frac{{\partial }^{2}u}{\partial x^2}+F(u,p),$$

(3)

where

$$F(u,p)=W(u,p)-S(u,p),$$

$$W(u,p)=\frac{a_{1}up}{(u+c_1)(p+c_2)}M,S(u,p)=au+(s+ru)p,$$

and p is assumed to be a positive constant. Assuming that $s=0$ and scaling the time and space variables (keeping for them the same notation), we can write Equation (3) in the following form:

$$\frac{\partial u}{\partial t}=D\frac{{\partial }^{2}u}{\partial x^2}+\frac{\theta u(b+u)}{1+k_{1}u}(1+\kappa u-k_2{H}_i)-k_{3}u,$$

(4)

where $i=1,2,3$, the parameters of this equation can be easily expressed through the parameters of Equation (3), and $H_1=u,H_2=J\left(u\right),H_3=I\left(u\right)$. To reduce the number of parameters and to simplify the calculations, we set in what follows $\theta =1,\kappa =0$.

For $i=1$, the nonlinearity

$$F_1\left(u\right)=\frac{u(b+u)}{1+k_{1}u}(1-k_{2}u)-k_{3}u$$

has from one to three non-negative zeros, $u=0$, and up to two positive zeros, $u_{*}$ and $u^{*}$, $u_{*}<u^{*}$. They can be found as roots of a second order polynomial. Clearly, $F_1\left(u\right)<0$ for $u>1/k_2$. We will denote the largest positive solution $u^{*}$. It is stable as solution of the equation $du/dt=F_1\left(u\right)$. If $u=0$ is also a stable solution, then there is an additional unstable solution $u_{*}$. If $u=0$ is unstable, then the positive solution $u^{*}$ is unique.

Considered on the whole axis, Equation (4) can have mono-stable or bistable waves with the limits $u=0$ and $u^{*}$ at infinity, and a pulse solution. The latter is a stationary solution of this equation with the zero limits at infinity. It exists in the bistable case if and only if ${\int }_0^{u^{*}}{F}_1\left(u\right)du>0$, and this solution is unstable. These are well-known results [30], and we will not discuss them in detail here.

3.1. Nonlocal Consumption of Resources

In the case of nonlocal consumption of resources, we consider the equation

$$\frac{\partial u}{\partial t}=D\frac{{\partial }^{2}u}{\partial x^2}+F_2(u,J\left(u\right)),$$

(5)

on the whole axis, where

$$F_2(u,J\left(u\right))=f\left(u\right)(1-k_{2}J\left(u\right))-k_{3}u,$$

$$f\left(u\right)=\frac{u(b+u)}{1+k_{1}u},J\left(u\right)={\int }_{-\infty }^{\infty }\varphi (x-y)u(y,t)dy.$$

Without loss of generality, we can assume that ${\int }_{-\infty }^{\infty }\varphi (x-y)dx={\int }_{-\infty }^{\infty }\varphi (x-y)dy=1$. Then this equation has the same homogeneous-in-space stationary solution $u=0,u_{*},u^{*}$ as Equation (4) in the local case $i=1$.

We will study the stability of the homogeneous-in-space stationary solution $u=u^{*}$. This solution satisfies the equation

$$f\left(u^{*}\right)(1-k_{2}J\left(u^{*}\right))-k_3{u}^{*}=0,(J\left(u^{*}\right)=u^{*}).$$

(6)

Consider the operator linearized about this solution:

$$Lv=D{v}^{″}+f^{\prime }\left(u^{*}\right)(1-k_2{u}^{*})v-k_{2}f\left(u^{*}\right)J\left(v\right)-k_{3}v.$$

Taking into account equality (6), we can write it as follows:

$$Lv=D{v}^{″}+av-k_{2}f\left(u^{*}\right)J\left(v\right),$$

where

$$a=\left(f^{\prime }\left(u^{*}\right)-\frac{f\left(u^{*}\right)}{u^{*}}\right)(1-k_2{u}^{*}).$$

Applying the Fourier transform to the equation $Lv=\lambda v$, we obtain explicit expressions for the eigenvalues:

$$\lambda \left(\xi \right)=-D{\xi }^2+a-k_{2}f\left(u^{*}\right)\tilde{\varphi }\left(\xi \right),\xi \in \mathbb{R},$$

where $\tilde{\varphi }\left(\xi \right)$ is the Fourier transform of the function $\varphi \left(x\right)$.

Consider the function

$$\Phi \left(\xi \right)=a-k_{2}f\left(u^{*}\right)\tilde{\varphi }\left(\xi \right).$$

If it is positive for some $\xi$, then $\lambda \left(\xi \right)$ is also positive for D sufficiently small. Since $u^{*}<1/k_2$, then the sign of a is determined by the expression $f^{\prime }\left(u^{*}\right)-\frac{f\left(u^{*}\right)}{u^{*}}$ that can be positive or negative.

Let us recall that $\Phi \left(0\right)=a-k_{2}f\left(u^{*}\right)<0$. If the function $\varphi \left(x\right)$ is even, then $\tilde{\varphi }\left(\xi \right)={\int }_{-\infty }^{\infty }\varphi \left(x\right)cos\left(\xi x\right)dx<1$ for $\xi \ne 0$. Hence $\Phi \left(\xi \right)>\Phi \left(0\right)$ for $\xi \ne 0$.

Examples

For a special choice of coefficients, $b=k_1=1$, $f\left(u\right)=u$. Then $a=0$, and $\lambda \left(\xi \right)$ can become positive for some $\xi$ only if $\tilde{\varphi }\left(\xi \right)$ becomes positive. If we replace $\varphi \left(x\right)$ by the $\delta$-function (local case), then $\tilde{\varphi }\left(\xi \right)=1$, and $\lambda \left(\xi \right)<0$ for all $\xi$. Set

$${\varphi }_1\left(x\right)=\frac{1}{2N}·\left\{\begin{array}{ccc}1& ,& \left|x\right|\le N\\ 0& ,& \left|x\right|>N\end{array}\right..$$

(7)

Then ${\tilde{\varphi }}_1\left(\xi \right)=sin\left(N\xi \right)/\left(N\xi \right)$. There are positive eigenvalues $\lambda \left(\xi \right)$ if $k_{2}f\left(u^{*}\right)sin\left(N\xi \right)/\left(N\xi \right)>D{\xi }^2$ for some $\xi$. This inequality has solutions if the diffusion coefficient D is sufficiently small. This example is considered in [31] in the analysis of the model of nonlocal consumption of resources in population dynamics. Positive eigenvalues lead to bifurcations of spatial structures.

In another special case $b=k_1=0$, $f\left(u\right)=u^2$, then $a=2u^{*}(1-k_2{u}^{*})$, and the function $\Phi \left(\xi \right)=2u^{*}(1-k_2{u}^{*})-k_2{\left(u^{*}\right)}^2\tilde{\varphi }\left(\xi \right)$ can become positive for a positive $\tilde{\varphi }\left(\xi \right)$. Set

$${\varphi }_2\left(x\right)=\frac{\mu }{2}{e}^{-\mu \left|x\right|}\Rightarrow {\tilde{\varphi }}_2\left(\xi \right)=\frac{{\mu }^2}{{\mu }^2+{\xi }^2}.$$

(8)

The corresponding functions $\lambda \left(\xi \right)$ are shown in Figure 2 (left) for different values of $\mu$. In the limit of large $\mu$ we obtain the local problem for which the solution $u^{*}$ is stable, and $\lambda \left(\xi \right)$ is negative for all real $\xi$. For decreasing $\mu$, nonlocal interaction becomes more essential, $\lambda \left(\xi \right)$ becomes positive for some $\xi$, and the solution $u^{*}$ loses its stability. Let us note that strictly speaking, the curves $\lambda \left(\xi \right)$ correspond to the essential spectrum of the operator L. They become the points of the discrete spectrum (eigenvalues) if we consider the eigenvalue problem on a bounded interval. Figure 2 (right) shows a stationary solution of Equation (5) obtained by direct numerical simulations.

Larger Efficiency of Nonlocal Economy

Let us consider Equation (5) in a bounded interval $0<x<L$ with periodic boundary conditions. The linear stability analysis of the homogeneous-in-space solution $u^{*}$ of Equation (5) shows that it can lose its stability resulting in the emergence of stationary solutions periodic with respect to the space variable x. Consider the diffusion coefficient D as a bifurcation parameter (Figure 3, right). It corresponds to the constant solution for $D>D_c=0.8$. If the diffusion coefficient is less the critical value, then this solution becomes unstable, and a periodic solution bifurcates from it. Denote the corresponding solution by $u_D\left(x\right)$, and determine the total wealth by the formula

$$W\left(D\right)={\int }_0^L{u}_D\left(x\right)dx.$$

Clearly, $W\left(D\right)=L{u}^{*}$ for $D\ge D_c$. Numerical simulations show that $W\left(D\right)$ is a decreasing function with a shallow minimum for some $D<D_c$ (Figure 4, left). Hence, the total wealth of a periodic solution is larger than that for the constant solution, and the total wealth increases when the redistribution of wealth by diffusion decreases. This result can be interpreted as a larger efficiency of nonlocal economy in comparison with the local economy from the point of view of wealth production. A similar result was obtained in [30] (p. 540) for a different model. We will return to this question in the next section for a more complete wealth–population model. It is interesting to note that different stationary regimes can coexist for the same values of parameters (Figure 4, right).

3.2. Global Consumption of Resources

We now consider the global consumption of resources with the equation

$$\frac{\partial u}{\partial t}=D\frac{{\partial }^{2}u}{\partial x^2}+F_3(u,I\left(u\right)),$$

(9)

on the whole axis, where

$$F_3(u,I\left(u\right))=f\left(u\right)(1-k_{2}I\left(u\right))-k_{3}u,$$

$$f\left(u\right)=\frac{u(b+u)}{1+k_{1}u},I\left(u\right)={\int }_{-\infty }^{\infty }u(y,t)dy.$$

The properties of solution in this case are different in comparison with the nonlocal consumption of resources. Instead of periodic solutions, in this case we observe pulse solutions characterized by a single peak of the wealth distribution. Therefore, global economy can lead to the emergence of a single center of wealth accumulation.

Existence and Stability of Pulses

We study the existence of a pulse solution, a positive stationary solution $w\left(x\right)$ of Equation (9) with the zero limits at infinity. We begin with a particular case where the existence of pulses can be studied analytically. We set $b=k_1=0$, and keep the same notation for the integral $I\left(w\right)$:

$$D{w}^{″}+w^2(1-k_{2}I\left(w\right))-k_{3}w=0,w(\pm \infty )=0.$$

(10)

Set $\gamma =1-k_{2}I\left(w\right)$ and consider the equation

$$D{w}^{″}+\gamma w^2-k_{3}w=0.$$

(11)

It is known that it has a positive solution $w_{\gamma }\left(x\right)$ decaying at infinity for any positive $D,\gamma$, and $k_3$. The proof of this assertion follows from the elementary phase plane analysis of the corresponding first-order system of equation. This solution can also be obtained by the integration of the equation. Let $w_{*}\left(x\right)$ be a pulse solution of the equation

$$w^{″}+w^2-w=0.$$

Then it can be directly verified that $w_{\gamma }\left(x\right)=a_1{w}_{*}\left(a_{2}x\right)$ with $a_1=k_3/\gamma$, $a_2=\sqrt{k_3/D}$ is a solution of Equation (10). Then, from the definition of $\gamma$,

$$\gamma =1-k_{2}I\left(w_{\gamma }\right)=1-\frac{\tau }{\gamma },\tau =k_2{k}_3{\int }_{-\infty }^{\infty }{w}_{*}\left(\sqrt{\frac{k_3}{D}}x\right)dx=k_2\sqrt{k_{3}D}{\int }_{-\infty }^{\infty }{w}_{*}\left(y\right)dy.$$

The equation ${\gamma }^2-\gamma +\tau =0$ with respect to $\gamma$ has a solution for $\tau \le 1/4$.

Thus, Equation (10) has two pulse solutions for $\tau <1/4$, a unique solution for $\tau =1/4$, and there are no solutions for $\tau >1/4$. The condition of the existence of solutions implies that the values of the parameters $k_2,k_3$, and D are sufficiently small. In the case of two solutions, numerical simulations show that one of them is stable and another one unstable. The pulse solutions for the single reaction–diffusion equation without nonlocal terms are unstable. Figure 3 shows an example of pulse solutions obtained by numerical simulations of Equation (9).

Bifurcation of Pulses

To analyze the bifurcation of pulses from a homogeneous-in-space solution, let us consider the equation

$$D{w}^{″}+f\left(w\right)(1-k_2{I}_L\left(w\right))-k_{3}w=0$$

(12)

on the bounded interval $0\le x\le L$ with the boundary conditions $w^{\prime }\left(0\right)=w^{\prime }\left(L\right)=0$. Here $I_L\left(w\right)={\int }_0^{L}w\left(x\right)dx$. Homogeneous-in-space solutions of this equation can be found from the equation

$$f\left(w\right)(1-k_{2}Lw)=k_{3}w.$$

(13)

It has two solutions $w_1$ and $w_2$, $w_1<w_2$, if $k_3<k_3^{*}$ for some positive critical value $k_3^{*}$. Considered to be solution of the ordinary differential equation

$$\frac{du}{dt}=F_3(u,I_L\left(u\right)),$$

(14)

$w_1$ is unstable and $w_2$ is stable since

$$f^{\prime }\left(w_2\right)(1-k_{2}L{w}_2)-k_{2}f\left(w_2\right)L-k_3<0.$$

(15)

We study stability of the solution $w_2$ with respect to Equation (9) taking into account spatial distribution. We obtain the following eigenvalue problem:

$$D{v}^{″}+f^{\prime }\left(w_2\right)(1-k_{2}L{w}_2)v-k_{3}v-k_{2}f\left(w_2\right)I_L\left(v\right)=\lambda v,v^{\prime }\left(0\right)=v^{\prime }\left(L\right)=0.$$

(16)

We search solution of this problem in the form $v\left(x\right)=cos(\pi nx/L)$, $n=0,1,2,\dots$ Then

$${\lambda }_0=f^{\prime }\left(w_2\right)(1-k_{2}L{w}_2)-k_{2}f\left(w_2\right)L-k_3,$$

$${\lambda }_n=f^{\prime }\left(w_2\right)(1-k_{2}L{w}_2)-k_3-D{\left(\pi n\right)}^2/L^2,n=1,2,\dots$$

By virtue of inequality (15), ${\lambda }_0<0$. Taking into account (13), the eigenvalue ${\lambda }_1$ can be written as follows:

$${\lambda }_1=\left(f^{\prime }\left(w_2\right)-\frac{f\left(w_2\right)}{w_2}\right)(1-k_{2}L{w}_2)-D{\pi }^2/L^2\equiv \rho -D{\pi }^2/L^2,$$

where $\rho$ denotes the first term in the right-hand side.

If the function $f\left(w\right)$ is convex, then $\rho >0$, and ${\lambda }_1$ can be positive or negative depending on D and L. If ${\lambda }_1>0$, then the homogeneous-in-space solution $w_2$ loses its stability with respect to spatial perturbation. The eigenfunction $v_1\left(x\right)=cos(\pi x/L)$ corresponds to a half-periodic solution. It can be extended to a periodic solution in the case of a fixed L. If we increase L, we obtain in the limit a single pulse on the whole axis, as discussed above.

Let us discuss the role of the integral term in the emergence of pulses. If we replace $I\left(u\right)$ in (9) by $Lu$, then the stationary solutions $w_1$ and $w_2$ do not change, and the eigenvalue ${\lambda }_0$ of the linearized problem remains also the same. However, all other eigenvalues

$${\tilde{\lambda }}_n=f^{\prime }\left(w_2\right)(1-k_{2}L{w}_2)-k_{2}f\left(w_2\right)L-k_3-D{\left(\pi n\right)}^2/L^2,n=1,2,\dots$$

are different, and they are negative together with the principal eigenvalue. Hence, the integral term does not change the principal eigenvalue, which is negative, but it can make the second eigenvalue positive leading to the bifurcation of a pulse solution.

4. Periodic Structures and Pulses for the Wealth-Population System

4.1. Periodic Structures for Nonlocal Consumption

Now we will study the system of two equations

$$\frac{\partial u}{\partial t}=D_u\frac{{\partial }^{2}u}{\partial x^2}+f\left(u\right)p(1-k_{2}J\left(u\right))-k_{3}u,$$

(17)

$$\frac{\partial p}{\partial t}=D_p\frac{{\partial }^{2}p}{\partial x^2}+\alpha p(\psi \left(u\right)-p)-{\sigma }_{0}p$$

(18)

on the whole real axis, where

$$f\left(u\right)=\frac{u(b+u)}{1+k_{1}u},\psi \left(u\right)=\frac{a_{2}u}{u^2+c_2^2}+\frac{{\sigma }_{1}u}{\alpha (1+u)}+{\psi }_0.$$

It is similar to system (1) and (2) with the only simplification that the wealth production in Equation (17) is considered to be a linear function of the population density p. This approximation is justified if $c_2$ is large enough.

The homogeneous-in-space stationary solutions of this system, different from $u=0$, $p=0$ can be found from the equations:

$$f\left(u\right)(\psi \left(u\right)-{\sigma }_0/\alpha )(1-k_{2}u)=k_{3}u,$$

(19)

and $p=\psi \left(u\right)-{\sigma }_0/\alpha$. As before, we assume here that ${\int }_{-\infty }^{\infty }\varphi \left(x\right)dx=1$. Equation (19) has a positive solution if $k_3$ is less than some critical value $k_3^{*}$. Denote by $u_{*}$ the maximal solution and assume that it is simple, and set $p_{*}=\psi \left(u_{*}\right)-{\sigma }_0/\alpha$. We suppose that this solution is stable with respect to the ODE system

$$\frac{du}{dt}=f\left(u\right)p(1-k_{2}u)-k_{3}u,$$

(20)

$$\frac{dp}{dt}=\alpha p(\psi \left(u\right)-p)-{\sigma }_{0}p.$$

(21)

Therefore, the eigenvalues of the linearized matrix

$$A_0=\left(\begin{array}{ccc}{f}^{\prime }\left(u_{*}\right)p_{*}(1-k_2{u}_{*})-k_{2}f\left(u_{*}\right)p_{*}-k_3& & f\left(u_{*}\right)(1-k_2{u}_{*})\\ \alpha {\psi }^{\prime }\left(u_{*}\right)p_{*}& & -\alpha p_{*}\end{array}\right)$$

have negative real parts.

To study the stability of the solution $(u_{*},p_{*})$ with respect to system (17) and (18), we consider the eigenvalue problem

$$D_u{v}^{″}+\left(f^{\prime }\left(u_{*}\right)p_{*}(1-k_2{u}_{*})-k_3\right)v-k_{2}f\left(u_{*}\right)p_{*}J\left(v\right)+f\left(u_{*}\right)(1-k_2{u}_{*})z=\lambda v,$$

(22)

$$D_p{z}^{″}+\alpha p_{*}{\psi }^{\prime }\left(u_{*}\right)v-\alpha p_{*}z=\lambda z.$$

(23)

The eigenfunction with the components $(v,z)$ corresponds to the perturbation of the solution $(u_{*},p_{*})$. Here we use that $J\left(u_{*}\right)=u_{*}$. Applying the Fourier transform, we obtain the algebraic system of equations

$$A\left(\xi \right)\left(\begin{array}{c}\tilde{v}\\ \tilde{z}\end{array}\right)=\lambda \left(\begin{array}{c}\tilde{v}\\ \tilde{z}\end{array}\right),$$

where tilde denotes the Fourier transform, and

$$A\left(\xi \right)=\left(\begin{array}{ccc}{f}^{\prime }\left(u_{*}\right)p_{*}(1-k_2{u}_{*})-k_{2}f\left(u_{*}\right)p_{*}\tilde{\varphi }\left(\xi \right)-k_3-D_u{\xi }^2& & f\left(u_{*}\right)(1-k_2{u}_{*})\\ \alpha {\psi }^{\prime }\left(u_{*}\right)p_{*}& & -\alpha p_{*}-D_v{\xi }^2\end{array}\right).$$

Here $\tilde{\varphi }\left(\xi \right)$ is the Fourier transform of the kernel $\varphi \left(x\right)$ of the integral $J\left(v\right)$. Let us note that $A\left(0\right)=A_0$. Therefore, the eigenvalues of this matrix have negative real parts. We search the values of parameters for which the matrix $A\left(\xi \right)$ has a positive eigenvalue for some $\xi$. The critical values of parameters (dispersion relation) is determined by the equality $detA\left(\xi \right)=0$. Hence, we obtain the following equation:

$$D_u{D}_v{\xi }^4-T\left(A_0\right){\xi }^2+D\left(A_0\right)=k_{2}f\left(u_{*}\right)p_{*}(1-\tilde{\varphi }\left(\xi \right))(\alpha p_{*}+D_v{\xi }^2),$$

(24)

where $T\left(A_0\right)=D_v{a}_{11}+D_u{a}_{22}$, $a_{11}$ and $a_{22}$ are the diagonal elements, and $D\left(A_0\right)$ is the determinant of the matrix $A_0$. We note that $T\left(A_0\right)<0$ and $D\left(A_0\right)>0$. Therefore, the left-hand side of this equation is positive for all real $\xi$. Graphical solution of this equation is shown in Figure 5 (left).

The direct analysis of this equation is cumbersome, and we begin with a particular case where ${\psi }^{\prime }\left(u_{*}\right)=0$. This condition is satisfied if $\psi \left(u\right)$ is identically constant. Then the matrix $A\left(\xi \right)$ has two eigenvalues,

$${\lambda }_1\left(\xi \right)=\left(f^{\prime }\left(u_{*}\right)-\frac{f\left(u_{*}\right)}{u_{*}}\right)p_{*}(1-k_2{u}_{*})-k_{2}f\left(u_{*}\right)p_{*}\tilde{\varphi }\left(\xi \right)-D_u{\xi }^2$$

(25)

and a negative eigenvalue ${\lambda }_2\left(\xi \right)=-(\alpha p_{*}+D_v{\xi }^2)$. The eigenvalue ${\lambda }_1\left(\xi \right)$ is similar to the eigenvalue $\lambda \left(\xi \right)$ in Section 2.1, and the conditions which determine the equality ${\lambda }_1\left(\xi \right)=0$ are also similar.

Example of numerical simulation of system (17) and (18) is shown in Figure 6. Each horizontal cross section of this figure corresponds to a stationary solution with a given value of $D_u$. The stationary solutions are periodic in space, and the population distribution has a minimum at the maximum of wealth distribution because the function $\psi \left(u\right)$ decreases for u sufficiently large (Figure 5). It is also interesting to note that decrease of the parameter $D_u$ increases the variation of wealth and the total wealth, while the total population decreases. If $D_p$ decreases, then both total wealth and total population increase.

4.2. Pulses for Global Consumption

In this section, we study the case of global consumption for the system of equations

$$\frac{\partial u}{\partial t}=D_u\frac{{\partial }^{2}u}{\partial x^2}+f\left(u\right)p(1-k_{2}I\left(u\right))-k_{3}u,$$

(26)

$$\frac{\partial p}{\partial t}=D_p\frac{{\partial }^{2}p}{\partial x^2}+\alpha p(\psi \left(u\right)-p)-{\sigma }_{0}p,$$

(27)

where, as before,

$$f\left(u\right)=\frac{u(b+u)}{1+k_{1}u},\psi \left(u\right)=\frac{a_{2}u}{u^2+c_2^2}+\frac{{\sigma }_{1}u}{\alpha (1+u)}.$$

This system is similar to system (17) and (18) where the integral $J\left(u\right)$ is replaced by $I\left(u\right)$. We consider it on a bounded interval $0<x<L$ with the Neumann boundary condition. In the stationary case, we have the following boundary value problem:

$$D_u{u}^{″}+f\left(u\right)p(1-k_{2}I\left(u\right))-k_{3}u=0,$$

(28)

$$D_p{p}^{″}+\alpha p(\psi \left(u\right)-p)-{\sigma }_{0}p=0,$$

(29)

$$u^{\prime }\left(0\right)=u^{\prime }\left(L\right)=0,p^{\prime }\left(0\right)=p^{\prime }\left(L\right)=0,$$

(30)

$I\left(u\right)={\int }_0^{L}u\left(x\right)dx$ (the subscript L in the integral is omitted for brevity). Similar to the nonlocal case considered in the previous section, we consider the maximal positive constant solution $(u_{*},p_{*})$ and suppose that it is stable with respect to the corresponding ODE system. This solution satisfies the equations

$$f\left(u\right)p(1-k_{2}Lu)-k_{3}u=0,\alpha (\psi \left(u\right)-p)-{\sigma }_0=0,$$

(31)

where we take into account that $I\left(u\right)=Lu$.

Linearizing problem (28)–(30) about the solution $(u_{*},p_{*})$, we obtain the eigenvalue problem:

$$D_u{v}^{″}+f^{\prime }\left(u_{*}\right)p_{*}(1-k_{2}L{u}_{*})v-k_{3}v-k_{2}f\left(u_{*}\right)p_{*}I\left(v\right)+f\left(u_{*}\right)(1-k_{2}L{u}_{*})z=\lambda v,$$

(32)

$$D_p{z}^{″}+\alpha p_{*}{\psi }^{\prime }\left(u_{*}\right)v-\alpha p_{*}z=\lambda z,$$

(33)

$$v^{\prime }\left(0\right)=v^{\prime }\left(L\right)=0,z^{\prime }\left(0\right)=z^{\prime }\left(L\right)=0.$$

(34)

We search a solution of this problem in the form $v\left(x\right)=c_{1}cos(\pi nx/L)$, $z\left(x\right)=c_{2}cos(\pi nx/L)$, $n=0,1,2,\dots$ Then we obtain the linear algebraic system of equations:

$$B\left(\begin{array}{c}{c}_1\\ c_2\end{array}\right)=\lambda \left(\begin{array}{c}{c}_1\\ c_2\end{array}\right),$$

where

$$B_0=\left(\begin{array}{ccc}{f}^{\prime }\left(u_{*}\right)p_{*}(1-k_{2}L{u}_{*})-k_3-k_{2}Lf\left(u_{*}\right)p_{*}& & f\left(u_{*}\right)(1-k_{2}L{u}_{*})\\ \alpha p_{*}{\psi }^{\prime }\left(u_{*}\right)& & -\alpha p_{*}\end{array}\right),$$

$$B_n=\left(\begin{array}{ccc}{f}^{\prime }\left(u_{*}\right)p_{*}(1-k_{2}L{u}_{*})-k_3-D_u{\left(\pi n\right)}^2/L^2& & f\left(u_{*}\right)(1-k_{2}L{u}_{*})\\ \alpha p_{*}{\psi }^{\prime }\left(u_{*}\right)& & -\alpha p_{*}-D_p{\left(\pi n\right)}^2/L^2\end{array}\right),$$

$n=1,2,\dots$ Since $(u_{*},p_{*})$ is a stable solution of the ODE system, then the eigenvalues of the matrix $B_0$ have negative real parts. If one of the eigenvalues of the matrix $B_1$ is positive, then this solution can lose its stability resulting in the bifurcation of a space-dependent solution. Denote the elements of this matrix by $b_{ij}$, $i,j=1,2$. Since $b_{12},b_{21}>0$ and $b_{22}<0$, then $det{B}_1<0$ if $b_{11}>0$. Hence, inequality

$$f^{\prime }\left(u_{*}\right)p_{*}(1-k_{2}L{u}_{*})-k_3-D_u{\pi }^2/L^2>0$$

(35)

provides a sufficient condition for the existence of a positive eigenvalue of the matrix $B_1$. This condition is similar to the condition obtained in Section 2.2 for the single equation. It is satisfied for $D_u$ sufficiently small or L sufficiently large if

$$f^{\prime }\left(u_{*}\right)p_{*}(1-k_{2}L{u}_{*})-k_3=\left(f^{\prime }\left(u_{*}\right)-\frac{f\left(u_{*}\right)}{u_{*}}\right)p_{*}(1-k_{2}L{u}_{*})>0$$

(36)

(see (31)). Condition (35) is not satisfied if the diffusion coefficient $D_u$ is sufficiently large. Hence, diffusion stabilizes the constant solution.

If this condition is not satisfied, $det{B}_1$ can be still negative. However, if $b_{11}<0$ and $det{B}_1>0$ for $D_u=D_p=0$, then the determinant remains positive in the case of any positive diffusion coefficients. Thus, the origin of the instability is different in comparison with the dissipative (Turing) instability.

Examples

If $f\left(u\right)=u$, then condition (36) is not satisfied. If $f\left(u\right)=u^2$, then Equation (31) have a solution for $k_3$ sufficiently small. In this case, condition (36) is satisfied. An example of wealth and population distributions is shown in Figure 7. Since the function $\psi \left(u\right)$ decreases for u sufficiently large, then the population distribution has a minimum at the maximum of wealth distribution. Dependence of stationary solutions on the coefficient $D_u$ is represented in Figure 8 where each horizontal cross section shows the levels of the corresponding solution. The stationary solution is homogeneous-in-space for large values of $D_u$. Decrease of this parameter leads to the emergence of a non-homogeneous distribution (Figure 7) accompanied by the increase of the total wealth and the decrease of the total population.

5. Discussion and Conclusions

5.1. On the Mechanisms of Pattern Formation

Turing Instability

There are numerous applications where a constant solution is stable as a solution of the corresponding ODE system but it loses its stability with respect to spatial perturbations. The most well-known mechanism of such instability is suggested by A. Turing. It is applicable for reaction–diffusion systems of two or more equations. From the mathematical point of view, this instability can be explained in terms of matrices corresponding to the eigenvalue problem obtained due to linearization about the constant solution. Consider the matrix $A\left(\xi \right)=A_0-D{\xi }^2$, where the matrix $A_0$ corresponds to the reaction part of the system, D is the diagonal diffusion matrix with positive diagonal elements, and $\xi$ is the wavenumber of the spatial perturbation. Assuming that all eigenvalues of the matrix $A_0$ have negative real parts, what can be concluded about the eigenvalues of the matrix $A\left(\xi \right)$? Clearly, if all diagonal elements of the matrix D are equal to each other, then the eigenvalues of the matrix $A\left(\xi \right)$ have negative real parts for any real $\xi$. A similar conclusion holds for matrices $A_0$ with non-negative off-diagonal elements and any diagonal matrix D. However, if off-diagonal elements of the matrix $A_0$ have variable signs and diagonal elements of the matrix D are different from each other, then the matrix $A\left(\xi \right)$ can have eigenvalues with positive real parts for some values of $\xi$. This counterintuitive results is the mathematical basis of the Turing (or dissipative, diffusive) instability. From the physical point of view, this means that diffusion can destabilize a homogeneous-in-space solution, which is stable without diffusion. Turing instability is widely studied in relation to various chemical and biological processes.

Nonlocal Consumption of Resources

Similar mechanisms of pattern formation can manifest themselves for nonlocal reaction–diffusion equations. In this work we consider two of such mechanisms. We will explain them with a model example where the location of the spectrum of the linearized problem is given by the formula

$$\lambda \left(\xi \right)=-D{\xi }^2+a-\tilde{\varphi }\left(\xi \right),$$

which represents a simplified version of the spectral curves considered in this work. The first term in the right-hand side of this expression appears in the Fourier transform of the diffusion term, a constant a remains from the reaction term, and $\tilde{\varphi }\left(\xi \right)$ is the Fourier transform of the integral kernel $\varphi \left(x\right)$. To make the calculations even more explicit, we set $\tilde{\varphi }\left(\xi \right)=sin\left(N\xi \right)/\left(N\xi \right)$ for the function ${\varphi }_1\left(x\right)$ given by (7). Condition $\lambda \left(0\right)=a-1<0$ implies the stability of solution without spatial perturbations. Under which conditions $\lambda \left(\xi \right)$ can become positive for some real $\xi \ne 0$?

Denote by ${\xi }_0$ the point of minimum of the function $\tilde{\varphi }\left(\xi \right)$, $\tilde{\varphi }\left({\xi }_0\right)<0$ and consider the value $\lambda \left({\xi }_0\right)=-D{\xi }_0^2+a-\tilde{\varphi }\left({\xi }_0\right)$. Since $a-\tilde{\varphi }\left({\xi }_0\right)>a-1$, then $\lambda \left({\xi }_0\right)$ can be positive even if $a-1<0$, and the solution becomes unstable with respect to spatial perturbations. Such situation is realized in the case of nonlocal consumption of resources in Section 2.1 and Section 3.1.

The value ${\xi }_0$ determines the spatial periodicity of the perturbation. If we consider the problem on a bounded interval, then the real variable $\xi$ takes a discrete set of values ${\xi }_k$ determined by the length of the interval. The frequency ${\xi }_k$ closest to ${\xi }_0$ gives the perturbation with the fastest growth.

Global Consumption of Resources

Suppose that

$$N{\xi }_k=\pi k,k=1,2,\dots ,$$

(37)

where N is the same as in the previous paragraph (see also Equation (7)). Then $\varphi \left({\xi }_k\right)=0$, $\lambda \left({\xi }_k\right)=-D{\xi }_k^2+a$, and the maximal eigenvalue $\lambda \left({\xi }_1\right)=-D{\pi }^2+a$ can be positive if $a>0$ and D is sufficiently small. Since ${\xi }_k=\pi k/L$ to satisfy the boundary conditions, then condition (37) is satisfied for $N=L$, i.e., in the case of global consumption. In this case, the integral does not change stability without spatial perturbations since $\lambda \left(0\right)=a-1<0$ but it can change the stability with respect to the spatial perturbations since the integral vanishes on all spatial modes.

There are two differences here in comparison with the previous case: the instability can occur only for $a>0$, and the fastest growing mode corresponds to $k=1$. Therefore, this instability leads to the emergence of pulses contrary to the periodic structures in the nonlocal case.

5.2. Properties of the Nonlocal Economy

The dynamics of human populations are strongly influenced by the wealth distribution through several ways. Increasing wealth can increase the natality, e.g., due to an increase in the household income and hence larger available resources to bring children, but it can also decrease the natality due to a cultural shift; it can decrease the mortality rate provided better health care; and a non-uniform wealth distribution can influence directions and intensity of human migration. Therefore, all these factors should be considered in their interaction. In our previous work, we introduced a novel coupled wealth–population model [23], the properties of its non-spatial counterpart studied in [32] and the emergence of spatial patterns due to the Turing instability analyzed in [24].

In this paper, we further develop wealth–population models by taking into account the dynamic nature of the available resources. Note that by resources we understand not only natural resources such as gas, oil, minerals, etc., but also the production components that results from human activity, e.g., equipment and knowledge. The latter becomes particularly important for modern high-tech products where the intellectual component (software, patents, and so on) can constitute the largest part of the cost of the final product. Moreover, intellectual resources, contrary to natural resources, do not have fixed geographical location. They represent global resources in the sense that the access to them is not spatially dependent, and the cost of their transportation is not a limiting factor for the global economy. Hence, we need to take into account nonlocal effects in the consumption of resources. This is the main objective of this work, and we observe that nonlocal and global economies have quite different properties in comparison with traditional local economies.

An important characteristic of the economy with nonlocal consumption of resources, where the access to resources remains space-dependent but can be sufficiently extended, is the emergence of non-uniform in space wealth distribution. Such structures arise due to the competition for resources if wealth redistribution (diffusion) is sufficiently small. The total wealth in the case of a non-uniform distribution is larger than in the uniform case signifying that nonlocal economy is more efficient than the local one.

In the case of global consumption of resources, where the access to resources is space-independent, a non-uniform distribution can also emerge. Contrary to the previous case, it manifests itself as a single peak in wealth distribution and not as a periodic structure. Hence, global economy can be characterized by a single place of wealth concentration.

5.3. Conclusions and Future Work

Distribution of wealth in space is known to be strongly heterogeneous and this occurs on different spatial scales, ranging from the small scale of separate cities, boroughs, towns and villages to the large scale of states, countries and continents. Although it is not our goal in this paper to provide any quantitative comparison between data and theory, we have been motivated by a few examples where the distribution of wealth across space (e.g., see Figure 1, right) and the distribution of population [24] exhibits features that can be regarded as signs of periodicity. In particular, in North America this suspected periodicity in both wealth and population distribution apparently occurs on the same spatial scale of 1500–2000 km (cf. Figure 1 (right) and Figure 3a in [24]).

In this paper, (see also [24]), we have shown that the emergence of a pattern in the spatial distributions of population and wealth is an inherent property of the corresponding dynamical system. Arguably, this is in a quantitative agreement with some available data (see the previous paragraph). Any quantitative comparison between the theory and the data would require the knowledge of realistic parameter values, which is currently not available; this should become a focus of future research.

We mention here that our model is conceptual and, as such, omits many factors and processes that can, in principle, affect the mechanisms of pattern formation. Transportation costs is one such factor [11]. Another one is the density-dependent diffusion and/or cross-diffusion [23], to allow the population to migrate along the wealth gradient rather than randomly. As the wealth distribution on the small scale is apparently stochastic (cf. Figure 1, left), an inclusion into the model an explicit stochasticity is yet another mathematical challenge. These will become a focus of future research.