An economic cross-diffusion mutualistic model for cities emergence

TL;DR

This paper models city emergence through an economic cross-diffusion system with mutualistic interactions, analyzing stability and pattern formation, supported by numerical simulations.

Contribution

It introduces a novel cross-diffusion model for urban development based on profit-driven mutualism, with stability analysis and numerical validation.

Findings

Conditions for pattern formation identified

Cities modeled as high and low concentration regions

Numerical simulations confirm theoretical predictions

Abstract

We study an evolution cross-diffusion problem with mutualistic Lotka-Volterra reaction term to modelize the long-term spatial distribution of labor and capital. The mutualistic behavior is deduced from the gradient flow associated to profits maximization. We perform a linear and weakly nonlinear stability analysis and find conditions under which the uniform optimum of profits becomes unstable, leading to pattern formation. The patterns alternate regions of high and low concentrations of both labor and capital, which may be interpreted as cities. Finally, numerical simulations based on the weakly nonlinear analysis, as well as in a finite element approximation, are provided.

Tables

Table 1. Table 1: Parameters of the CES production function, resulting (scaled) Lotka-Volterra coefficients, rescaled diffusion coefficients, equilibrium, critical bifurcation and wave numbers, time to stationary state and mean square error between the FEM solution and the WNL solution.

	Experiments
	1	2	3	4
$w^{\ast },r^{\ast }$	1, 0.3	0.4, 0.3	0.95, 0.95	0.4, 0.5
$\alpha ,\beta$	0.3, 0.6	0.3, 0.6	0.29, 0.3	0.3, 0.6
$ϵ,\eta$	0.5, 0.2	0.5, 0.2	0.75, 0.1	0.5, 0.2
${\alpha }_1,{\alpha }_2$	0.5, 0.15	0.2, 0.15	0.24, 0.24	0.05, 0.15
${\beta }_1,{\beta }_2$	2.35, 2.47	0.48, 8.38	0.66, 1.97	4.5e-2, 5.4
$c_1,c_2$	0.01, 0.01	0.01, 0.01	0.1, 0.1	0.01, 0.01
$a_1,a_2$	0.3, 3e-4	7e-3, 0.01	0.41, 2.4	3.72, 2.7e-5
$(L^{\ast },K^{\ast })$	0.29, 0.18	0.6, 0.09	2.35, 1.31	1.96, 0.03
$b_c,k_c$	1.56, 3.99	0.42, 6	1.14, 4.02	205.9, 3.51
$T_s,\text{MSE}$	1557, 1.7e-2	1630, 1.4e-2	240, 4.7e-2	5000

Equations

$$\displaystyle\partial_{t}L-\operatorname{div}\big{(}(c_{1}+a_{11}L)\nabla L-a_{12}L\nabla K\big{)}=L(\alpha_{1}-\beta_{11}L+\beta_{12}K)$$

$$\displaystyle\partial_{t}K-\operatorname{div}\Big{(}\big{(}c_{2}-a_{22}\frac{K}{K_{s}^{2}+K^{2}}\big{)}\nabla K\Big{)}=K(\alpha_{2}+\beta_{21}L-\beta_{22}K)$$

$$\nabla L\cdot \nu =\nabla K\cdot \nu =0$$

$$L(0,\cdot )=L_0,K(0,\cdot )=K_0$$

$$\Pi (K,L)=pY(K,L)-wL-rK,$$

$${\partial }_{L}Y(K,L)-w=0,{\partial }_{K}Y(K,L)-r=0.$$

$$\displaystyle\partial_{t}L=L\Big{(}\tilde{\alpha}_{1}+\tilde{\beta}_{11}L+\tilde{\beta}_{12}K\Big{)},$$

$$\displaystyle\partial_{t}K=K\Big{(}\tilde{\alpha}_{2}+\tilde{\beta}_{21}L+\tilde{\beta}_{22}K\Big{)},$$

$$\tilde{B}:=({\tilde{\beta }}_{ij})=H({\Pi }_e^{\ast }),({\tilde{\alpha }}_1,{\tilde{\alpha }}_2{)}^t=-\tilde{B}(L_e^{\ast },K_e^{\ast }{)}^t.$$

$${\tilde{\beta }}_{11}<0,{\tilde{\beta }}_{22}<0,\text{and }{\tilde{\beta }}_{12}={\tilde{\beta }}_{21}>0,$$

$$\mathrm{tr}(J\mathbf{G}(\mathbf{w}))=\sum _{i=1}^2{w}_i{\partial }_i{\mathbf{F}}_i(\mathbf{w}),\det (J\mathbf{G}(\mathbf{w}))=w_1{w}_2\det (J\mathbf{F}(\mathbf{w})),$$

$$Y(K,L)=A{\left(\alpha K^{\eta }+\beta L^{\eta }\right)}^{\frac{ϵ}{\eta }},$$

$$\alpha ,\beta ,ϵ\in (0,1),\eta <1,\eta \ne 0,\text{and }ϵ\ge \eta .$$

$${\partial }_{L}Y(K,L)=Aϵ[\alpha K^{\eta }+\beta L^{\eta }{]}^{\frac{ϵ}{\eta }-1}\beta L^{\eta -1},$$

$${\partial }_{K}Y(K,L)=Aϵ[\alpha K^{\eta }+\beta L^{\eta }{]}^{\frac{ϵ}{\eta }-1}\alpha K^{\eta -1},$$

$$H(\Pi )=\left(\begin{array}{cc}\frac{\beta (ϵ-\eta )Q^{-1}w{L}^{\eta }-(1-\eta )w}{L}& \alpha (ϵ-\eta )Q^{-1}w{K}^{\eta -1}\\ \beta (ϵ-\eta )Q^{-1}r{L}^{\eta -1}& \frac{\alpha (ϵ-\eta )Q^{-1}r{K}^{\eta }-(1-\eta )r}{K}\end{array}\right),$$

$$\displaystyle(L_{e}^{*})^{1-\varepsilon}=\frac{\rho}{w^{*}}\beta^{\epsilon/\eta},\quad(K_{e}^{*})^{1-\varepsilon}=\frac{\rho}{r^{*}}\alpha^{\epsilon/\eta},\quad\rho=A\epsilon\Big{(}1+\big{(}\frac{w^{*}}{r^{*}}\big{)}^{\frac{\eta}{1-\eta}}\big{(}\frac{\alpha}{\beta}\big{)}^{\frac{1}{1-\eta}}\Big{)}^{\frac{\epsilon-\eta}{\eta}}.$$

$${\tilde{\alpha }}_1>0,{\tilde{\alpha }}_2>0,{\tilde{\beta }}_{11}<0,{\tilde{\beta }}_{22}<0,{\tilde{\beta }}_{12}={\tilde{\beta }}_{21}\ge 0,\det (\tilde{B})>0.$$

$$K_e^{\ast }={\left(\frac{r^{\ast }\beta }{w^{\ast }\alpha }\right)}^{\frac{1}{\eta -1}}{L}_e^{\ast }.$$

$$\det (\tilde{B})=(1-ϵ)(1-\eta )\frac{w^{\ast }{r}^{\ast }}{L_e^{\ast }{K}_e^{\ast }}>0.$$

$$\displaystyle\alpha\beta\frac{(L_{e}^{*})^{\eta}(K_{e}^{*})^{\eta}}{(\alpha(L_{e}^{*})^{\eta}+\beta(K_{e}^{*})^{\eta})^{2}}>\Big{(}\frac{1-\eta}{\epsilon-\eta}\Big{)}^{2}.$$

$$\bar{L}={\beta }_{21}L,\bar{K}={\beta }_{12}K,{\bar{L}}_0={\beta }_{21}{L}_0,{\bar{K}}_0={\beta }_{12}{K}_0,{\bar{K}}_s={\beta }_{12}{K}_s,$$

$$a_1=\frac{a_{11}}{{\beta }_{21}},a_2=a_{22}{\beta }_{12},b=\frac{a_{12}}{{\beta }_{12}},{\beta }_1=\frac{{\beta }_{11}}{{\beta }_{21}},{\beta }_2=\frac{{\beta }_{22}}{{\beta }_{12}},$$

$$\bar{t}=\frac{t}{\gamma },\bar{x}=\frac{x}{\sqrt{\gamma }},$$

$$\displaystyle\partial_{t}L-\partial_{x}\big{(}(c_{1}+a_{1}L)\partial_{x}L-bL\partial_{x}K\big{)}=\gamma L(\alpha_{1}-\beta_{1}L+K)$$

$$\displaystyle\partial_{t}K-\partial_{x}\Big{(}\big{(}c_{2}-a_{2}g(K)\big{)}\partial_{x}K\Big{)}=\gamma K(\alpha_{2}+L-\beta_{2}K)$$

$${\partial }_{x}L\cdot \nu ={\partial }_{x}K\cdot \nu =0$$

$$L(0,\cdot )=L_0,K(0,\cdot )=K_0$$

$$c_1+a_1>0,c_2>a_{2}g(K_s).$$

$$\displaystyle(L^{*},K^{*})=\Big{(}\frac{\alpha_{2}+\alpha_{1}\beta_{2}}{\beta_{1}\beta_{2}-1},\frac{\alpha_{1}+\alpha_{2}\beta_{1}}{\beta_{1}\beta_{2}-1}\Big{)}.$$

$${\partial }_t\mathbf{w}-Q\Delta \mathbf{w}=\gamma R\mathbf{w}$$

$$\nabla u_1\cdot \nu =\nabla u_2\cdot \nu =0$$

$$u_1(0,\cdot )=\tilde{L}(0,\cdot )-L_0,u_2(0,\cdot )=\tilde{K}(0,\cdot )-K_0$$

$$R=\left(\begin{array}{cc}-{\beta }_1{L}^{\ast }& L^{\ast }\\ K^{\ast }& -{\beta }_2{K}^{\ast }\end{array}\right),Q=\left(\begin{array}{cc}{c}_1+a_1{L}^{\ast }& -b{L}^{\ast }\\ 0& c_2-a_{2}g(K^{\ast })\end{array}\right).$$

$$A_k\mathbf{w}=\lambda \mathbf{w},\text{for }{A}_k=\gamma R-k^{2}Q\text{ and }k=|\mathbf{k}|.$$

$$\displaystyle\lambda_{jk}=\frac{1}{2}\Big{(}\operatorname{tr}(A_{k})+(-1)^{j}\sqrt{\operatorname{tr}(A_{k})^{2}-4\det(A_{k})}\Big{)},\quad j=1,2,$$

Full text

*Keywords: * Cross-diffusion, mutualism, Turing instability, weakly nonlinear analysis, labor, capital, city.

1 Introduction

The first human settlements that can be regarded as cities had one element in common: the presence of a natural formation that endowed the location with a characteristic making small amounts of capital particularly productive. This is, for instance, the case of the alluvial plains of Lower Mesopotamia, where a natural levee made it very easy to irrigate large extensions of land with few tools. Moreover, the recent invention of the plow led to an extraordinary gain in productivity.

In this article, we study a mathematical model to explore the type of relation between labor and capital that leads to the formation of cities, understood as accumulation of labor and capital. Unlike other models, where this interaction is assumed to take place in a discrete spatial domain formed by patches of employment [2], we consider a continuous time-space model where the relations between labor and capital leading to global growth or decay are derived solely from a production function, and where the spatial relocation is induced by attraction-repulsion mechanisms among the units of labor and capital and, secondarily, by random motion.

Current economic geography models like those introduced by Krugman [17] or by Tabuchi and Thisse [20] are focused in modelling industry concentration or in explaining the emergence of central places. These models are set in a discrete spatial domain, being the main forces explaining concentration the decreasing transportation costs and the increasing returns to scale in the production of manufactured goods. A production function is said to have increasing returns to scale when the increase of one unit of all the inputs, induce an increase of output larger than one. See [18] for an extended discussion about returns to scale in a production function similar to ours.

In the model we introduce, cities emergence is a possible outcome of the action of general interacting forces that relocate labor and capital, and which is tidily connected to growth. Thus, we do not need to specify the actual forces involved in the process, that can include those of Krugman or Tabuchi and Thisse. Furthermore, we do not make use of the assumption of increasing return to scale in the production of goods.

Our model is related to that introduced by Volpert, Petrovskii, and Zincenko [22, 23] for the interaction of human migration (labor, in our model) and wealth distribution (capital). Both models are formulated in terms of reaction-diffusion partial differential equations in which the reaction term establishes the growth-decline relationship between humans and capital, and where the diffusion term stands for the spatial relocation of both densities. However, there are important differences between both models, being perhaps the most important the role played by the cross diffusion mechanism, which is weak in Volpert’s model while crucial in ours. Besides, Volpert’s analysis is focused in the existence of travelling fronts, this is, in the existence of particular solutions of the problem, while our focus is in Turing’s bifurcation analysis.

1.1 A historical example

Before moving into the description of the mathematical model, let us motivate it by examining the economic events that happened five thousand years ago in the levees of Mesopotamia [8].

Previously to the arising of the first civilizations, the Humanity, for hundreds of thousands of years, was nearly uniformly distributed in geographical areas where the environmental conditions were favorable. This uniform state suffered a perturbation caused by the Neolithic revolution. In the levees of Mesopotamia, people concentrated along the natural ridge because it was easy to irrigate, and with the use of plows and animal force, people found their labor very productive. There was growth.

The increase of productivity attracted more people who relocated in longer stretches of the ridge producing more and more until the ridge was crowded. Because, the ridge was on limited supply, more workers and more capital faced decreasing returns to scale due to competence.

At some moment, the extra capital and the extra workers were repelled by the overcrowding, finding themselves more productive somewhere else. Thus, they relocated to a, possibly, less well endowed place than the ridge.

This process continued along the time, with the formation of new capital and labor agglomerations that led to the distribution that we observe nowadays, far from the initial prehistorical uniform equilibrium.

1.2 The mathematical model

After the work of Turing [21], it is well-known that the introduction of spatial mechanisms through diffusion terms in reaction systems of differential equations may change the stability of uniform equilibria, developing new ones with a non-uniform pattern-like spatial profile.

Diffusion terms allow to modelize random relocation as well as redistribution due to attractive or repulsive pressures. An important example of reaction-diffusion models is the chemotaxis model, where bacteria relocates due to attraction forces caused by a chemical stimulus. Starting with the works of Keller and Segel [16], this model received much attention due to its modelling simplicity, this is, the easy identification of the biological mechanisms implied by each term of the equations; its analytical and numerical tractability, which allows to prove key mathematical properties like the existence, uniqueness and nonnegativity of the solutions; and, finally, the ability of the model to capture essential aspects of the population behavior [15].

The chemotaxis model is a particular example of reaction cross-diffusion models [19, 3, 9, 6, 7, 10, 13, 14, 1]. The reaction term induces the global growth or decline of the substances, while the diffusive term redistributes them in space. Correspondingly, our model is constructed in two steps.

For the reaction terms, we investigated which functional form should be considered for the labor-capital relationship. As shown in (1)-(2), we found that this relation is mutualistic, in contrast with previous literature, where a predator-pray model was used [4]. We delay the justification of this assumption to Section 2. 2. 2.

For the diffusion part of the system, we assume the following general qualitative laws:

1. (a)

Labor and capital relocate randomly moving from higher to lower density regions. 2. (b)

Labor is repelled by high labor density regions. 3. (c)

Labor and capital are attracted by high capital density regions. However, the attraction felt by capital is limited, due to increasing costs of relocation to such regions.

With these laws operating, the mathematical problem is formulated as follows. Let $L$ and $K$ denote the labor and capital concentrations. Find $L,K:[0,T]\times\bar{\Omega}\subset\mathbb{R}^{n}\to\mathbb{R}_{+}$, with $n\in\mathbb{N}$, such that

[TABLE]

where $Q_{T}=(0,T)\times\Omega$, $\Gamma_{T}=(0,T)\times\partial\Omega$, and the coefficients are non-negative constants.

For the diffusion part of (1)-(2), we adopted the simplest mathematical functional forms for random diffusion (terms containing $c_{1},c_{2}$) and for the intra-repulsion and inter-attraction experienced by labor (terms containing $a_{11}$ and $a_{12}$). However for the intra-attraction experienced by capital (term containing $a_{22}$), we incorporate a saturation function which results on a more complex nonlinear interaction. The reason is the following: high density capital regions are attractive for labor and capital due to work demand and capital synergies, respectively. However, this trend has a limit for the capital due to rising prices of access to real estate and labor, increment of taxes, etc. In this limit, that we capture through the saturation constant $K_{s}$, the attraction function reaches a global maximum at the given geographical region. Thus, injecting additional capital in the place only makes it less attractive. See Remark 3 for further explanations on the saturation term.

In relation to the boundary conditions (3), we assume that labor and capital are isolated within the spatial domain $\Omega$, that is, there are no labor or capital flows through the boundary of $\Omega$. Concerning the initial data (4), we assume that $L_{0}$ and $K_{0}$ are smooth non-negative given functions.

The rest of the paper is organized as follows. In Section 2, we investigate the kind of reaction term which should be adopted in the model. Starting from a general production function, we show that the gradient flow of the corresponding profits function may be approximated by a Lotka-Volterra system. We show that standard economic assumptions imply that this system is of mutualistic type. We illustrate this fact by considering the CES production function, see (9).

In Section 3, we perform a Turing stability analysis to find conditions under which the uniform equilibrium becomes unstable and a new pattern-like equilibrium arises. In geographical economy terms, we check that the spatially uniform coexistence equilibrium, which was prevalent during the long early period of human development, becomes unstable and produces a new non-uniform equilibrium in which geographical patterns arise. These geographical patterns, accumulations of labor and capital, may be interpreted as cities.

Although the linear stability theory is a useful step for understanding pattern formation, it only gives a rough indication of the patterns we should expect. Thus, in the Appendix, we perform a weakly nonlinear analysis based on the method of multiple scales that, in addition to the unstable wave-numbers provided by the linear theory, allows also to approximate the amplitude of such instabilities.

Finally, in Section 4, we illustrate the theoretical results by performing several numerical experiments in which the data lead to the instability of the uniform equilibrium. We approximate the full nonlinear system (1)-(4) with a finite elements scheme and, when possible, compare this approximation with that provided by the weakly nonlinear analysis.

2 The reaction term: parameters of competence and mutualism

A production function, $Y(K,L)$, is an abstraction that summarizes a technical relation between the inputs of labor and capital used in a productive process and the resulting quantity of output.

Firms, in order to operate a production function, must pay the use of productive factors. As a compensation, they obtain income through sales, but at the cost of remunerating the worker hours used and the capital hired. The difference between the revenue sales and the factor costs, constitutes the profits:

[TABLE]

where $p,~{}w,~{}r$ are functions of time describing, respectively, the sale price of one unit of output, the wages paid for a unit of work, and the rental rate of capital used in production. Function $p$ can be taken as numeraire, normalizing it to $1$.

The optimal trajectory $(K(t),L(t))$ is determined by the maximization of profits. The first optimality condition is

[TABLE]

To compute the optimal trajectory, one must prescribe wages and rental rates or, otherwise, state suitable relationships among $L,~{}K,~{}w$ and $r$ that close the system of equations. A simplifying assumption consists on supposing that $w$ and $r$ are close to some stable constant values, say $w^{*},~{}r^{*}$. Then, the profits optimum for these values, $(L_{e}^{*},K_{e}^{*})$, is obtained from the optimality conditions (5).

A dynamics approximating this optimal value is given by the nonlinear problem corresponding to the gradient flow $\partial_{t}(K,L)^{t}=\nabla\Pi(K,L)$. In a further simplification, we may linearize the gradient flow by performing the substitution $\nabla\Pi(K,L)\approx\frac{1}{2}H(\Pi_{e}^{*})(L-L_{e}^{*},K-K_{e}^{*})^{t}$, where $H(\Pi_{e}^{*})$ is the Hessian matrix of $\Pi$ at $(L_{e}^{*},K_{e}^{*})$.

In order to accommodate the gradient flow to the Lotka-Volterra terms, we use a modified gradient flow, see Remark 1, that leads to the system

[TABLE]

with

[TABLE]

In particular, notice that the profits optimum $(L_{e}^{*},K_{e}^{*})$ coincides with the non-trivial equilibrium of the Lotka-Volterra system. In addition, $\tilde{\beta}_{12}=\tilde{\beta}_{21}$ and the second optimality condition for the profits maximization, i.e. that $H(\Pi_{e}^{*})$ is positive definite, imposes the necessary condition $\det(\tilde{B})>0$.

Since, by definition, $H(Y_{e}^{*})=H(\Pi_{e}^{*})$, the signs of $\tilde{\beta}_{ij}$ are determined by the production function. In economic terms, the Law of decreasing marginal productivities states that the more one factor is used, the smaller is the increase in output obtained. This implies that $\partial_{LL}Y$ and $\partial_{KK}Y$ are negative. However, the crossed marginal productivities are increasing, i.e., $\partial_{LK}Y=\partial_{KL}Y>0$. The reason is that if, for instance, labor is abundant, an increase in capital will increase the productivity of workers, and as a consequence it will improve the output, as the farmer with naked hands is much less productive than the farmer with a plow. Therefore, one expects the following signs for the quadratic coefficients of the Lotka-Volterra term:

[TABLE]

capturing a mutualistic model.

Remark 1**.**

Suppose that the dynamical system $\partial_{t}\mathbf{u}=\mathbf{F}(\mathbf{u})$, for some differentiable function $\mathbf{F}:\mathbb{R}^{2}\to\mathbb{R}^{2}$, has a nontrivial asymptotically stable equilibrium $\mathbf{u}^{*}$ with non-negative components. Let $J\mathbf{F}(\mathbf{u}^{*})$ denote its Jacobian matrix at $\mathbf{u}^{*}$, so that $\operatorname{tr}(J\mathbf{F}(\mathbf{u}^{*}))<0$ and $\det(J\mathbf{F}(\mathbf{u}^{*}))>0$.

Consider the system $\partial_{t}\mathbf{w}=\mathbf{G}(\mathbf{w})$, with $\mathbf{G}_{i}(\mathbf{w})=w_{i}\mathbf{F}_{i}(\mathbf{w})$, for $i=1,2$, having the equilibria $\mathbf{0}$ and $\mathbf{u}^{*}$. We have

[TABLE]

and, consequently, $\mathbf{u}^{*}$ is its only stable equilibrium.

2.1 Example: The CES production function

The production function of constant elasticity of substitution (CES) [18] has the following functional form:

[TABLE]

where $A$ is the total factor productivity, a non-negative function of time which is a scaling parameter used to accommodate the measures of $K$ and $L$ in relation to $Y$, and where $\alpha,\beta$ are the income share parameters. We assume the usual parameters range

[TABLE]

For this choice of the production function, we have

[TABLE]

and

[TABLE]

where $Q=\alpha K^{\eta}+\beta L^{\eta}$.

In the following theorem, we show that for a CES production function the signs of the reaction part of the system (6)-(7) are those of a mutualistic Lotka-Volterra model. This gives support to our assumption on the form of the reaction part of the system (1)-(2). We also check that the equilibrium of this model coincides with the optimum of the profits function.

Theorem 1**.**

Let $w^{*},~{}r^{*}$ be positive numbers, and consider the CES production function (9) with parameters satisfying (10). Then, the optimum of the profits function, $\Pi$, is $(L_{e}^{*},K_{e}^{*})$, with

[TABLE]

Moreover, $(K_{e}^{*},L_{e}^{*})$ is also the coexistence equilibrium of the Lotka-Volterra system (6)-(7) with coefficients given by (8), which satisfy

[TABLE]

In addition, if $\epsilon>\eta$ then $\tilde{\beta}_{12}=\tilde{\beta}_{21}>0$.

*Proof. * Dividing the first optimality condition by the second, see (5), and using the explicit expressions (11)-(12) yields

[TABLE]

Replacing in (5) and rearranging terms we obtain the expression of $(L_{e}^{*},K_{e}^{*})$ stated in the proposition. The optimum $(L_{e}^{*},K_{e}^{*})$ is also the coexistence equilibrium of (6)-(7) by construction, see (8). From the definition of the Hessian, see (13), and the condition $\epsilon\geq\eta$ (resp. $\epsilon>\eta$) we easily deduce that $\beta_{12}=\beta_{21}\geq 0$ (resp. $\beta_{12}=\beta_{21}>0$). Moreover, after some computations, we obtain

[TABLE]

Therefore, necessarily $\tilde{\beta}_{11}\tilde{\beta}_{22}>0$. Suppose that both numbers are positive. Then, using their definition, we deduce

[TABLE]

Set $x=\alpha(L_{e}^{*})^{\eta}$, $y=\beta(K_{e}^{*})^{\eta}$ and $\delta=(1-\eta)^{2}/(\epsilon-\eta)^{2}$. The above inequality translates to $\delta(x^{2}+y^{2})+(2\delta-1)xy<0$, which is not possible if $2\delta>1$. And this is indeed the case, due to the conditions $\varepsilon<1$ and $\eta<1$. The contradiction arises after the assumption $\tilde{\beta}_{11}>0$ and $\tilde{\beta}_{22}>0$. Therefore, these quantities must be both negative, since their product is positive. Finally, from the definition of $\tilde{\alpha}_{i}$, see (8), we get $\tilde{\alpha}_{1}=w^{*}(1-\epsilon)>0$, and $\tilde{\alpha}_{2}=r^{*}(1-\epsilon)>0$. $\Box$

3 Turing instability analysis

Before starting the stability analysis, we reduce slightly the profusion of parameters in our model by introducing the assumption that the mutualistic coefficients $\beta_{12}$, $\beta_{21}$ of the reaction terms are positive. This is, for instance, the case for a CES production function with $\epsilon>\eta$, see Theorem 1. Thus, we perform the following change of variables and parameters:

[TABLE]

for some time-space re-scaling factor $\gamma>0$. This change renders problem (1)-(4) to the following form (omitting bars)

[TABLE]

where $(0,T)$ and $\Omega$ are redefined according to the scaling factor $\gamma$. We collect here some assumptions and necessary conditions for the well-posedness of the problem:

The saturation function $g\in C^{2}(\mathbb{R}_{+})$ is non-negative and reaches its maximum at $K_{s}>0$. 2. 2.

The coefficients $c_{i}$, $b$, and $a_{i}$, for $i=1,2$, are non-negative. In particular, we assume the ellipticity conditions

[TABLE] 3. 3.

In agreement with Theorem 1, the coefficients $\alpha_{i}$ and $\beta_{i}$, for $i=1,2$, are non-negative. Moreover, dominance of competitiveness over mutualism, captured by $\beta_{1}\beta_{2}>1$, is also assumed.

Notice that after the change of unknowns leading to the system (14)-(17), the coexistence equilibrium is expressed as

[TABLE]

In the following theorem we establish that diffusion induced instability occurs when the attraction felt by labor toward regions with a high density of capital is large enough.

Theorem 2**.**

The coexistence equilibrium (18), which is stable for the reaction system associated to (14)-(17), becomes unstable for the whole reaction-diffusion system (14)-(17) when $b>b_{c}$, see (26). This instability leads to pattern formation if $\gamma$ is large enough.

*Proof. * Let $(\tilde{L},\tilde{K})$ be a small perturbation of this equilibrium satisfying (14)-(16) and set $(w_{1},w_{2})=(\tilde{L}-L^{*},\tilde{K}-K^{*})$. The linearized problem for $\mathbf{w}=(w_{1},w_{2})$ is given in matrix form by

[TABLE]

where

[TABLE]

Taking into account the boundary conditions, we search for particular solutions of (19)-(21) of the form $\mathbf{w}(t,x)\propto e^{\lambda t+i\mathbf{k}\cdot x}$, where $\lambda$ represents the linear growth rate and $\mathbf{k}$ is the wave number of the perturbation. Introducing this expression of $\mathbf{w}$ in (19)-(21) we are led to the eigenvalue problem

[TABLE]

The corresponding characteristic equation yields the eigenvalues

[TABLE]

so that, for Turing instability to occur, one or more of the eigenvalues $\lambda_{jk}$ must have a positive real part.

For $k=0$ (no diffusion), we have $\operatorname{tr}(A_{0})=\gamma\operatorname{tr}(R)<0$ and $\det(A_{0})=\gamma\det(R)>0$, implying $\operatorname{Re}(\lambda_{j0})<0$ for $j=1,2$. This is, the mode corresponding to the wave number $k=0$ is always stable.

For $k>0$ we have $\operatorname{tr}(A_{k})=\gamma\operatorname{tr}(R)-k^{2}\operatorname{tr}(Q)<0$, and therefore, the only way which may lead to a positive $\operatorname{Re}(\lambda_{jk})$ involves $\det(A_{k})$ being negative. We introduce the convex quadratic polynomial

[TABLE]

with $q=m_{2}-bm_{1}$, and $m_{1},m_{2}$ given by the positive constants

[TABLE]

The minimum of $h$ is attained at

[TABLE]

which is a real root if $q<0$. Hence, a necessary condition for instability is

[TABLE]

The corresponding minimum value is

[TABLE]

which is negative for q<-2\big{(}\det(R)\det(Q)\big{)}^{1/2}, or, in terms of the bifurcation parameter $b$, if $b>b_{c}$, with the critical bifurcation value given by

[TABLE]

Observe that, in particular, $b>b_{c}$ implies the necessary condition (25). The critical wavenumber corresponding to the critical bifurcation values is found by replacing $q^{c}=m_{2}-b_{c}m_{1}$ in (24), yielding

[TABLE]

Thus, for $b>b_{c}$ the system has a range $(k_{1}^{2},k_{2}^{2})$ of unstable wave-numbers, where $k_{1}^{2},k_{2}^{2}$ are the roots of $h$. Pattern formation will arise if the spatial spread of $\Omega$ is large enough so that at least one of the modes admitted by the boundary conditions lies in the interval $(k_{1}^{2},k_{2}^{2})$. Since the roots of $h$ depend linearly on $\gamma$, the time-space scale parameter, pattern formation in the original model will be observed if this parameter is large enough. $\Box$

To gain some insight into the bifurcation condition, we replace the expressions of $m_{1},~{}m_{2},~{}\det(R)$ and $\det(Q)$ in (26) to obtain that $b>b_{c}$ is equivalent to

[TABLE]

where $X^{2}=c_{2}-a_{2}g(K^{*})$ and $Y^{2}=c_{1}+a_{1}L^{*}$. Completing squares, and using the elementary inequality $(x+y)^{2}\leq 2(x^{2}+y^{2})$, we deduce the sufficient condition

[TABLE]

Since $\det(B)=\beta_{1}\beta_{2}-1$, we obtain the following, more strict, sufficient condition, providing a clear interpretation of the bifurcation condition in terms of the diffusion coefficients,

[TABLE]

As expected, we see that the diffusion coefficients promoting stability are those representing random relocation and labor intra-repulsion, while those promoting instability are the intra-attraction coefficients of both labor and capital.

In which respect to the relationship between instability and the kinematic terms, let us assume the simplification $\alpha_{1}=\alpha_{2}=:\alpha$. Replacing these values in the expressions of the equilibrium, and then in (27), we get the sufficient condition

[TABLE]

from where we interpret that if there is a high increase of labor and capital, expressed by a small intra-competitive behavior and a large intrinsic growth, instability will arise.

4 Numerical simulations

We approximate the solutions of the nonlinear problem (14)-(17) by employing two methods: the finite element method (FEM) and the weakly nonlinear approximation (WNL). The latter is explained in detail in the Appendix, but we give here a short summary of the main ideas. See also [11, 12].

Let $\varepsilon$ be a small control parameter representing the dimensionless distance from the critical threshold $\varepsilon^{2}=(b-b_{c})/b_{c}$. The idea of the weakly nonlinear analysis is to expand the perturbation, $\mathbf{w}$, the bifurcation parameter, $b$, and the time scale, $t$, in powers of $\varepsilon$,

[TABLE]

After substitution in the nonlinear system satisfied by $\mathbf{w}$, see (35), and equating with respect to the order of $\varepsilon$, we obtain a chain of linear problems for $\mathbf{w}_{i}$, for $i=1,2,\ldots$, whose solutions are of the form $\mathbf{w}_{i}(T_{1},T_{2},\ldots,x)=A(T_{1},T_{2},\ldots)\mathbf{u}_{i}(x)$, with $A$ denoting amplitude. In particular, the third order approximation leads to the cubic Stuart-Landau equation for the amplitude,

[TABLE]

where the growth rate coefficient $\sigma$ is positive. The dynamics of this equation is divided into two different cases according to the sign of the Landau constant $\ell$: the supercritical case ($\ell>0$), and the subcritical case ($\ell<0$).

In the supercritical case, the equilibrium solution $A_{\infty}=\sqrt{\sigma/\ell}$ is stable, and represents the asymptotic value of the amplitude $A$. The corresponding solution, corrected to satisfy the Neumann boundary conditions, is given by

[TABLE]

where $\bar{k_{c}}$ is the first integer or semi-integer to become unstable when $b$ passes the critical value $b_{c}$.

In the subcritical case, the cubic Stuart-Landau equation does not give any valid information for the amplitude, and a higher degree expansion on powers of $\varepsilon$ must be considered, leading to higher order Stuart-Landau equations.

For the FEM approximation, we used the open source software deal.II [5] to implement a time semi-implicit scheme with a spatial linear-wise finite element discretization. For the time discretization, we take in the experiments a uniform time partition of time step $\tau$. For the spatial discretization, we take a uniform partition of the interval $\Omega=[0,2\pi]$.

Let, initially, $t=t_{0}=0$ and set $(L^{0},K^{0})=(L_{0},K_{0})$. For $n\geq 1$, the discrete problem is: Find $L^{n},K^{n}\in S^{h}$ such that

[TABLE]

for every $\chi\in S^{h}$, the finite element space of piecewise $\mathbb{Q}_{1}$-elements. Here, $(\cdot,\cdot)^{h}$ stands for a discrete semi-inner product on $\mathcal{C}(\overline{\Omega})$.

Since (33)-(34) is a nonlinear algebraic problem, we use a fixed point argument to approximate its solution, $(L^{n},K^{n})$, at each time slice $t=t_{n}$, from the previous approximation $(L^{n-1},K^{n-1})$. Let $L^{n,0}=L^{n-1}$ and $K^{n,0}=K^{n-1}$. Then, for $k\geq 1$ the linear problem to solve is: Find $(L^{n,k},K^{n,k})$ such that for for all $\chi\in S^{h}$

[TABLE]

We use the stopping criteria

[TABLE]

for values of $\text{tol}_{FP}$ chosen empirically, and set $(L^{n},K^{n})=(L^{n,k},K^{n,k})$. Finally, we integrate in time until a numerical stationary solution, $(L^{S},K^{S})$, is achieved. This is determined by

[TABLE]

where $\text{tol}_{S}$ is chosen empirically too.

We run four simulations with data leading to pattern formation. Unless otherwise stated, the following parameters are fixed for all the experiments. The initial data is a perturbation of the stable equilibrium given by $L_{0}=L^{*}(1+0.05\sin(10\pi x))$ and $K_{0}=K^{*}(1+0.05\cos(10\pi x))$. The time step is $\tau=0.01$, and the spatial domain is $\Omega=(0,2\pi$), with 256 spatial nodes. The tolerances for the fixed point method and for the stationary state are taken as $\text{tol}_{FP}=1.e-6$ and $\text{tol}_{S}=1.e-8$. The total factor productivity is $A=1$, the time-scale parameter is $\gamma=1$, the saturation constant is $K_{s}=10K^{*}$, and the bifurcation parameter is $b=1.01b_{c}$. In Table 1, we show the parameters which vary in each experiment and additional information on the numerical results. In Figures 1-3, we show the approximated solutions obtained in the experiments and the determinant of the matrix $A_{k}$, see (23), showing the unstable wave numbers.

Remark 2**.**

For developed economies, the values $\alpha\approx 0.3$ and $\beta\approx 0.6$ are well established. The values of $\epsilon$ and $\beta$ are more difficult to estimate. They can be obtained from the definition of the CES production function (9), from setting the scaling variable $A$ to a known value, and from data indicating that the ratio $K/Y\approx 3.5$ holds in a variety of situations. Of course, the choice of these parameters is not unique. In any case, we found that their role in the stability analysis is limited so, for the examples, we fixed them in order to obtain interesting visualizations.

In the first two experiments, the main reason because instability arises, see (27) and (28), is that the diffusion parameters promoting stability, $c_{1},~{}c_{2}$, and $a_{1}$ are small in comparison with the ratios $\alpha_{i}/\det(B)$, determining the growth capacity of the system. In the Experiment 1, this happens in a symmetric way for the competition coefficients, $\beta_{1}\approx\beta_{2}$ and in an asymmetric way for the intrinsic growth coefficients, with $\alpha_{1}\gg\alpha_{2}$. In the second experiment, the opposite relation is considered. We see that the critical bifurcation parameter is, at least, one order higher than the other diffusion coefficients, implying that, under these situations, instability arises if the attraction felt by labor for capital is high in comparison with random or repulsion effects.

In the third experiment we set larger stability promoting diffusion coefficients and smaller values of the growth and competition coefficients. However, it is clear that the important relationship for instability is that of the size of $\alpha_{i}/\det(B)$ against the size of the diffusion coefficients. Thus, instability emerges as well in this case. However, notice that here the critical bifurcation parameter is of the same order as the other diffusion coefficients. This is, it is not necessary a extreme difference between attraction and repulsion for instability to arise: growth plays a fundamental role. In this experiment, we set $A=100$ and $\gamma=5$ for a better visualization.

In the Experiment 4 we set parameters leading to large amplitude instabilities containing regions abandoned by labor and with few remaining of capital. Unlike the other experiments, in this case we experimented with a subcritical bifurcation so that the third order Stuart-Landau equation for the amplitude does not give any useful information. The stabilizing factor $\beta_{2}a_{1}$, see (27), is large, forcing a large bifurcation coefficient for instabilities to arise. Observe that $\beta_{2}$ is the intra-competition coefficient for $K$, and $a_{1}$ is the intra-repulsion coefficient for $L$. Thus, even if these intra-population stabilizing mechanisms are intense, a large enough attraction of labor for capital is still able to cause severe instabilities. In this experiment, we set $\gamma=10$.

Finally, let us notice the good agreement between the FEM and the NWL approximations in the Experiments 1 to 3, for which the bifurcation is supercritical and the third order Stuart-Landau equation provides an estimate for the steady state amplitude. In all of them, the mean square error between the FEM and the NWL approximations is of the order $\text{MSE}\approx 1.e-2$.

Remark 3**.**

The best place for understanding the contribution of the saturation term to instability is the sufficient condition (27). We have that: (i) the saturation term promotes instability, and (ii) it regulates the size needed for the bifurcation parameter to overpass the critical threshold.

In the experiments 1 and 4 this term is small in comparison to the other diffusion terms. In the experiment 2, we have $a_{2}\sim c_{1},c_{2}$ and $a_{2}\gg a_{1}$. Finally, in the experiment 3 we have $a_{2}\gg a_{1},c_{1},c_{2}$ . However, instability arises in all these experiments. Therefore, it does not seem that the saturation term plays any special role for the emergence of instabilities. However, this term is crucial for the existence theory. Indeed, if this term were not uniformly bounded then the well-posedness of the problem could not be guaranteed due to the lack of parabolicity.

Appendix: Weakly nonlinear analysis in 1D

Let $\mathbf{w}$ be the random perturbation around the coexistence equilibrium introduced in Section 3. We recast the nonlinear system (14)-(15) as

[TABLE]

where $\mathcal{L}^{b}=\gamma R+Q^{b}\partial_{xx}$, and $\mathcal{N}^{b}$ is a nonlinear operator. Here, we introduced in the notation the superscript $b$ to stress the dependence of these operators on the bifurcation parameter. We consider the decomposition

[TABLE]

with the bilinear symmetric forms

[TABLE]

and the nonlinear operator

[TABLE]

Introducing the parameter $\varepsilon^{2}=(b-b_{c})/b_{c}$, and the expansions (29)-(31), replacing them in the system (35), and equating in terms of the order of $\varepsilon$, leads to the following systems of equations, that we set in $(0,T)\times(0,2\pi/k_{c})$ and complement with non-flow boundary conditions:

[TABLE]

where $\mathcal{L}^{b_{c}}=\gamma R+Q^{b_{c}}\partial_{xx}$, and the elements of the matrix $B_{j}$ are zero, with the exception of $(B_{j})_{12}=-b_{j}L^{*}$. We compute the solutions of (36)-(38).

**Order $\varepsilon$: **The solution of the linear problem (36) is given by \mathbf{w}_{1}=A(T_{1},T_{2})\mbox{\boldmath\rho}\cos(k_{c}x), with \mbox{\boldmath\rho}\in\ker(\gamma R-k_{c}^{2}Q^{b_{c}}), where $A$ is the amplitude of the pattern, unknown at the moment. Observe that $\gamma R-k_{c}^{2}Q^{b_{c}}=A_{k_{c}}^{b_{c}}$, with $A_{k}$ defined in (23), and with $k_{c}$ determined to yield $\det(A_{k_{c}}^{b_{c}})=0$. Therefore, the vector $\rho$ is defined up to a multiplicative constant, that we shall fix later, see (41).

**Order $\varepsilon^{2}$: ** We have B_{1}\partial_{xx}\mathbf{w}_{1}=-Ak_{c}^{2}\cos(k_{c}x)B_{1}\mbox{\boldmath\rho}. On noting that \mathcal{Q}_{U}(\mathbf{w}_{1},\mathbf{w}_{1})=A^{2}\mathcal{Q}_{U}(\mbox{\boldmath\rho},\mbox{\boldmath\rho})\cos^{2}(k_{c}x), for $U=R,~{}Q$, and using standard trigonometric identities, we find that

[TABLE]

with $\mathcal{M}_{j}=\mathcal{Q}_{R}-j^{2}k_{c}^{2}\mathcal{Q}_{Q}.$ Computing the other terms of (37) yields the problem

[TABLE]

By Fredholm’s alternative, (37) admits a solution if and only if \langle\mathbf{F},\mbox{\boldmath\psi}\rangle_{L^{2}}=0, where $\langle\cdot,\cdot\rangle_{L^{2}}$ denotes the scalar product in $L^{2}(0,2\pi/k_{c})$, and \mbox{\boldmath\psi}\in\ker((\mathcal{L}^{b_{c}})^{*}) is of the form

[TABLE]

For similar reasons than $\rho$, $\eta$ is defined up to a multiplicative constant. We fix $\eta$ in (41), where it is also shown that \langle\mbox{\boldmath\rho},\mbox{\boldmath\eta}\rangle\neq 0.

The compatibility condition implies that the terms in $\varepsilon$ of the expansions of $b$ and $T$ are secular terms, and thus we impose $T_{1}\equiv 0$ and $b_{1}\equiv 0$, implying $A\equiv A(T_{2})$. With these restrictions, the Fredholm’s alternative is satisfied, and we look for a solution of (37) of the form $\mathbf{w}_{2}=A^{2}\sum_{j=0,2}\mathbf{w}_{2j}\cos(jk_{c}x)$, for which

[TABLE]

Then, $\mathcal{L}^{b_{c}}\mathbf{w}_{2}=\mathbf{F}$ if the vectors $\mathbf{w}_{2j}$ are the solutions of the linear systems

[TABLE]

**Order $\varepsilon^{3}$: ** Replacing $\mathbf{w}_{1},\mathbf{w}_{2}$ in the terms of $\mathbf{G}$, see (38), yields

[TABLE]

Using that $\mathcal{Q}_{R}$ and $\mathcal{Q}_{Q}$ are bilinear and recalling the definition of $\mathcal{M}$, we get

[TABLE]

Computing the other nonlinear terms and replacing them in the definition of $\mathbf{G}$ given in (38), we obtain

[TABLE]

The solvability condition for problem (38) is \langle\mathbf{G},\mbox{\boldmath\psi}\rangle_{L^{2}}=0, with \mbox{\boldmath\psi}=\mbox{\boldmath\eta}\cos(k_{c}x) given by (39). This condition leads to the differential equation

[TABLE]

where \mathbf{G}_{1}=k_{c}^{2}B_{2}\mbox{\boldmath\rho}=-k_{c}^{2}L^{*}b_{2}\rho_{2}\mathbf{e}_{1}, and

[TABLE]

Thus, we deduce the cubic Stuart-Landau equation for the amplitude

[TABLE]

with \sigma=-\langle\mathbf{G}_{1},\mbox{\boldmath\eta}\rangle/\langle\mbox{\boldmath\rho},\mbox{\boldmath\eta}\rangle and \ell=\langle\mathbf{G}_{3},\mbox{\boldmath\eta}\rangle/\langle\mbox{\boldmath\rho},\mbox{\boldmath\eta}\rangle. We, finally, fix the vectors \mbox{\boldmath\rho}\in\ker(\gamma R-k_{c}^{2}Q^{b_{c}}) , and \mbox{\boldmath\eta}\in\ker((\gamma R-k_{c}^{2}Q^{b_{c}})^{*}). Since $\gamma R_{21}-k_{c}^{2}Q_{21}^{b_{c}}=\gamma K^{*}>0$, see (22), we choose the forms \mbox{\boldmath\rho}=(M,1)^{t}, and \mbox{\boldmath\eta}=(1,M^{*})^{t}, where

[TABLE]

With this election, we have \langle\mathbf{G}_{1},\mbox{\boldmath\eta}\rangle=-k_{c}^{2}b_{2}L^{*}<0, and \langle\mbox{\boldmath\rho},\mbox{\boldmath\eta}\rangle>0, and therefore, the growth rate coefficient $\sigma$ is always positive. As mentioned at the beginning of Section 4, the dynamics of (40) in the supercritical case, $\ell>0$, has the stable equilibrium solution $A_{\infty}=\sqrt{\sigma/\ell}$, representing the asymptotic value of the amplitude. The corresponding solution may be approximated by (32).