The Stability of Steady-State Hot-Spot Patterns for a Reaction-Diffusion Model of Urban Crime

The existence and stability of localized patterns of criminal activity are studied for the reaction-diffusion model of urban crime that was introduced by Short et. al. [Math. Models. Meth. Appl. Sci., 18, Suppl. (2008), pp. 1249--1267]. Such patterns, characterized by the concentration of criminal activity in localized spatial regions, are referred to as hot-spot patterns and they occur in a parameter regime far from the Turing point associated with the bifurcation of spatially uniform solutions. Singular perturbation techniques are used to construct steady-state hot-spot patterns in one and two-dimensional spatial domains, and new types of nonlocal eigenvalue problems are derived that determine the stability of these hot-spot patterns to ${\mathcal O}(1)$ time-scale instabilities. From an analysis of these nonlocal eigenvalue problems, a critical threshold $K_c$ is determined such that a pattern consisting of $K$ hot-spots is unstable to a competition instability if $K>K_c$. This instability, due to a positive real eigenvalue, triggers the collapse of some of the hot-spots in the pattern. Furthermore, in contrast to the well-known stability results for spike patterns of the Gierer-Meinhardt reaction-diffusion model, it is shown for the crime model that there is only a relatively narrow parameter range where oscillatory instabilities in the hot-spot amplitudes occur. Such an instability, due to a Hopf bifurcation, is studied explicitly for a single hot-spot in the shadow system limit, for which the diffusivity of criminals is asymptotically large. Finally, the parameter regime where localized hot-spots occur is compared with the parameter regime, studied in previous works, where Turing instabilities from a spatially uniform steady-state occur.


Introduction
Recently, Short et. al. [29,30,31] introduced an agent-based model of urban crime that takes into account repeat or near-repeat victimization. In dimensionless form, the continuum limit of this agent-based model is the two-component reaction-diffusion PDE system x ∈ Ω ; ∂ n A = 0 , x ∈ ∂Ω , (1.1 a) where the positive constants ε 2 , D, α, γ and τ , are all assumed to be spatially independent. In this model, P (x, t) represents the density of the criminals, A(x, t) represents the "attractiveness" of the environment to burglary or other criminal activity, and the chemotactic drift term −2D∇ · P ∇A A represents the tendency of criminals to move towards sites with a higher attractiveness. In addition, α is the baseline attractiveness, while (γ − α)/τ represents the constant rate of re-introduction of criminals after a burglary. For further details on the model see [29].
In [29], the reaction-diffusion system (1.1) with chemotactic drift term was derived from a continuum limit of a lattice-based model. It was then analyzed using linear stability theory to determine a parameter range for the Turing instability leads to a formation of three hot-spots; one is annihilated almost immediately due to a fast-time instability, while the second hot-spot is annihilated after a long time. (b) D = 0.5 with all other parameters as in (a). Two hot-spots remain stable. (c) Numerical solution of (1.1) in a two-dimensional square of width 4. Parameters are α = 1, γ = 2, ε = 0.08, τ = 1, D = 1. Initial conditions are P (x, 0) = 1 − α/γ and A(x, 0) = γ(1 + rand * 0.001) where rand generates a random number between 0 and 1.
existence of a Turing instability of the spatially uniform steady-state. A weakly nonlinear theory, based on a multiscale expansion valid near the Turing bifurcation point, was developed in [30,31] for (1.1) for both one and twodimensional domains. This theoretical framework is very useful to explore the origins of various patterns that are observed in full numerical solutions of the model. However, the major drawback of a weakly nonlinear theory is that the parameters must be tuned near the bifurcation point of the Turing instability. When the parameters values are at an O(1) distance from the bifurcation point, an instability of the spatially homogeneous steady-state often leads to patterns consisting of localized structures. Such localized patterns for the crime model (1.1), consisting of the concentration of criminal activity in localized spatial regions, are referred to as either hot-spot or spike-type patterns. A localized hot-spot solution, not amenable to an analytical description by a weakly nonlinear analysis, was observed in the full numerical solutions of [30].
As an illustration of localization behavior, in Fig. 1(a) we plot the numerical solution to (1.1) in the one-dimensional domain Ω = [0, 1] with parameter values α = 1, γ = 2, ε = 0.02, τ = 1, and D = 1. The initial conditions, consisting of a small mode-three perturbation of the spatially homogeneous steady-state A e = γ and P e = (γ − α)/γ are first amplified due to linear instability. Shortly thereafter, nonlinear effects become significant and the solution quickly becomes localized leading to the formation of three hot-spots, as shown at t ≈ 14. Subsequently, one of the hot-spots appears to be unstable and is quickly annihilated. The remaining two hot-spots drift towards each other over a long time, until finally around t ≈ 180, 000, another hot-spot is annihilated. The lone remaining hot-spot then drifts towards the center of the domain where it then remains. Next, in Fig. 1(b) we re-run the simulation when D is decreased to D = 0.5 with all other parameters the same as in Fig. 1(a). For this value of D, we observe that the final state consists of two hot-spots. Similar complex dynamics of hot-spots in a two-dimensional domain are shown in Fig. 1(c).
It is the goal of this paper to give a detailed study of the existence and stability of steady-state localized hot-spot patterns for (1.1) in both one and two-dimensional domains in the singularly perturbed limit (1. 2) The assumption that ε 2 ≪ D implies that the length-scale associated with the change in the attractiveness of potential burglary sites is much smaller than the length-scale over which criminals explore new territory to commit crime. In this limit, a singular perturbation methodology will be used to construct steady-state hot-spot solutions and to derive new nonlocal eigenvalue problems (NLEP's) governing the stability of these solutions. From an analysis of the spectrum of these NLEP's, explicit stability thresholds in terms of D and τ for the initiation of O(1) time-scale instabilities of these patterns are obtained. In a one-dimensional domain, an additional stability threshold on D for the initiation of slow translational instabilities of the hot-spot pattern is derived. Among other results, we will be able to explain both the fast and slow instabilities of the localized hot-spots patterns as observed in Fig. 1(a).
In related contexts, there is now a rather large literature on the stability of spike-type patterns in two-component reaction-diffusion systems with no drift terms. The theory was first developed in a one-dimensional domain to analyze the stability of steady-state spike patterns for the Gierer-Meinhardt model (cf. [10,3,35,37,38,34,43,46]) and, in a parallel development, the Gray-Scott model (cf. [4,5,15,23,24,18,1]). The stability theory for these two models was extended to two-dimensional domains in [39,41,40,44,43,2]. Related studies for the Schnakenburg model are given in [11,36,42]. The dynamics of quasi-equilibrium spike patterns is studied for one-dimensional domains in [9,6,7,33,22], and in a multi-dimensional context in [13,14,16,2]. More recently, in [19] the stability of spikes was analyzed for a reaction-diffusion model of species segregation with cross-diffusion. A common feature in all of these studies, is that an analysis of the spectrum of various classes of NLEP's is central for determining the stability properties of localized patterns. A survey of NLEP theory is given in [46], and in, a broader context, a survey of phenomena and results for far-from-equilibrium patterns is given in [25].
In contrast, for reaction-diffusion systems with chemotactic drift terms, such as the crime model (1.1), there are only a few studies of the existence and stability of spike solutions. These previous studies have focused mainly on variants of the well-known Keller-Segel model (cf. [8,12,28,32]).
We now summarize and illustrate our main results. In §2.1 we construct a multi hot-spot steady state solution to (1.1) on a one-dimensional interval of length S. We refer to a symmetric hot-spot steady-state solution as one for which the hot-spots are equally spaced and, correspondingly, each hot-spot has the same amplitude. In §2.2 asymmetric steady-state hot-spot solutions, characterized by unevenly spaced hot-spots, are shown to bifurcate from the symmetric branch of hot-spot solutions at a critical value of D.
In §3 we study the stability of steady-state K-hot-spot solutions on an interval of length S when τ = O(1). A singular perturbation approach is used to derive a NLEP that determines the stability of these hot-spot patterns to O(1) time-scale instabilities. In contrast to the NLEP's arising in the study of spike stability for the Gierer-Meinhardt model (cf. [35]), this NLEP is explicitly solvable. In this way, a critical threshold K c+ is determined such that a pattern consisting of K hot-spots with K > 1 is unstable to a competition instability if and only if K > K c+ . This instability, which develops on an O(1) time scale as ε → 0, is due to a positive real eigenvalue, and it triggers the collapse of some of the hot-spots in the pattern. This critical threshold K c+ > 0 is the unique root of (see Principal Result 3.2 below) In addition, from the location of the bifurcation point associated with the birth of an asymmetric hot-spot equilibrium, a further threshold K c− is derived that predicts that a K-hot-spot steady-state with K > 1 is stable with respect to slow translational instabilities of the hot-spot locations if and only if K < K c− . This threshold is given explicitly by (see (3.20) below) Since K c− < K c+ , the stability properties of a K-hot-spot steady-state solution with K > 1 and τ = O(1) are as follows: stability when K < K c− ; stability with respect to O(1) time-scale instabilities but unstable with respect to slow translation instabilities when K c− < K < K c+ ; a fast O(1) time-scale instability dominates when K > K c+ .
As an illustration of these results consider again Fig. 1 Fig. 1(a). We then predict from the threshold K c− that the two-hot-spot solution will become unstable on a very long time interval. This is also confirmed by the full numerical solutions shown in Fig. 1(a). In contrast, if we decrease D to D = 0.5 as in Fig. 1(b) then we calculate from (1.3) and (1.4) that K c+ ≈ 2.612 and K c− ≈ 2.372. Our prediction is that the three hot-spot solution that emerges from initial data will be unstable on an O(1) time-scale, but that a two-hot-spot steady-state will be stable. These predictions are again corroborated by the full numerical results.
In §4 we examine oscillatory instabilities of the amplitudes of the hot-spots in terms of the bifurcation parameter τ in (1.1). From an analysis of a new NLEP with two separate nonlocal terms, we show that an oscillatory instability of the hot-spot amplitudes as a result of a Hopf bifurcation is not possible on the regime τ ≤ O(ε −1 ). This non-existence result for a Hopf bifurcation is in contrast to the results obtained in [35] for the Gierer-Meinhardt model showing the existence of oscillatory instabilities of the spike amplitudes in a rather wide parameter regime. However, for the asymptotically larger range of τ with τ = O(ε −2 ), in §4.1 we study oscillatory instabilities of a single hot-spot in the simplified system corresponding to letting D → ∞ in (1.1). In this shadow system limit, we show for a domain of length one that low frequency oscillations of the spot amplitude due to a Hopf bifurcation will occur when τ > τ c where In §5 we extend our results to two dimensional domains. We first construct a quasi-equilibrium multi hot-spot pattern, and then derive an NLEP governing O(1) time-scale instabilities of the spot pattern. As in the analyses of [39,40,41,42,43,44] for the Gierer-Meinhardt and Gray-Scott models, our existence and stability theory for localized hot-spot solutions is accurate only to leading-order in powers of −1/ log ε. In §5.1, we show from an analysis of a certain NLEP problem that for τ = O(1) a K-hot-spot solution with K > 1 is unstable when K > K c where as ε → 0. Here w is the radially symmetric ground-state solution of ∆w − w + w 3 = 0 in R 2 and |Ω| is the area of Ω. As an example, consider the parameter values as in Fig. 1(c), for which |Ω| = 16. Then, from (1.5) we get K c ≈ 44.48. Starting with random initial conditions, we observe from Fig. 1(c) that at t = 9000 we have K = 7.5 hot-spots, where we count boundary spots having weight 1/2 and corner spots having weight 1/4. Since K < K c , this is in agreement with the stability theory.
Finally, in §5, we contrast results for Turing instabilities and Turing patterns with our results for localized hot-spots.
We also propose a few open problems.
2 Asymptotic Analysis of Steady-State Hot-Spot Solutions in 1-D so that (2.1) transforms to To motivate the ε-dependent re-scaling of V that facilitates the analysis below, we suppose that D ≫ l 2 and we integrate the steady-state of (2 . In addition, from the steady-state of (2.3 a), we conclude that in the inner region near a hot-spot centered at . This implies that −3p + (3p − 1) = −p, so that p = 1. Therefore, for D ≫ l 2 we conclude that V = O(ε 2 ) globally on −l < x < l, while A = O(ε −1) in the inner region near a hot-spot. Finally, in the outer region we must have A = O(1), so that from the steady-state of (2.3 b), we conclude that , this balance requires that D = O(ε −2 ). Since V = O(ε 2 ) globally, while A = O(ε −1 ) in the core of a hot-spot, we conclude that within a hot-spot of criminal activity the density P = V A 2 of criminals is O(1).
In summary, this simple scaling analysis motivates the introduction of new O(1) variables v and D 0 defined by In terms of (2.4), (2.3) transforms to

A Single Steady-State Hot-Spot Solution
We will now construct a steady-state hot-spot solution on the interval −l < x < l with a peak at the origin. In order to construct a K-hot-spot pattern on a domain of length S, with evenly spaced spots, we need only set l = S/(2K) and perform a periodic extension of the results obtained below on the basic interval −l < x < l. As such, the fundamental problem considered below is to asymptotically construct a one-hot-spot steady-state solution on −l < x < l.
In the inner region, near the center of the hot-spot at x = 0, we expand A and v as From (2.5 a) we obtain, in terms of y, that A j (y) for j = 0, 1 satisfy In contrast, from (2.5 b), we obtain that v j for j = 0, 1 satisfy In order to match to an outer solution, we require that v 0 and v 1 are bounded as |y| → ∞. In this way, we then obtain that v 0 and v 1 must both be constants, independent of y.
We look for a solution to (2.7) for which the hot-spot has a maximum at y = 0. The homoclinic solution to (2.7 a) with A ′ 0 (0) = 0 is written as where w is the unique solution to the ground-state problem given explicitly by w = √ 2 sech y. Next, we decompose the solution A 1 to (2.7 b) as where w 1 (y) satisfies with w ′ 1 (0) = 0 and w 1 → 0 as |y| → ∞. A key property of the operator L 0 , which relies on the cubic exponent in (2.10), is the remarkable identity that (2.12) The proof of this identity is a straightforward manipulation of (2.10) and the operator L 0 in (2.11). This property plays an important role in an explicit analysis of the spectral problem in §3. Here this identity is used to provide an explicit solution to (2.11) in the form In this way, in the inner region the two-term expansion for A in terms of the unknown constants v 0 and v 1 is In the outer region, defined for ε ≪ |x| ≤ l, we have that v = O(1) and that A = O(1). From (2.5), we obtain that subject to the matching condition that h 0 → v 0 as x → 0 ± . The solution to this problem gives the outer expansion (2.14) Next, we must calculate the constants v 0 and v 1 appearing in (2.13) and (2.14). We integrate (2.5 b) over −l < x < l and use v x = 0 at x = ±l to get In (2.15), we emphasize that the first two terms on the left-hand side arise solely from the inner expansion, whereas the O(ε 2 ) term would be obtained from both the inner and outer expansions. By equating coefficients of ε in (2.15), Then, upon using (2.13) for A 0 and A 1 , together with w = √ 2 sech y, (2.16) We summarize our result for a single steady-state hot-spot solution as follows: Principal Result 2.1: Let ε → 0, and consider a one-hot-spot solution centered at the origin for (2.5) on the interval |x| ≤ l. Then, in the inner region y = x/ε = O(1), we have In addition, in the inner region, the leading-order steady-state criminal density P from (2.1) is P ∼ w 2 . Here w =   Table 1. Comparison of numerical and asymptotic results for the amplitude A max ≡ A(0) and for v(0) of a one-hotspot solution on [−1, 1] with D 0 = 1, γ = 1, and α = 2. The 1-term and 2-term asymptotic results for A max and v(0) are obtained from (2.17).
w(y) = √ 2 sech y is the homoclinic of (2.10), while v 0 and v 1 are given in (2.16). In the outer region, O(ε) < |x| ≤ l, Note that to get a solution for A which is uniformly valid in both inner and outer region, we can combine the formulas (2.17) and (2.18). The resulting first-order composite solution is given explicitly by For a specific parameter set, a comparison of the full numerical steady-state solution of (2.5) with the composite asymptotic solution (2.19) is shown in Fig. 2. A comparison of numerical and asymptotic values for A(0) and v(0) at various ε is shown in Table 1. From this table we note that the two-term asymptotic expansion for v(0) agrees very favorably with full numerical results.

Asymmetric Steady-State K-Hot-Spot Solutions
In the limit ε → 0, we now construct an asymmetric steady-state K-hot-spot solution to (2.5) in the form of a sequence of hot-spots of different heights. This construction will be used to characterize the stability of symmetric steady-state K-hot-spot solutions with respect to the small eigenvalues λ = o(1) in the spectrum of the linearization. Since the asymmetric solution is shown to bifurcate from the symmetric branch, the point of the bifurcation corresponds to a zero eigenvalue crossing along the symmetric branch. To determine this bifurcation point, we compute v(l) for the one-hot-spot steady-state solution to (2.5) on |x| ≤ l, where l > 0 is a parameter. This canonical problem is shown to have two different solutions. A K-hot-spot asymmetric solution to (2.5) is then obtained by using translates of these two local solutions in such a way to ensure that the resulting solution is C 1 continuous. Since the details of the construction of the asymmetric solution is very similar to that in [36] for the Schnakenburg model, we will only give a brief outline of the analysis.
The key quantity of interest is the critical value D s 0K of D for which an asymmetric K-hot-spot solution branch bifurcates off of the symmetric branch. To this end, we first calculate from (2.14) that where the function B(z) on 0 < z < ∞ is defined by The function B(z) > 0 in (2.20 b) has a unique global minimum point at z = z c = 1, and it satisfies B ′ (z) < 0 on [0, z c ) and B ′ (z) > 0 on (z c , ∞). Therefore, given any z ∈ (0, z c ), there exists a unique pointz ∈ (z c , ∞) such that . This shows that given any l, with l/q < z c = 1, there exists a uniquel, withl/q ≡z > z c = 1, such We refer to solutions of length l andl as A-type and B-type hot-spots. Now consider the interval To construct a K-hot-spot steady-state solution to (2.5) on this interval with K 1 ≥ 0 hot-spots of type A and K 2 = K − K 1 ≥ 0 hot-spots of type B, arranged in any order across the interval, we must solve the At this value of the parameters, a steady-state K-hot-spot asymmetric solution branch bifurcates off of the symmetric K-hot-spot branch. This critical value of D 0 determines the small eigenvalue stability threshold in the linearization of the symmetric K-hot-spot steady-state solution. For a symmetric configuration of K hot-spots on an interval of length S we have 2Kl = S so that the critical value D 0 = D S 0K , as defined by (2.21), can be written as A more detailed construction of the asymmetric solution branches parallels that done in [36] for the Schnakenburg model and is left to the reader.

The NLEP Stability of Steady-State 1-D Hot-Spot Patterns
We now study the stability of the K-hot-spot steady-state solution to (2.5) that was constructed in §2. The analysis for the "large" O(1) eigenvalues in the spectrum of the linearization is done in several distinct steps.
Firstly, we let A e , v e denote the one-hot-spot quasi-steady-state solution to (2.5) on the basic interval |x| ≤ l, which was given in Principal Result 2.1. Upon introducing the perturbation we obtain from the linearization of (2.5) that We consider (3.2 a) and (3.2 b) on |x| ≤ l subject to the Floquet-type boundary conditions where z is a complex parameter.
For simplicity, in this section we will set τ = 0 in (3.2 b). The analysis of the possibility of Hopf bifurcations induced by taking τ = 0 is studied in §4.
After formulating the NLEP associated with solving (3.2) for arbitrary z, we then must determine z so that we and hence to obtain periodic boundary conditions on an interval of length 2Kl we require that z K = 1, so that By using these values of z j in the NLEP problem associated with (3.2), we obtain the stability threshold of a Khot-spot solution on a domain of length 2Kl subject to periodic boundary conditions. The last step in the analysis is then to extract the stability thresholds for the corresponding Neumann problem from the thresholds for the periodic problem, and to choose l appropriately so that the Neumann problem is posed on [−1, 1]. This is done below. This Floquet-based approach to determine the NLEP problem of a K-hot-spot steady-state solution for the Neumann problem has been used previously for reaction-diffusion systems exhibiting mesa patterns [22], for the Gierer-Meinhardt model [34], and for a cross-diffusion system [19].
We now implement the details of this calculation. The asymptotic analysis for ε → 0 of (3.2) proceeds as follows.
In the inner region with y = ε −1 x, we use A e = O(ε −1 ) and v ex ≪ 1, to obtain from (3.2 b) that to leading order w 2 ψ y y = 0 in the inner region, where w is the homoclinic satisfying (2.10). To prevent exponential growth for ψ as |y| → ∞, we must take ψ = ψ 0 where ψ 0 is a constant to be determined. Then, for (3.2 a) we look for a localized inner eigenfunction in the form Upon using the leading-order approximation , we obtain to leading order that Φ(y) satisfies Here v 0 is given in (2.16).
In the outer region, away from the hot-spot centered at yields the leading-order outer problem subject to the Floquet-type boundary conditions (3.2 c). The matching condition for the inner and outer represen- However, the problem for ψ is not yet complete, as it must be supplemented by appropriate jump conditions for ψ x across x = 0.
We now proceed to derive this jump condition. We first define an intermediate scale η satisfying ε ≪ η ≪ 1, and We use the limiting behavior as x → 0 ± of the outer expansion to calculate the terms on the left hand-side of (3.7).
From A e ∼ α, (2.14) to calculate v ex (0 ± ), and (3.5) to calculate φ(0 ± ), we obtain that Next, since η ≫ O(ε), we can estimate the integrals on the right-hand side of (3.7) by their contributions from the Upon substituting (3.8) and (3.9) into (3.7), we obtain the following jump condition for ψ x across x = 0: For the range λ > −1, we can neglect the negligible O(ε 2 ) term in the jump condition (3.10). In this way, the problem for the outer eigenfunction ψ(x) is to solve subject to the continuity condition ψ(0 + ) = ψ(0 − ) = ψ 0 and the following jump condition across x = 0: Upon calculating ψ(0) = ψ 0 from this problem, the NLEP is then obtained from (3.4).
Upon solving (3.11) for ψ(x), and evaluating the result at Upon substituting these results into (3.4), we obtain the following NLEP for a K-hot-spot steady-state on a domain of length 2Kl subject to periodic boundary conditions: 14 a) The final step in the analysis is extract the NLEP for the Neumann problem from the NLEP (3.14) for the periodic problem. More specifically, the stability thresholds for a K-hot-spot pattern with Neumann boundary conditions can Therefore, to obtain the NLEP problem governing the stability of an steady-state K-hot-spot pattern on an interval of length S subject to Neumann boundary conditions, we simply replace cos(2πj/K) with cos(πj/K) in (3.14) and then set l = S/(2K) in the NLEP of (3.14). In this way, we obtain the following main result: Principal Result 3.1: Consider a K-hot-spot solution to (2.5) on an interval of length S subject to Neumann boundary conditions. For ε → 0, and τ = O(1), the stability of this solution with respect to the "large" eigenvalues λ = O(1) of the linearization is determined by the spectrum of the NLEP where w(y) is the homoclinic solution satisfying w ′′ − w + w 3 = 0.

Lemma 3.2:
Consider the NLEP problem

16)
for an arbitrary constant c corresponding to eigenfunctions for which Then, on this range there is only one element in the point spectrum, and it is given explicitly by To prove this we consider only the region Re(λ) > −1, where we can guarantee that |Φ| → 0 exponentially as |y| → ∞. The continuous spectrum for (3.16) is λ < −1, with λ real. To establish (3.17) we use Green's identity on w 2 and Φ, which is written as In obtaining the last equality in (3.18) we calculated the integrals using w = √ 2 sech y. Since χ 0 = 3 < χ 1 < χ 2 < . . . χ K−1 , a one-hot-spot solution is stable for all D 0 , while the instability threshold for a multi hot-spot pattern is set by χ K−1 . In this way, we obtain the following main stability result.
Principal Result 3.2: Consider a K-hot-spot solution to (2.5) on an interval of length S with K > 1 subject to Neumann boundary conditions. For τ = 0, and in the limit ε → 0, this solution is stable on an O(1) time-scale In terms of the original diffusivity D, given by D = ε −2 D 0 , the stability threshold is D L K = ε −2 D L 0K when K > 1. Alternatively, a one-hot-spot solution is stable for all D 0 > 0, provided that D 0 is independent of ε.
Although we have not calculated the stability threshold for the small eigenvalues for which λ → 0 as ε → 0 in the spectrum of the linearization (3.2), we conjecture that this stability threshold is the same critical value D S 0K of D 0 , given in (2.22), for which an asymmetric K-hot-spot steady-state branch bifurcates off of the symmetric K-hot-spot branch. This simple approach to calculate the small eigenvalue stability threshold, which avoids the lengthy matrix manipulations of [10], has been validated for the Gierer-Meinhardt, Gray-Scott, and Schnakenburg reaction-diffusion models in [37], [36], and [18]. Since D S 0K < D L 0K for K ≥ 2, we conclude that a symmetric K-hot-spot steady-state solution is stable with respect to both the large and the small eigenvalues only when D < D S 0K . We make two remarks. Firstly, for the case of a single hot-spot where K = 1 we expect that the stability threshold for D 0 will be exponentially large in 1/ε, and similar to that derived in [14] for the Gierer-Meinhardt model in the near-shadow limit. Secondly, we remark that the possibility of stabilizing multiple hot-spots for (2.1) is in direct contrast to the result obtained in the analysis of [32] of spike solutions for a Keller-Segel-type chemotaxis model with a logarithmic sensitivity function for the drift term. For this chemotaxis problem of [32], only a one-spike solution can be stable.

Numerical Results
We now compare our stability predictions with results from full numerical solutions of (2.3). As derived above, under Neumann boundary conditions the thresholds on D for the stability of a symmetric K-hot-spot pattern on a domain of length 2KL are .

Hopf Bifurcation of K-Hot-Spot Steady-State Solutions
In this section we study the spectrum of (3.2) for τ > 0. This is done by first deriving an NLEP similar to (3.15).
Since the analysis leading to the new NLEP is very similar to that in §3, we will only outline it here briefly.
For τ ≪ O(ε 2 ), we get ψ ≈ ψ 0 in the inner region x = O(ε), and hence (3.4) for the inner approximation Φ(y) for φ remains valid. For τ ≪ O(ε 2 ), we get to leading-order that ψ xx = 0 in the outer region 0 < |x| ≤ l and so (3.6) still holds. However, for τ = 0, the jump conditions (3.7)-(3.9) must be modified. In place of (3.7), we get The left hand-side of (4.1) was estimated in (3.8), while the first two terms on the right-hand side of (4.1) were estimated in (3.9). We then use and v e ∼ v 0 , to estimate the last term on the right hand-side of (4.1) as Upon substituting (3.8), (3.9), and (4.2), into (4.1), we obtain that With this modification of the coefficients in (3.11 b), the outer problem for ψ is still (3.11).

Principal Result 4.2:
There is no value of τ c > 0 for which the NLEP of (4.7) has a Hopf bifurcation.
We note that there is a key step in the derivation of Principal   Derivation of Principal Result 4.2: We use the notation h dy ≡ ∞ ∞ h dy. Upon using Green's identity w 2 LΦ dy = ΦLw 2 dy and L 0 w 2 = 3w 2 , together with (4.7 a), we obtain 3 − 3χ j w 5 dy w 3 dy w 2 Φ dy − χ 1j 2 w 5 dy wΦ dy = λ w 2 Φ dy .
Upon solving for w 2 Φ dy in terms of wΦ dy, and then substituting into (4.7 a), we get (4.8 a) We then simplify f (λ) by using (4.7 b), together with w 5 = (3/2) w 3 , to obtain .
Next, we observe that the eigenvalues λ of the NLEP problem (4.8) are the roots of the transcendental equation Upon recalling that L 0 w = 2w 3 , we calculate Substituting this result together with (4.8 b) and w 2 = 4 into (4.9), we obtain that λ is a root of . (4.10) To determine whether a Hopf bifurcation is possible we set λ = iλ I in (4.10) and replace (L 0 − λ) −1 by (L 0 − λ) −1 = L 2 0 + λ 2 I −1 (L 0 + iλ I ). Then, upon comparing the real and imaginary parts in the resulting expression, we obtain that µ ≡ λ 2 I and any Hopf bifurcation threshold τ c must be the roots of the coupled system Upon eliminating τ c from (4.11), we obtain a transcendental equation solely for µ = λ 2 I > 0, By using the identities L 0 w = 2w 3 and L −1 0 w = (w + yw ′ ) /2, the limiting behavior for ρ(µ) is readily calculated as Moreover, direct numerical computations of ρ(µ) show that it is an increasing function of µ (see Fig. 4). On the other hand, for µ > 0 we have 3 − 2β j − 2βj 9 µ < 1 since β j ≥ 1. It follows that (4.12) cannot have any solution with µ > 0. Consequently, there is no Hopf bifurcation on the parameter regime τ = O(ε −1 ).
Such a non-existence result for Hopf bifurcations for the crime model when τ = O(ε −1 ) is qualitatively very different than for the Gierer-Meinhardt and Gray-Scott models, analyzed in [35] and [15], where Hopf bifurcations occur in wide parameter regimes.

A Hopf Bifurcation for the Shadow Limit
Principal Result 4.2 has shown that there is no Hopf bifurcation for the regime τ = O(ε −1 ). However, a Hopf bifurcation can and does appear when τ = O(ε −2 ). As will be shown below, in such a regime the amplitude of the hot-spot becomes oscillatory with an asymptotically large temporal period, due to an eigenvalue that is dominated, to leading-order in ε, by its pure imaginary part. To illustrate this phenomenon, in this section we analytically derive the condition for a Hopf bifurcation of a single boundary spot on a domain of length one. To further simplify our computations, we will assume that D 0 in (2.5 b) is taken sufficiently large such that v(x, t) = v(t) can be approximated by a time-dependent constant. The limit D → ∞ is called the shadow-limit (cf. [38]). Our main result is the following: Principal Result 4.3: Suppose that D 0 ≫ 1 and consider a half hot-spot of (2.5) located at the origin on the domain Let τ = τ 0 /ε 2 . Then, there is a Hopf bifurcation at τ 0 = τ 0c . That is, the hot spot is stable for τ 0 < τ 0c and is unstable for τ 0 > τ 0c . Destabilization takes place via a Hopf bifurcation. More precisely, when τ 0 = O(1), the related stability problem has an eigenvalue near the origin with the following asymptotic behavior as ε → 0: (4.14) Numerical example. To illustrate Principal Result 4.3 we take γ = 4, α = 1, ε = 0.05 and D 0 = 1000. Then, (4.13) yields τ 0c ≈ 1.07378. Now take τ 0 = 0.95 so that (4.14) yields the eigenvalue λ ≈ 0.4082i − 0.00529. We then expect the single hot-spot to be stable, although it will exhibit long transient oscillations. From the eigenvalue, we can estimate the period of the oscillation to be P = 2π 0.4082 ≈ 15.39. This agrees with full numerical solutions of (4.15) as shown in Fig. 5(a).
Next, we increase τ 0 to 1.15, while keeping the other parameters the same. In this case, τ 0 > τ 0c so that the hot-spot is unstable in the limit ε → 0. However, τ 0 = 1.15 is very close to the threshold value, and with ε = 0.05, we expect even longer transients with the final state still unclear at t = 300. This behavior is shown in Fig. 5(b).
Finally, as shown in Fig. 5(c), when we increase τ 0 to 1.35, we clearly observe oscillations of an increasing amplitude. we then obtain to leading order that v(x, t) ∼ v(x). We then integrate (2.5 b) and use the no-flux boundary conditions to obtain the following shadow-limit system on x ∈ [0, 1]: where we have defined µ by The shadow problem (4.15) is the starting point of our analysis. The corresponding steady-state system is As shown below, in order to analyze the Hopf bifurcation it is necessary to construct the steady-state solution to two orders in ε. To do so, we let y = x/ε and we expand where the operator L 0 is defined by L 0 u ≡ u yy − u + 3w 2 u. This operator has several key readily derivable identities, which allows us to determine the solution u 1 to (4.18) as Next, to determine v 0 and v 1 we must calculate the integral in (4.16) for ε ≪ 1. This yields that By equating coefficients of ε, we get that v 0 and v 1 satisfy Upon using the solution (4.19) for u 1 , the unknown v 1 can be determined in terms of a quadrature as The integrals defining v 0 and v 1 are then calculated explicitly by using which yields the explicit formulae Next, we study the stability of this solution. For convenience, we extend the problem to the interval [−1, 1] by even reflection. We linearize (4.15) around the steady-state solution to obtain the eigenvalue problem where the constant ψ denotes the perturbation in v. Upon solving for ψ we obtain (4.21) To motivate the analysis below, we first suppose that τ λε ≫ 1. Then, by using u ∼ w/ √ v 0 and v ∼ v 0 , (4.21) reduces to leading order to the NLEP Here and below, f denotes ∞ −∞ f dy. This problem has a zero eigenvalue corresponding to the eigenfunction φ = w. All other discrete eigenvalues satisfy Re(λ) < 0 (cf. [46]). Therefore, the critical eigenvalue will be a perturbation of the zero eigenvalue.
A posteriori computations shows that the correct anzatz is in fact The analysis below shows that λ 0 is purely imaginary, and hence determines the frequency of the oscillation, but not its stability. Therefore, a two-term expansion in λ must be obtained in order to determine the stability of the oscillations. As such, we must expand all quantities in the shadow problem up to O(ε). The delicate part in the calculation is to note that where the last integral is in fact O(ε) as a result of Thus, this term "jumps" an order and is comparable in magnitude to ε 2u 0 u 1 . The remaining part of the analysis is more straightforward. We let τ = τ 0 ε −2 and expand ψ in (4.21) as so that the eigenvalue problem for φ from (4.21) becomes (4.22) After tedious but straightforward computations, the three operators in (4.22) are given by in terms of the coefficients c 0 , . . . , c 9 defined by ; Finally, we expand φ in (4.22) as and equate powers of ε in (4.22). This yields the following problems for φ 1 and φ 2 : To determine λ 0 and φ 1 we must formulate the appropriate solvability condition based on the adjoint operator L * 1 of L 1 defined by Since L 1 admits a zero eigenvalue of multiplicity one, then so does L * 1 . In fact, w * defined by w * ≡ (yw y + w) /2 , (4.28) is the unique element in the kernel of L * 1 , i.e. L * 1 w * = 0, owing to the following two readily derived identities: Next, we impose a solvability condition on (4.25) in the usual way. We multiply (4.25) by w * and integrate by parts to derive that where we used the integral identity in (4.29) together with w * w = w 2 /4. From the formulae for the coefficients c 0 and c 1 in (4.24), (4.30) determines λ 0 as Since λ 0 is purely imaginary, the next order term λ 1 needs to be computed to determine stability.
The problem (4.25) for φ 1 can be written by using (4.23 a) and (4.23 b) as Since L 0 w * = w and L 0 w = 2w 3 , we can write φ 1 as Upon substituting φ 1 into the definition of d, we can then solve for d by using (4.30) for λ 0 to get Finally, by using w * w = w 2 /4, the expression above simplifies to d = c 0 w 2 + c 1 w 3 so that φ 1 is given explicitly from (4.32) as (4.33) With φ 1 explicitly known, we impose the solvability condition on the problem (4.26) for φ 2 to determine λ 2 as Upon using (4.33) for φ 1 , L 3 w from (4.23 c), and (L 2 − λ 0 ) from (4.23 b), the integrals above are evaluated as Finally, upon substituting c 0 , . . . , c 9 from (4.24) into this expression, we determine λ 1 explicitly as The two-term expansion for λ given in (4.14) follows from (4.31) and (4.34). The Hopf bifurcation threshold is obtained by setting λ 1 = 0. This occurs precisely as τ 0 is increased past τ 0c , where τ c is given by (4.13).

Hot-Spot Patterns in 2-D: Equilibria and Stability
In this section we construct a K-spot quasi-steady-state solution to (1.1) in an arbitrary 2-D domain with spots centered at x 1 , . . . , x K . To leading-order in σ = −1/ log ε, we then derive a threshold condition on the diffusivity D for the stability of the K-spot quasi-steady-state solution to instabilities that develop on an O(1) time-scale.
As in the analysis of hot-spot patterns in one spatial dimension, we set V = P/A 2 (see (2.2)) into (1.1) to obtain We first motivate the ε-dependent re-scalings of V and A that are needed for the 2-D case. We suppose that D ≫ 1, so that V is approximately constant. By integrating the steady-state equation of (5.1 b) over Ω we get . In addition, from the steady-state of (5.1 a), we must have in globally on Ω, while A = O(ε −2 ) in the inner region near a hot-spot. Finally, in the outer region we must have . Finally, in the core of a hot-spot we conclude that the density P of criminals, given by Although this simple scaling analysis correctly identifies the algebraic factors in ε, there are more subtle logarithmic terms of the form σ ≡ −1/ log ε that are needed in the construction of the quasi-steady-state hot spot solution.
The scaling analysis above motivates the introduction of new variables v, u, and D defined by In terms of (5.2), (5.1) transforms exactly to Owing to the non-uniformity in the behavior as |y| → ∞ of the solution to this core problem near a spot (see below), the construction of a quasi-steady-state K-spot solution for (5.3) is more intricate than that for the Gierer-Meinhardt, Schnakenburg, or Gray-Scott problems analyzed in [39]- [44], [13], [16], and [1]. As such, we will only develop a theory that is accurate to leading order in σ ≡ −1/ log ε, similar to that undertaken in [39]- [44], and [13]. This is in contrast to the recent approach in [16] and [1] that used a hybrid asymptotic-numerical method to construct quasi-steady-state spot patterns to the Schnakenburg and Gray-Scott systems, respectively, with an error that is beyond-all-orders with respect to σ.
In the outer region, away from the spots centered at x 1 , . . . , x K , we expand u and v as where σ = −1/ log ε and D = D 0 /σ, where D 0 = O(1). From the steady-state of (5.3 b) we obtain that h 0 is constant, and that h 1 satisfies As shown below, this problem must be augmented by certain singularity conditions that are obtained by matching the outer solution for v to certain inner solutions, one in the neighborhood of each spot.
In the inner region near the j-th spot centered at x j we introduce the inner variables y, u j , and v j by In terms of these inner variables, and with D = σ −1 D 0 , the steady-state of (5.3) transforms exactly on y ∈ R 2 to where σ ≡ −1/ log ε. We will construct a radially symmetric solution u j = u j (ρ), v j = v j (ρ) to this problem, where The complication in analyzing (5.7) is that u j = O(1) for |y| = O(1), whereas u j = O(ε 2 ) for |y| ≫ 1. Therefore, the "diffusivity" u 2 j in the operator for v in (5.7 b) ranges from O(1) when |y| = O(1) to O(ε 4 ) when |y| ≫ 1. For this reason, we cannot simply neglect the second term on the left-hand side of (5.7 b) for all |y|. However, as we show below, we can neglect the O(ε 6 ) term on the right-hand side of (5.7 b).
Next, we calculate the far-field behavior, valid for |y| ≫ 1, for the solution v j to (5.7 b). To do so, we define the ball B δ = {y | |y| ≤ δ}, where 1 ≪ δ ≪ O(ε −1 ). Therefore, this ball is defined in the intermediate matching region between the inner and outer scales y and x, respectively. Upon integrating (5.7 b) over B δ , and using the divergence theorem, we obtain that , the integral on the right hand-side of (5.11) can be estimated by using u j ∼ v −1/2 j0 w. In contrast, on the left hand-side of (5.11) we use u j ∼ αε 2 on ρ = δ ≫ 1. In this way, we obtain that (5.11) becomes Since δ ≪ O(ε −1 ), we can neglect the last term on the right hand-side of (5.12), which is equivalent to neglecting the right-hand side of (5.7 b) for the inner problem.
From (5.12) we obtain that the far-field behavior for |y| ≫ 1 for the solution to (5.7 where S j is defined by Therefore, the appropriate core problem determining the asymptotic shape of the hot-spot profile is to seek a radially symmetric solution to 14) The next step in the construction of the multi hot-spot quasi-steady-state pattern is to match the inner and outer solutions for v in order to determine v j0 . We let y = ε −1 (x − x j ) in (5.13 a) to obtain that the outer solution for v must have the singularity behavior Upon comparing (5.15) with the outer expansion v ∼ h 0 + σh 1 + · · · from (5.4), we conclude that

a)
and that h 1 satisfies (5.5) subject to the singularity behaviors Upon using the divergence theorem, the problem for h 1 has a solution only when the solvability condition is satisfied, where |Ω| is the area of Ω. When this condition is satisfied, the solution for h 1 can be written as whereh 1 is a constant to be determined, and where G(x; x i ) is the Neumann Green's function satisfying In summary, the asymptotic matching provides the following algebraic system for determining v j0 for j = 1, . . . , K: A symmetric K-hot-spot quasi-steady-state solution corresponds to a solution of (5.19) for which v j0 = v 0 for all j.
From (5.9), this solution is characterized by the fact that the hot-spot profile is, to leading-order, the same for each j. The result for such symmetric quasi-equilibria is summarized as follows: regime D = ε −4 D 0 /σ with σ = −1/ log ε, is characterized as follows: In the inner region near the j-th hot-spot, where 20) where b = ∞ 0 w 3 ρ dρ, and w(ρ) with ρ = |y| is the ground-state solution of (5.10). Alternatively, in the outer region where |x − x j | ≫ O(ε) for j = 1, . . . , K, then Here h 1 is given in (5.17) in terms of the Neumann Green's function, and h 0 is a constant given by For a symmetric hot-spot pattern, the source strength S j in (5.17) is the same for each j, and is given in terms of To illustrate Principal Result 5.1, we let Ω be the unit disk centered at the origin, and we choose α = 1, γ = 2, ε = 0.05, and D = 500ε −2 . We consider a single hot-spot at the center of the disk. Upon numerically computing the ground-state solution w satisfying (5.10) we obtain that w(0) ≈ 2.2062 , Alternatively, from the full numerical solution of the radially symmetric steady-state solution of (5.3), we compute that v(0) ≈ 23.017 and u(0) ≈ 0.455. The error is about 0.5% for v(0) and about 0.8% for u(0). A comparison of asymptotic and numerical results is shown in Fig. 6.
As a remark, the general shape of the function F (v) defined in (5.19) also shows that there can be asymmetric K-spot quasi-equilibria corresponding to spots of two distinct heights. Similar asymmetric patterns in 2-D have been constructed for the Gierer-Meihnardt and Gray-Scott systems in [43]) and ( [44]. Let K s and K b be non-negative integers denoting the number of small and large spots, respectively, with K = K s + K b . Then, from (5.19), a K-spot asymmetric pattern is constructed by determining two distinct values v 0s and v 0b satisfying where v min = 3 (c/2) 3/2 and c is defined in (5.19). The spatial profile of the hot-spot of large and small amplitude is w, respectively.
We will not investigate the solvability with respect to D 0 of the algebraic system (5.23) governing asymmetric spot patterns. Instead, in the next subsection we will study the stability properties of the symmetric K-hot-spot quasi-steady-state solution given in Principal Result 5.1.

The Stability of Hot-Spot Quasi-Steady-State Patterns
We linearize (5.3) around the quasi-steady-state K-hot-spot pattern to obtain the eigenvalue problem In our stability analysis we will consider the range of D where D = ε −4 D 0 /σ and D 0 = O(1). It is on this parameter range of D that a stability threshold occurs. Our main stability result is as follows: Principal Result 5.2: Consider the K-spot quasi-steady-state solution of (5.3) as constructed in Principal Result 5.1 for ε ≪ 1. Assume that τ = τ 0 /ε 2 where τ 0 = O(1). Then, for ε → 0, and to leading order in σ = −1/ log ε, the stability on an O(1) time-scale of a K-hot-spot quasi-steady-state solution with K ≥ 2 is determined by the spectrum of the two distinct NLEP's with Φ → 0 as |y| → ∞. Here κ i for i = 1, 2 are defined by in terms of the constants β and D 0c defined by For a one-hot-spot solution, for which K = 1, there is only a single NLEP with κ i in (5.25 a) replaced by κ 1 . For and is unstable when D 0 > D 0c . In terms of the unscaled D, this yields the stability threshold For K = 1 and τ 0 ≪ 1, a single hot-spot is stable for all D 0 independent of ε.
To show that D 0c is the stability threshold of D 0 for K ≥ 2, we let τ 0 → 0 in (5.25) to obtain the limiting NLEP's where κ can assume either of the two values For (5.27), the rigorous result in Theorem 3 of [45] establishes the existence of an eigenvalue with Re(λ) > 0 whenever κ < 2. Thus, since κ 20 < 2 when D 0 > D 0c , we conclude that a K-hot-spot quasi-equilibria with K > 2 is unstable when D 0 > D 0c . In addition, the result in Theorem 1 of [45] (see the remark following Theorem 1) proves that all eigenvalues satisfy Re(λ) < 0 when 2 < κ ≤ 6. Since κ 10 = 3 and 2 < κ 20 < 3 for all D 0 < D 0c , it follows that the NLEP for a K-hot-spot quasi-steady-state solution with K ≥ 2 does not have unstable eigenvalues when D 0 < D 0c .
For K = 1, which corresponds to a one-hot-spot solution, the only choice for κ is κ = κ 10 = 3, and so we have stability for any D 0 > 0 independent of ε.
We make two remarks. Firstly, we anticipate that a single hot-spot solution where K = 1 will be stable provided that D 0 is not exponentially large in 1/ε. The analysis to determine this stability threshold in the near-shadow limit should be similar to that done for the Gierer-Meinhardt model in [14]. Secondly, it is an open question to analyze is much more difficult than in 1-D since the identity L 0 w 2 = 3w 2 for the local operator L 0 Φ ≡ ∆ y Φ − Φ + 3w 2 Φ in R 1 no longer holds in R 2 . Moreover, since this identity does not hold in the 2-D case, we cannot determine λ explicitly when τ 0 ≪ 1 as was done for the 1-D case in Lemma 3.2.
The stability threshold in Principal Result 5.2 follows once we derive the NLEP (5.25). The derivation of this NLEP is done in several distinct steps.
solution for ψ as We then obtain that ψ 0 is a constant, given by and that ψ 1 satisfies The solvability condition for (5.40) is that K j=1 B j = 0. Upon summing (5.39) from j = 1, . . . , K, we obtain that where the constants c i for i = 1, . . . , 4 are defined by
By solving (5.42) forΨ 0 , and substituting the resulting expression into (5.30), we obtain the vector NLEP where the matrix M is defined by The matrix M −1 is a rank-one update of a scalar multiple of the identity matrix. As such, its spectrum Mω = κω can readily be calculated as , ω T j e = 0 , j = 2, . . . , K .

Discussion
We have studied localized hot-spot solutions of (1.1) in one and two spatial dimensions in the regime ε 2 ≪ 1 with D ≫ 1. In this large D limit, steady-state multi hot-spot solutions have been constructed and their stability properties investigated with respect to D and τ from the analysis of certain nonlocal eigenvalue problems. An open problem is to characterize the dynamics of multi hot-spot patterns by reducing (1.1) to a finite dimensional dynamical system for the locations of the hot-spots in a quasi-steady-state pattern as was done for various two-component reaction-diffusion models without drift terms in [6,7,9,33,16,1].
We now remark on how the localized states constructed in this paper are related to the weakly nonlinear Turing patterns studied in [29,30,31]. In contrast to the theory developed in [30,31], our parameter values for (1.1) are not restricted to lie close to the Turing bifurcation point. In the limit ε → 0, the spatially homogeneous steady-state solution for (1.1) is A e = γ and P e = (γ − α)/γ. For ε → 0, it is linearly unstable when (see equation (2.8) of [30]) γ > 3 2 α , as ε → 0 , and a spatially heterogeneous solution bifurcates off the homogeneous steady state at γ ∼ 3 2 α as ε → 0. In contrast, a localized hot-spot exists for the wider parameter range γ > α. In Fig. 7 we plot the numerically computed bifurcation diagram for A(0) vs. γ on a one-dimensional interval, with other parameters as indicated in the figure caption. The Turing bifurcation at γ = 3 2 α is subcritical when ε ≪ 1 (cf. [30]), which is consistent with the stability of the constant state when γ < 3 2 α. However, the bifurcation curve quickly turns around as it enters the localized regime. This is consistent with the existence of localized states when γ > α.

(5.2)
Although l turing and l c are not related, they are both of the same asymptotic order O(D 1/4 ε 1/2 ), which implies that the number of stable localized hot-spots corresponds roughly to the most unstable Turing mode. In fact, for a large parameter range, the inequality l c < l turing holds. For example, when τ = 1, α = 1, γ = 2 we calculate that l c /l turing = 0.7 < 1. This was already observed empirically in Fig. 6 of [29].
Formula (1.4) provides an upper bound K c on the number K of stable hot-spots in the regime where D ≫ 1.
On the other hand, numerical evidence shows the existence of a lower bound that occurs when D is sufficiently small. In this regime an instability occurs when there are too few hot-spots. In fact, if the hot-spot inter-distance exceeds some critical length, then a hot-spot insertion phenomena is observed (see Fig. 8). From Fig. 8, it appears that hot-spot insertion takes place every time that D is quartered. A similar insertion phenomenon occurs in other reaction-diffusion systems such as a chemotaxis model [26], the Gray-Scott model [27], the ferrocyanide-iodide-sulfite  system [20], the Brusselator [17], the Schnakenburg model [16], and a model of droplet breakup [21]. The detailed analysis of hot-spot insertion phenomena in (1.1) will be considered in future work.