# The field-road diffusion model: fundamental solution and asymptotic behavior

Matthieu Alfaro, Romain Ducasse, Samuel Tréton

► To cite this version:

Matthieu Alfaro, Romain Ducasse, Samuel Tréton. The field-road diffusion model: fundamental solution and asymptotic behavior. Journal of Differential Equations, 2023, 10.1016/j.jde.2023.05.002. hal-03712723

HAL Id: hal-03712723

https://hal.science/hal-03712723

Submitted on 4 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# The field-road diffusion model: fundamental solution and asymptotic behavior

Matthieu Alfaro, Romain Ducasse, Samuel Tréton

<figure><img src="image_1.png" /></figure>

## To cite this version:

Matthieu Alfaro, Romain Ducasse, Samuel Tréton. The field-road diffusion model: fundamental solution and asymptotic behavior. Journal of Differential Equations, 2023, 10.1016/j.jde.2023.05.002. hal-03712723v2

HAL Id: hal-03712723

https://hal.science/hal-03712723v2

Submitted on 3 May 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The field-road diffusion model: fundamental solution and asymptotic behavior

Matthieu ALFARO¹, Romain DUCASSE², and Samuel TRÉTON¹

¹Université de Rouen Normandie, CNRS, Laboratoire de Mathématiques Raphaël Salem (LMRS), Saint-Étienne-du-Rouvray, France.

²Université Paris Cité and Sorbonne Université, CNRS, Laboratoire Jacques-Louis Lions (LJLL), F-75006 Paris, France.

## Abstract

We consider the linear field-road system, a model for fast diffusion channels in population dynamics and ecology. This system takes the form of a system of PDEs set on domains of different dimensions, with exchange boundary conditions. Despite the intricate geometry of the problem, we provide an explicit expression for its fundamental solution and for the solution to the associated Cauchy problem. The main tool is a Fourier (on the road variable)/Laplace (on time) transform. In addition, we derive estimates for the decay rate of the $L^\infty$ norm of these solutions.

**Key Words:** field-road model, diffusion, exchange boundary conditions, fundamental solution, decay rate.

**AMS Subject Classifications:** 35K40, 35C15, 92D25, 35B40.

## 1 Introduction

In this work we consider the solutions $(v = v(t, x, y), u = u(t, x))$ to the linear field-road model

$$ egin{cases} 
\partial_t v = d\Delta v, & t > 0, x \in \mathbb{R}^{N-1}, y > 0, \\ 
-d\partial_y v|_{y=0} = \mu u - \nu v|_{y=0}, & t > 0, x \in \mathbb{R}^{N-1}, \\ 
\partial_t u = D\Delta u + \nu v|_{y=0} - \mu u, & t > 0, x \in \mathbb{R}^{N-1}. 
\end{cases} 
\quad (1.1) $$

Here $N \ge 2$ and $d, D, \mu$ and $\nu$ are positive constants. This system is actually the linear (purely diffusive) part of the field-road model, introduced by Berestycki, Roquejoffre and Rossi [10] to describe the spread of invasive species in presence of networks with fast propagation. We review the original system later in this section.

Let us mention that phenomena of spatial spread are highly relevant to the understanding of biological invasions, spreads of emergent diseases, as well as spatial shifts in distributions in the context of climate change. There is a wide literature dedicated to these topics, let us refer for instance to [21, 22] and the reference therein for more details. In the recent years, there has been a growing recognition of the importance of fast diffusion channels on biological propagations: for instance, an accidental transportation via human activities
# The field-road diffusion model: fundamental solution and asymptotic behavior

Matthieu ALFARO¹, Romain DUCASSE², and Samuel TRÉTON¹

¹Université de Rouen Normandie, CNRS, Laboratoire de Mathématiques Raphaël Salem (LMRS), Saint-Étienne-du-Rouvray, France.

²Université Paris Cité and Sorbonne Université, CNRS, Laboratoire Jacques-Louis Lions (LJLL), F-75006 Paris, France.

## Abstract

We consider the linear field-road system, a model for fast diffusion channels in population dynamics and ecology. This system takes the form of a system of PDEs set on domains of different dimensions, with exchange boundary conditions. Despite the intricate geometry of the problem, we provide an explicit expression for its fundamental solution and for the solution to the associated Cauchy problem. The main tool is a Fourier (on the road variable)/Laplace (on time) transform. In addition, we derive estimates for the decay rate of the $L^\infty$ norm of these solutions.

**Key Words:** field-road model, diffusion, exchange boundary conditions, fundamental solution, decay rate.

**AMS Subject Classifications:** 35K40, 35C15, 92D25, 35B40.

## 1 Introduction

In this work we consider the solutions $(v = v(t, x, y), u = u(t, x))$ to the linear field-road model

$$ egin{cases} 
\partial_t v = d\Delta v, & t > 0, x \in \mathbb{R}^{N-1}, y > 0, \\ 
-d\partial_y v|_{y=0} = \mu u - \nu v|_{y=0}, & t > 0, x \in \mathbb{R}^{N-1}, \\ 
\partial_t u = D\Delta u + \nu v|_{y=0} - \mu u, & t > 0, x \in \mathbb{R}^{N-1}. 
\end{cases} 
\quad (1.1) $$

Here $N \ge 2$ and $D, \mu$ and $\nu$ are positive constants. This system is actually the linear (purely diffusive) part of the field-road model, introduced by Berestycki, Roquejoffre and Rossi [10] to describe the spread of invasive species in presence of networks with fast propagation. We review the original system later in this section.

Let us mention that phenomena of spatial spread are highly relevant to the understanding of biological invasions, spreads of emergent diseases, as well as spatial shifts in distributions in the context of climate change. There is a wide literature dedicated to these topics, let us refer for instance to [21, 22] and the reference therein for more details. In the recent years, there has been a growing recognition of the importance of fast diffusion channels on biological propagations: for instance, an accidental transportation via human activities
of some individuals towards northern and eastern France may be the cause of accelerated propagation of the pine processionary moth [26]. Also, in Canada, some GPS data revealed that wolves travel faster along seismic lines (i.e. narrow strips cleared for energy exploration), thus increasing their chances to meet a prey [20]. It is also acknowledged that fast diffusion channels (roads, airlines, etc.) play a central role in the propagation of epidemics. As is well known, the spread of the black plague, which killed about a third of the European population in the 14th century, was favoured by the trade routes, especially the Silk Road, see [28]. More recently, some evidences of the radiation of the COVID epidemic along highways and transportation infrastructures were found [17].

In order to capture the phenomena induced by such fast diffusion channels, the so-called field-road reaction diffusion system

$$ \left\{ \begin{array}{ll} \partial_t v = d\Delta v + f(v), & t > 0, x \in \mathbb{R}^{N-1}, y > 0, \\ -d\partial_y v|_{y=0} = \mu u - \nu v|_{y=0}, & t > 0, x \in \mathbb{R}^{N-1}, \\ \partial_t u = D\Delta u + \nu v|_{y=0} - \mu u, & t > 0, x \in \mathbb{R}^{N-1}, \end{array} \right. \quad (1.2) $$

was proposed by Berestycki, Roquejoffre and Rossi [10]. The mathematical problem then amounts to describing survival and propagation in a non-standard physical space: the geographical domain consists in the half-space (the “field”) $x \in \mathbb{R}^{N-1}$, $y > 0$, bordered by the hyperplane (the “road”) $x \in \mathbb{R}^{N-1}$, $y = 0$. In the field, individuals diffuse with coefficient $d > 0$ and their density is given by $v = v(t, x, y)$. In particular $\Delta v$ has to be understood as $\Delta_x v + \partial_{yy} v$. On the road, individuals typically diffuse faster ($D > d$) and their density is given by $u = u(t, x)$. In particular $\Delta u$ has to be understood as $\Delta_x u$. The exchanges of population between the road and the field are described by the second equation in system (1.2), where $\mu > 0$ and $\nu > 0$. These Robin type boundary conditions link the field and the road equations and, in some sense, are the core of the model.

In a series of works [10, 9, 12, 11], Berestycki, Roquejoffre and Rossi studied the field-road system with $N = 2$ and $f$ a Fisher-KPP nonlinearity. They shed light on an acceleration phenomenon: when $D > 2d$, the road enhances the global diffusion and the spreading speed exceeds the standard Fisher-KPP invasion speed. Since then, many generalizations or related problems have been studied. Field-road models with (periodic) heterogeneities are considered in [18], [33], [1]. Situations where the field is not a half-space but a more general domain are studied in [31], [27], [15], [13]. The setting where the Laplace operators are replaced with non-local operators, accounting for more complex diffusion processes, is tackled in [6, 5]. The papers [7, 8] consider the interaction between fast-diffusion networks and ecological phenomena such as climate change. In [23, 25, 24], the author introduces long range exchanges so that the road may catch individuals from anywhere in the field and vice versa.

Despite these results on (1.2) and its variations, it turned out that the linear field-road system — obtained by letting $f \equiv 0$ in (1.2) — was not fully understood: there were no complete expressions for the fundamental solution, and the decay rate of the $L^\infty$ norm of the solutions to (1.1) as $t \to +\infty$ was not known.

In the present paper, we thus consider (1.1) as a starting point and intend at filling this gap. It turns out that we can reach an explicit expression for both the fundamental solution and the solution to the associated Cauchy problem. To do so, we perform a “Fourier in $x$/Laplace in $t$” transform, which seems to be the good way to understand the intricate exchange boundary conditions. These explicit expressions enable, in particular, to provide a sharp (possibly up to a logarithmic term) decay rate of the $L^\infty$ norm of the solution.
of some individuals towards northern and eastern France may be the cause of accelerated propagation of the pine processionary moth [26]. Also, in Canada, some GPS data revealed that wolves travel faster along seismic lines (i.e. narrow strips cleared for energy exploration), thus increasing their chances to meet a prey [20]. It is also acknowledged that fast diffusion channels (roads, airlines, etc.) play a central role in the propagation of epidemics. As is well known, the spread of the black plague, which killed about a third of the European population in the 14th century, was favoured by the trade routes, especially the Silk Road, see [28]. More recently, some evidences of the radiation of the COVID epidemic along highways and transportation infrastructures were found [17].

In order to capture the phenomena induced by such fast diffusion channels, the so-called field-road reaction diffusion system

$$ \left\{ \begin{array}{ll} \partial_t v = d\Delta v + f(v), & t > 0, x \in \mathbb{R}^{N-1}, y > 0, \\ -d\partial_y v|_{y=0} = \mu u - \nu v|_{y=0}, & t > 0, x \in \mathbb{R}^{N-1}, \\ \partial_t u = D\Delta u + \nu v|_{y=0} - \mu u, & t > 0, x \in \mathbb{R}^{N-1}, \end{array} \right. \quad (1.2) $$

was proposed by Berestycki, Roquejoffre and Rossi [10]. The mathematical problem then amounts to describing survival and propagation in a non-standard physical space: the geographical domain consists in the half-space (the “field”) $x \in \mathbb{R}^{N-1}$, $y > 0$, bordered by the hyperplane (the “road”) $x \in \mathbb{R}^{N-1}$, $y = 0$. In the field, individuals diffuse with coefficient $d > 0$ and their density is given by $v = v(t, x, y)$. In particular $\Delta v$ has to be understood as $\Delta_x v + \partial_{yy} v$. On the road, individuals typically diffuse faster ($D > d$) and their density is given by $u = u(t, x)$. In particular $\Delta u$ has to be understood as $\Delta_x u$. The exchanges of population between the road and the field are described by the second equation in system (1.2), where $\mu > 0$ and $\nu > 0$. These Robin type boundary conditions link the field and the road equations and, in some sense, are the core of the model.

In a series of works [10, 9, 12, 11], Berestycki, Roquejoffre and Rossi studied the field-road system with $N = 2$ and $f$ a Fisher-KPP nonlinearity. They shed light on an acceleration phenomenon: when $D > 2d$, the road enhances the global diffusion and the spreading speed exceeds the standard Fisher-KPP invasion speed. Since then, many generalizations or related problems have been studied. Field-road models with (periodic) heterogeneities are considered in [18, 33, 1]. Situations where the field is not a half-space but a more general domain are studied in [31, 27, 15, 13]. The setting where the Laplace operators are replaced with non-local operators, accounting for more complex diffusion processes, is tackled in [6, 5]. The papers [7, 8] consider the interaction between fast-diffusion networks and ecological phenomena such as climate change. In [23, 25, 24], the author introduces long range exchanges so that the road may catch individuals from anywhere in the field and vice versa.

Despite these results on (1.2) and its variations, it turned out that the linear field-road system — obtained by letting $f \equiv 0$ in (1.2) — was not fully understood: there were no complete expressions for the fundamental solution, and the decay rate of the $L^\infty$ norm of the solutions to (1.1) as $t \to +\infty$ was not known.

In the present paper, we thus consider (1.1) as a starting point and intend at filling this gap. It turns out that we can reach an explicit expression for both the fundamental solution and the solution to the associated Cauchy problem. To do so, we perform a “Fourier in $x$/Laplace in $t$” transform, which seems to be the good way to understand the intricate exchange boundary conditions. These explicit expressions enable, in particular, to provide a sharp (possibly up to a logarithmic term) decay rate of the $L^\infty$ norm of the solution.
Such a linear rate estimate is a key stone information to tackle difficult nonlinear issues. For instance, it is expected, see [16] or [3], that such a decay rate determines the so-called Fujita exponent $p_F$ for the system (1.2) when $f(v) = v^{1+p}$: this exponent separates “systematic blow-up” (when $0 < p \le p_F$) from “possible extinction” (when $p > p_F$). Also, understanding the Fujita blow-up phenomena when $f(v) = v^{1+p}$ typically enables, see [2], to solve the issue of “extinction vs propagation” in the population dynamics model (1.2) with $f(v) = v^{1+p}(1 - v)$. Such a nonlinearity serves a model for the so-called Allee effect [4], roughly meaning that the per capita growth rate of the population is not maximal at small density. We plan to address such issues in a future work.

## 2 Main results

Let $N \ge 2$ be an integer. A generic point $X \in \mathbb{R}^N$ will be written $X = (x, y)$ with $x \in \mathbb{R}^{N-1}$ and $y \in \mathbb{R}$. For the “variable of integration”, we reserve the notation $Z = (z, \omega) \in \mathbb{R}^N$, with $z \in \mathbb{R}^{N-1}$ and $\omega \in \mathbb{R}$. For the “$x$-Fourier variable”, we use the notation $\xi \in \mathbb{R}^{N-1}$, while, for the “$t$-Laplace variable”, we reserve the notation $s > 0$, see Appendix A.3. We denote the upper half-space (the “field”) $\mathbb{R}_+^N := \mathbb{R}^{N-1} \times (0, +\infty)$, its boundary being the hyperplane (the “road”) $\partial \mathbb{R}_+^N = \mathbb{R}^{N-1} \times \{0\}$.

For given $\mu, \nu > 0$ and $D > d > 0$ (we refer to Remark 4.1 for the case $D \le d$), we thus consider the linear field-road problem (1.1) supplemented with the initial datum

$$ \left\{ \begin{array}{lll} v|_{t=0} = v_0, & X \in \mathbb{R}_+^N, & v_0 \in L^\infty(\mathbb{R}_+^N), \\ u|_{t=0} = u_0, & x \in \mathbb{R}^{N-1}, & u_0 \in L^\infty(\mathbb{R}^{N-1}). \end{array} \right. \qquad (2.1) $$

Our first main contribution is to provide an explicit expression of the solution to this Cauchy problem.

**Theorem 2.1** (Solution to the linear field-road Cauchy problem) Assume $D > d$. Then the unique bounded solution to the Cauchy problem (1.1)-(2.1) is

$$ v(t, X) = V(t, X) + \frac{\mu}{\sqrt{d}} \int_{\mathbb{R}^{N-1}} \Lambda(t, z, y) u_0(x - z) dz + \frac{\mu \nu}{\sqrt{d}} \int_0^t \int_{\mathbb{R}^{N-1}} \Lambda(s, z, y) V|_{y=0}(t - s, x - z) dz ds, \quad (2.2) $$

$$ u(t, x) = e^{-\mu t} U(t, x) + \nu \int_{0}^{t} e^{-\mu(t-s)} \int_{\mathbb{R}^{N-1}} G_{R}^{D}(t-s, x-z) v|_{y=0}(s, z) dz ds, \quad (2.3) $$

where

* $V = V(t, X)$ is the solution to the Cauchy problem

$$ \left\{ \begin{array}{lll} \partial_t V = d\Delta V, & t > 0, x \in \mathbb{R}^{N-1}, y > 0, \\ \nu V|_{y=0} - d\partial_y V|_{y=0} = 0, & t > 0, x \in \mathbb{R}^{N-1}, \\ V|_{t=0} = v_0, & x \in \mathbb{R}^{N-1}, y > 0, \end{array} \right. \qquad (2.4) $$

* $U = U(t, x)$ is the solution to the Cauchy problem

$$ \left\{ \begin{array}{ll} \partial_t U = D\Delta U, & t > 0, x \in \mathbb{R}^{N-1}, \\ U|_{t=0} = u_0, & x \in \mathbb{R}^{N-1}, \end{array} \right. \qquad (2.5) $$
Such a linear rate estimate is a key stone information to tackle difficult nonlinear issues. For instance, it is expected, see [16] or [3], that such a decay rate determines the so-called Fujita exponent $p_F$ for the system (1.2) when $f(v) = v^{1+p}$: this exponent separates “systematic blow-up” (when $0 < p \le p_F$) from “possible extinction” (when $p > p_F$). Also, understanding the Fujita blow-up phenomena when $f(v) = v^{1+p}$ typically enables, see [2], to solve the issue of “extinction vs propagation” in the population dynamics model (1.2) with $f(v) = v^{1+p}(1 - v)$. Such a nonlinearity serves a model for the so-called Allee effect [4], roughly meaning that the per capita growth rate of the population is not maximal at small density. We plan to address such issues in a future work.

## 2 Main results

Let $N \ge 2$ be an integer. A generic point $X \in \mathbb{R}^N$ will be written $X = (x, y)$ with $x \in \mathbb{R}^{N-1}$ and $y \in \mathbb{R}$. For the “variable of integration”, we reserve the notation $Z = (z, \omega) \in \mathbb{R}^N$, with $z \in \mathbb{R}^{N-1}$ and $\omega \in \mathbb{R}$. For the “x-Fourier variable”, we use the notation $\xi \in \mathbb{R}^{N-1}$, while, for the “t-Laplace variable”, we reserve the notation $s > 0$, see Appendix A.3. We denote the upper half-space (the “field”) $\mathbb{R}_+^N := \mathbb{R}^{N-1} \times (0, +\infty)$, its boundary being the hyperplane (the “road”) $\partial \mathbb{R}_+^N = \mathbb{R}^{N-1} \times \{0\}$.

For given $\mu, \nu > 0$ and $D > d > 0$ (we refer to Remark 4.1 for the case $D \le d$), we thus consider the linear field-road problem (1.1) supplemented with the initial datum

$$ \left\{ \begin{array}{lll} v|_{t=0} = v_0, & X \in \mathbb{R}_+^N, & v_0 \in L^\infty(\mathbb{R}_+^N), \\ u|_{t=0} = u_0, & x \in \mathbb{R}^{N-1}, & u_0 \in L^\infty(\mathbb{R}^{N-1}). \end{array} \right. \qquad (2.1) $$

Our first main contribution is to provide an explicit expression of the solution to this Cauchy problem.

**Theorem 2.1** (Solution to the linear field-road Cauchy problem) Assume $D > d$. Then the unique bounded solution to the Cauchy problem (1.1)-(2.1) is

$$ v(t, X) = V(t, X) + \frac{\mu}{\sqrt{d}} \int_{\mathbb{R}^{N-1}} \Lambda(t, z, y) u_0(x - z) dz + \frac{\mu \nu}{\sqrt{d}} \int_0^t \int_{\mathbb{R}^{N-1}} \Lambda(s, z, y) V|_{y=0}(t - s, x - z) dz ds, \quad (2.2) $$

$$ u(t, x) = e^{-\mu t} U(t, x) + \nu \int_{0}^{t} e^{-\mu(t-s)} \int_{\mathbb{R}^{N-1}} G_{R}^{D}(t-s, x-z) v|_{y=0}(s, z) dz ds, \quad (2.3) $$

where

* $V = V(t, X)$ is the solution to the Cauchy problem

$$ \left\{ \begin{array}{lll} \partial_t V = d\Delta V, & t > 0, x \in \mathbb{R}^{N-1}, y > 0, \\ \nu V|_{y=0} - d\partial_y V|_{y=0} = 0, & t > 0, x \in \mathbb{R}^{N-1}, \\ V|_{t=0} = v_0, & x \in \mathbb{R}^{N-1}, y > 0, \end{array} \right. \qquad (2.4) $$

* $U = U(t, x)$ is the solution to the Cauchy problem

$$ \left\{ \begin{array}{ll} \partial_t U = D\Delta U, & t > 0, x \in \mathbb{R}^{N-1}, \\ U|_{t=0} = u_0, & x \in \mathbb{R}^{N-1}, \end{array} \right. \qquad (2.5) $$
• $G_R^D = G_R^D(t, x)$ denotes the $(N-1)$-dimensional Heat kernel with diffusion D, that is

$$G_R^D(t, x) := \frac{1}{(4\pi Dt)^{\frac{N-1}{2}}} e^{-\frac{\|x\|^2}{4Dt}}, \qquad (2.6)$$

• and $\Lambda = \Lambda(t, x, y)$ is defined as

$$\Lambda(t, x, y) := \frac{e^{-\frac{y^2}{4dt}}}{(2\pi)^{N-1}} \int_{\mathbb{R}^{N-1}} \left[ a\alpha\Phi_{\alpha} + b\beta\Phi_{\beta} + c\gamma\Phi_{\gamma} \right] (t, \xi, y) e^{-dt\|\xi\|^2 + i\xi \cdot x} d\xi, \quad (2.7)$$

with $(\alpha, \beta, \gamma) = (\alpha, \beta, \gamma)(\xi)$ being the three complex roots of the $\delta$-indexed polynomials

$$P_{\delta}(\sigma) := \sigma^3 + \frac{\nu}{\sqrt{d}}\sigma^2 + (\mu + \delta)\sigma + \frac{\nu\delta}{\sqrt{d}}, \quad \text{with } \delta := (D - d)\|\xi\|^2, \qquad (2.8)$$

$(a, b, c) = (a, b, c)(\xi)$ being given by

$$a := \frac{1}{(\alpha - \beta)(\alpha - \gamma)}, \quad b := \frac{1}{(\beta - \alpha)(\beta - \gamma)}, \quad c := \frac{1}{(\gamma - \alpha)(\gamma - \beta)}, \qquad (2.9)$$

and for $\bullet \in \{\alpha, \beta, \gamma\}$,

$$\Phi_{\bullet}(t, \xi, y) := \frac{\text{Erfc}}{\Gamma} \left( \frac{-2 \cdot \sqrt{d} t + y}{2\sqrt{dt}} \right), \qquad (2.10)$$

where $\Gamma(\ell) := e^{-\ell^2}$, and Erfc is the complementary error function, whose definition is recalled in Appendix A.1.

Let us give some remarks on Theorem 2.1. First, the theorem gives the unique bounded solution. Indeed, there could be non-physical solutions that change sign and grow very fast at infinity: this fact was proved for the Heat equation by Tychonoff [32], see also [29, Chapter 9]. To get rid of such solutions, and so, to ask uniqueness, one needs to impose very loose restriction on the growth at infinity; here we choose boundedness for the sake of simplicity.

Next, the Cauchy problem (2.5) is nothing else than the Heat equation in the hyperplane $\mathbb{R}^{N-1}$ and both the expression and the decay rate of its solution $U = U(t, x)$ are very well-known. On the other hand, the Cauchy problem (2.4) is the Heat equation in the half-space $\mathbb{R}_+^N$ with Robin boundary conditions, and the understanding of its solution $V = V(t, X)$ is an important preliminary, that we recall in Section 3.

Last, but not least, the “migration-kernel” $\Lambda = \Lambda(t, x, y)$ is the keystone to write the solution to the field-road diffusion model. Its main role is to describe the in-flux of individuals in the field from the road via the last two terms of $v$ in (2.2), whose form evokes a Duhamel’s formula. We may give a physical sense to $\Lambda$ by remarking that it is the “$v$-component” of the solution to (1.1) starting from $(v_0, u_0) \equiv (0, \frac{\sqrt{d}}{\mu} \delta_{x=0})$. Notice also that, having in mind the acceleration phenomenon [10, 9] when $D > 2d$, $\Lambda$ is the only term in the expression of $v$ that involves the constant $D$. However, the expression of $\Lambda(t, x, y)$ is quite intricate. As a consequence of our Fourier/Laplace transform approach, it involves an integral over $\xi \in \mathbb{R}^{N-1}$. Furthermore, we show in Appendix A.2 that, for almost all $\xi \in \mathbb{R}^{N-1}$, the three complex roots $\alpha = \alpha(\xi)$, $\beta = \beta(\xi)$ and $\gamma = \gamma(\xi)$ of the polynomials $P_{\delta}$ are distinct, which allows to define
• $G_R^D = G_R^D(t, x)$ denotes the $(N-1)$-dimensional Heat kernel with diffusion D, that is

$$G_R^D(t, x) := \frac{1}{(4\pi Dt)^{\frac{N-1}{2}}} e^{-\frac{\|x\|^2}{4Dt}}, \qquad (2.6)$$

• and $\Lambda = \Lambda(t, x, y)$ is defined as

$$\Lambda(t, x, y) := \frac{e^{-\frac{y^2}{4dt}}}{(2\pi)^{N-1}} \int_{\mathbb{R}^{N-1}} \left[ a\alpha\Phi_{\alpha} + b\beta\Phi_{\beta} + c\gamma\Phi_{\gamma} \right] (t, \xi, y) e^{-dt\|\xi\|^2 + i\xi \cdot x} d\xi, \quad (2.7)$$

with $(\alpha, \beta, \gamma) = (\alpha, \beta, \gamma)(\xi)$ being the three complex roots of the $\delta$-indexed polynomials

$$P_{\delta}(\sigma) := \sigma^3 + \frac{\nu}{\sqrt{d}}\sigma^2 + (\mu + \delta)\sigma + \frac{\nu\delta}{\sqrt{d}}, \quad \text{with } \delta := (D - d)\|\xi\|^2, \quad (2.8)$$

$(a, b, c) = (a, b, c)(\xi)$ being given by

$$a := \frac{1}{(\alpha - \beta)(\alpha - \gamma)}, \quad b := \frac{1}{(\beta - \alpha)(\beta - \gamma)}, \quad c := \frac{1}{(\gamma - \alpha)(\gamma - \beta)}, \quad (2.9)$$

and for $\bullet \in \{\alpha, \beta, \gamma\}$,

$$\Phi_{\bullet}(t, \xi, y) := \frac{\text{Erfc}}{\Gamma} \left( \frac{-2 \cdot \sqrt{d} t + y}{2\sqrt{dt}} \right), \qquad (2.10)$$

where $\Gamma(\ell) := e^{-\ell^2}$, and Erfc is the complementary error function, whose definition is recalled in Appendix A.1.

Let us give some remarks on Theorem 2.1. First, the theorem gives the unique bounded solution. Indeed, there could be non-physical solutions that change sign and grow very fast at infinity: this fact was proved for the Heat equation by Tychonoff [32], see also [29, Chapter 9]. To get rid of such solutions, and so, to ask uniqueness, one needs to impose very loose restriction on the growth at infinity; here we choose boundedness for the sake of simplicity.

Next, the Cauchy problem (2.5) is nothing else than the Heat equation in the hyperplane $\mathbb{R}^{N-1}$ and both the expression and the decay rate of its solution $U = U(t, x)$ are very well-known. On the other hand, the Cauchy problem (2.4) is the Heat equation in the half-space $\mathbb{R}_+^N$ with Robin boundary conditions, and the understanding of its solution $V = V(t, X)$ is an important preliminary, that we recall in Section 3.

Last, but not least, the “migration-kernel” $\Lambda = \Lambda(t, x, y)$ is the keystone to write the solution to the field-road diffusion model. Its main role is to describe the in-flux of individuals in the field from the road via the last two terms of $v$ in (2.2), whose form evokes a Duhamel’s formula. We may give a physical sense to $\Lambda$ by remarking that it is the “$v$-component” of the solution to (1.1) starting from $(v_0, u_0) \equiv (0, \frac{\sqrt{d}}{\mu}\delta_{x=0})$. Notice also that, having in mind the acceleration phenomenon [10, 9] when $D > 2d$, $\Lambda$ is the only term in the expression of $v$ that involves the constant $D$. However, the expression of $\Lambda(t, x, y)$ is quite intricate. As a consequence of our Fourier/Laplace transform approach, it involves an integral over $\xi \in \mathbb{R}^{N-1}$. Furthermore, we show in Appendix A.2 that, for almost all $\xi \in \mathbb{R}^{N-1}$, the three complex roots $\alpha = \alpha(\xi)$, $\beta = \beta(\xi)$ and $\gamma = \gamma(\xi)$ of the polynomials $P_\delta$ are distinct, which allows to define
$a = a(\xi)$, $b = b(\xi)$ and $c = c(\xi)$ in (2.9). Further details will appear in Section 4 and 5 and in Appendix A.2.

Our second main contribution is to estimate the decay rate of the solution to the linear field-road system (1.1) starting from the datum, say

$$ \left\{ \begin{array}{lll} v|_{t=0} = v_0, & X \in \mathbb{R}_+^N, & v_0 \in L^\infty(\mathbb{R}_+^N) \text{ is nonnegative, and compactly supported,} \\ u|_{t=0} = u_0, & x \in \mathbb{R}^{N-1}, & u_0 \in L^\infty(\mathbb{R}^{N-1}) \text{ is nonnegative, and compactly supported.} \end{array} \right. \qquad (2.11) $$

In the statement below, the notation $B \lesssim B'$ means there is a constant $k = k(N, d, D, \nu, \mu) > 0$ such that $B \le kB'$.

**Theorem 2.2** (Decay rate of the solutions to the linear field-road system) Assume $D > d$ and let $(v, u)$ be the bounded solution to the Cauchy problem (1.1)–(2.11). Then

$$ \|v(t, \cdot)\|_{L^\infty(\mathbb{R}_+^N)} \lesssim \frac{C_{v_0} \ln(1+t) + C_{v_0}^{u_0}}{(1+t)^{\frac{N}{2}}}, \quad \forall t > 0, \qquad (2.12) $$

$$ \|u(t, \cdot)\|_{L^\infty(\mathbb{R}^{N-1})} \lesssim \frac{C_{v_0} \ln(1+t) + C_{v_0}^{u_0}}{(1+t)^{\frac{N}{2}}}, \quad \forall t > 0, \qquad (2.13) $$

for some nonnegative constant $C_{v_0}$ depending only on $v_0$, and some nonnegative constant $C_{v_0}^{u_0}$ depending on both $v_0$ and $u_0$.

The starting point to prove this upper estimate is the explicit expression of $v$ and $u$ given in Theorem 2.1. The main difficulty is to control the “migration-kernel” $\Lambda$ defined in (2.7) — (2.10) which, as explained above, is the core of the field-road model. To do so, we shall rely on a rather technical result, namely Lemma 5.1, whose proof needs a full understanding of the $\delta$-indexed polynomials $P_\delta$ in (2.8), achieved in Appendix A.2. Let us emphasize that the expression of $\Lambda$ we own, namely (2.7) — (2.10), is highly conditioned by the fact that, at some point, we use a decomposition into partial fractions to match up with known forms of Laplace-transforms. This decomposition makes appear the artefacts $a$, $b$ and $c$ which blow up as the roots of $P_\delta$ collide for singular values of $\delta$, see (2.9). Fortunately, due to the compensating structure of $a\alpha\Phi_\alpha + b\beta\Phi_\beta + c\gamma\Phi_\gamma$ and as confirmed by the technical Lemma 5.1, these singularities are “artificial”.

By following the proof in Section 5, one can track the values of $C_{v_0}$ and $C_{v_0}^{u_0}$ and check, in particular, that $v_0 \equiv 0$ implies $C_{v_0} = 0$. As a result, if there is initially no individual in the field, the estimate of the decay rate becomes of the magnitude $O((1+t)^{-N/2})$, which corresponds to the decay rate of the Heat equation in the half space $\mathbb{R}_+^N$ with Neumann boundary conditions — see Theorem 3.1. This is consistent with the fact that the boundary conditions of the field-road model are of the exchange type so that, roughly speaking, individuals are “not lost” but “stored” in the road. On the other hand, if $v_0 \ne 0$, our estimate of the decay rate becomes of the magnitude $O((1+t)^{-N/2}\ln(1+t))$, which is very slightly “less good”. We believe that the actual decay rate is always of the magnitude $O((1+t)^{-N/2})$: to deal with the terms $P(t, x, y)$ and $Q(t, x, y)$ appearing in (5.1), we rely on the uniform estimate (5.4) of the rather intricate function $\Phi(t, \xi, y)$ appearing in Lemma 5.1. The use of this uniform estimate is the cause of the logarithmic term appearing later, but refined estimates seem delicate to reach. Notice however that, in order to study the Fujita blow-up phenomenon, the logarithmic term should only play a (bad) role in the so-called *critical case*.
$a = a(\xi)$, $b = b(\xi)$ and $c = c(\xi)$ in (2.9). Further details will appear in Section 4 and 5 and in Appendix A.2.

Our second main contribution is to estimate the decay rate of the solution to the linear field-road system (1.1) starting from the datum, say

$$ \left\{ \begin{array}{lll} v|_{t=0} = v_0, & X \in \mathbb{R}_+^N, & v_0 \in L^\infty(\mathbb{R}_+^N) \text{ is nonnegative, and compactly supported,} \\ u|_{t=0} = u_0, & x \in \mathbb{R}^{N-1}, & u_0 \in L^\infty(\mathbb{R}^{N-1}) \text{ is nonnegative, and compactly supported.} \end{array} \right. \qquad (2.11) $$

In the statement below, the notation $B \lesssim B'$ means there is a constant $k = k(N, d, D, \nu, \mu) > 0$ such that $B \le kB'$.

**Theorem 2.2** (Decay rate of the solutions to the linear field-road system) Assume $D > d$ and let $(v, u)$ be the bounded solution to the Cauchy problem (1.1)–(2.11). Then

$$ \|v(t, \cdot)\|_{L^\infty(\mathbb{R}_+^N)} \lesssim \frac{C_{v_0} \ln(1+t) + C_{v_0}^{u_0}}{(1+t)^{\frac{N}{2}}}, \quad \forall t > 0, \qquad (2.12) $$

$$ \|u(t, \cdot)\|_{L^\infty(\mathbb{R}^{N-1})} \lesssim \frac{C_{v_0} \ln(1+t) + C_{v_0}^{u_0}}{(1+t)^{\frac{N}{2}}}, \quad \forall t > 0, \qquad (2.13) $$

for some nonnegative constant $C_{v_0}$ depending only on $v_0$, and some nonnegative constant $C_{v_0}^{u_0}$ depending on both $v_0$ and $u_0$.

The starting point to prove this upper estimate is the explicit expression of $v$ and $u$ given in Theorem 2.1. The main difficulty is to control the “migration-kernel” $\Lambda$ defined in (2.7) — (2.10) which, as explained above, is the core of the field-road model. To do so, we shall rely on a rather technical result, namely Lemma 5.1, whose proof needs a full understanding of the $\delta$-indexed polynomials $P_\delta$ in (2.8), achieved in Appendix A.2. Let us emphasize that the expression of $\Lambda$ we own, namely (2.7) — (2.10), is highly conditioned by the fact that, at some point, we use a decomposition into partial fractions to match up with known forms of Laplace-transforms. This decomposition makes appear the artefacts $a$, $b$ and $c$ which blow up as the roots of $P_\delta$ collide for singular values of $\delta$, see (2.9). Fortunately, due to the compensating structure of $a\alpha\Phi_\alpha + b\beta\Phi_\beta + c\gamma\Phi_\gamma$ and as confirmed by the technical Lemma 5.1, these singularities are “artificial”.

By following the proof in Section 5, one can track the values of $C_{v_0}$ and $C_{v_0}^{u_0}$ and check, in particular, that $v_0 \equiv 0$ implies $C_{v_0} = 0$. As a result, if there is initially no individual in the field, the estimate of the decay rate becomes of the magnitude $O((1+t)^{-N/2})$, which corresponds to the decay rate of the Heat equation in the half space $\mathbb{R}_+^N$ with Neumann boundary conditions — see Theorem 3.1. This is consistent with the fact that the boundary conditions of the field-road model are of the exchange type so that, roughly speaking, individuals are “not lost” but “stored” in the road. On the other hand, if $v_0 \ne 0$, our estimate of the decay rate becomes of the magnitude $O((1+t)^{-N/2}\ln(1+t))$, which is very slightly “less good”. We believe that the actual decay rate is always of the magnitude $O((1+t)^{-N/2})$: to deal with the terms $P(t, x, y)$ and $Q(t, x, y)$ appearing in (5.1), we rely on the uniform estimate (5.4) of the rather intricate function $\Phi(t, \xi, y)$ appearing in Lemma 5.1. The use of this uniform estimate is the cause of the logarithmic term appearing later, but refined estimates seem delicate to reach. Notice however that, in order to study the Fujita blow-up phenomenon, the logarithmic term should only play a (bad) role in the so-called *critical case*.
The organization of the present paper is as follows. In Section 3, we recall some facts about the fundamental solution of the Heat equation in a half-space, that will prove useful in the sequel. In particular, we give a complete derivation of the expression of the fundamental solution in the case of Robin boundary conditions, thus introducing the strategy of proof we will employ to study the complete system (1.1) in Section 4, where we compute the explicit form of the solutions of the linear field-road system, while their decay rate is estimated in Section 5. Last, in Section 6, we present some numerical simulations which not only “validate” Theorem 2.2 but also explore some open problems related to the level sets of the solution, or to the sign of the flux of individuals from the road to the field.

### 3 The Heat equation in the half-space

As an important preliminary, we discuss in this section the Heat equation in the half-space, in particular the explicit expression and the decay rate of the $L^\infty$ norm of the solution to the associated Cauchy problem with different boundary conditions. These facts are very classical for the Neumann and Dirichlet boundary conditions, but the details for the case of Robin boundary conditions are not so easy to find. We therefore take this opportunity to present the main arguments, by using a Fourier/Laplace transform strategy, which we use again in the next section to deal with the whole field-road system (1.1).

Let $N \ge 1$ be an integer. We recall that a generic point $X \in \mathbb{R}^N$ will be written $X = (x, y)$ with $x \in \mathbb{R}^{N-1}$ and $y \in \mathbb{R}$. For the “variable of integration”, we reserve the notation $Z = (z, \omega) \in \mathbb{R}^N$, with $z \in \mathbb{R}^{N-1}$ and $\omega \in \mathbb{R}$. For $d > 0$, we denote the $N$-dimensional Heat kernel with diffusion $d$,

$$G(t, X) := \frac{1}{(4\pi dt)^{\frac{N}{2}}} e^{-\frac{\|X\|^2}{4dt}} = \frac{1}{(4\pi dt)^{\frac{N}{2}}} e^{-\frac{\|x\|^2+y^2}{4dt}}, \quad (3.1)$$

and finally, we denote the upper half-space $\mathbb{R}_+^N := \mathbb{R}^{N-1} \times (0, +\infty)$, its boundary being the hyperplane $\partial \mathbb{R}_+^N = \mathbb{R}^{N-1} \times \{0\}$.

For $0 \le \theta \le 1$, we thus consider the Cauchy problem

$$\left\{
\begin{array}{ll}
\partial_t v = d\Delta v, & t > 0, X \in \mathbb{R}_+^N, \\
\theta v + (1 - \theta)d\partial_n v = 0, & t > 0, X \in \partial \mathbb{R}_+^N, \\
v|_{t=0} = v_0, & X \in \mathbb{R}_+^N,
\end{array}
\right.
\quad (3.2)$$

where, say, $v_0 \in L^\infty(\mathbb{R}_+^N)$. Here $n := (0_{\mathbb{R}^{N-1}}, -1)$ and the boundary condition may be recast

$$\theta v - (1 - \theta) d\partial_y v = 0, \quad t > 0, X \in \partial \mathbb{R}_+^N.$$

In this section, the notation $B \lesssim B'$ means there is a constant $k = k(N, d, \theta) > 0$ such that $B \le kB'$.

**Theorem 3.1 (Heat equation in the half-space) The solution to the Cauchy problem (3.2) is**

$$v(t, X) = \int_{\mathbb{R}_+^N} H_\theta(t, X, Z) v_0(Z) dZ, \quad (3.3)$$

for some $H_\theta(t, X, Z)$ precised below.
The organization of the present paper is as follows. In Section 3, we recall some facts about the fundamental solution of the Heat equation in a half-space, that will prove useful in the sequel. In particular, we give a complete derivation of the expression of the fundamental solution in the case of Robin boundary conditions, thus introducing the strategy of proof we will employ to study the complete system (1.1) in Section 4, where we compute the explicit form of the solutions of the linear field-road system, while their decay rate is estimated in Section 5. Last, in Section 6, we present some numerical simulations which not only “validate” Theorem 2.2 but also explore some open problems related to the level sets of the solution, or to the sign of the flux of individuals from the road to the field.

### 3 The Heat equation in the half-space

As an important preliminary, we discuss in this section the Heat equation in the half-space, in particular the explicit expression and the decay rate of the $L^\infty$ norm of the solution to the associated Cauchy problem with different boundary conditions. These facts are very classical for the Neumann and Dirichlet boundary conditions, but the details for the case of Robin boundary conditions are not so easy to find. We therefore take this opportunity to present the main arguments, by using a Fourier/Laplace transform strategy, which we use again in the next section to deal with the whole field-road system (1.1).

Let $N \ge 1$ be an integer. We recall that a generic point $X \in \mathbb{R}^N$ will be written $X = (x, y)$ with $x \in \mathbb{R}^{N-1}$ and $y \in \mathbb{R}$. For the “variable of integration”, we reserve the notation $Z = (z, \omega) \in \mathbb{R}^N$, with $z \in \mathbb{R}^{N-1}$ and $\omega \in \mathbb{R}$. For $d > 0$, we denote the $N$-dimensional Heat kernel with diffusion $d$,

$$G(t, X) := \frac{1}{(4\pi dt)^{\frac{N}{2}}} e^{-\frac{\|X\|^2}{4dt}} = \frac{1}{(4\pi dt)^{\frac{N}{2}}} e^{-\frac{\|x\|^2+y^2}{4dt}}, \quad (3.1)$$

and finally, we denote the upper half-space $\mathbb{R}_+^N := \mathbb{R}^{N-1} \times (0, +\infty)$, its boundary being the hyperplane $\partial \mathbb{R}_+^N = \mathbb{R}^{N-1} \times \{0\}$.

For $0 \le \theta \le 1$, we thus consider the Cauchy problem

$$\left\{
\begin{array}{ll}
\partial_t v = d\Delta v, & t > 0, X \in \mathbb{R}_+^N, \\
\theta v + (1 - \theta)d\partial_n v = 0, & t > 0, X \in \partial \mathbb{R}_+^N, \\
v|_{t=0} = v_0, & X \in \mathbb{R}_+^N,
\end{array}
\right.
\quad (3.2)$$

where, say, $v_0 \in L^\infty(\mathbb{R}_+^N)$. Here $n := (0_{\mathbb{R}^{N-1}}, -1)$ and the boundary condition may be recast

$$\theta v - (1 - \theta) d\partial_y v = 0, \quad t > 0, X \in \partial \mathbb{R}_+^N.$$

In this section, the notation $B \lesssim B'$ means there is a constant $k = k(N, d, \theta) > 0$ such that $B \le kB'$.

**Theorem 3.1 (Heat equation in the half-space) The solution to the Cauchy problem (3.2) is**

$$v(t, X) = \int_{\mathbb{R}_+^N} H_\theta(t, X, Z) v_0(Z) dZ, \quad (3.3)$$

for some $H_\theta(t, X, Z)$ precised below.
(i) When $\theta = 0$ (Neumann case), we have

$$H_0(t, X, Z) = G(t, x-z, y-\omega) + G(t, x-z, y+\omega),$$

and

$$\|v(t, \cdot)\|_{L^{\infty}(\mathbb{R}_+^N)} \lesssim \frac{\int_{\mathbb{R}_+^N} |v_0(Z)| dZ}{t^{\frac{N}{2}}}, \quad \forall t > 0.$$

(ii) When $\theta = 1$ (Dirichlet case), we have

$$H_1(t, X, Z) = G(t, x-z, y-\omega) - G(t, x-z, y+\omega),$$

and

$$\|v(t, \cdot)\|_{L^{\infty}(\mathbb{R}_+^N)} \lesssim \frac{\int_{\mathbb{R}_+^N} \omega |v_0(Z)| dZ}{t^{\frac{N+1}{2}}}, \quad \forall t > 0.$$

(iii) When $0 < \theta < 1$ (Robin case), we have

$$H_{\theta}(t, X, Z) = H_{0}(t, X, Z) - 2\sqrt{\pi} A \sqrt{dt} G(t, x-z, y+\omega) \frac{\text{Erfc}\left(\frac{2A dt + y + \omega}{2\sqrt{dt}}\right)}{\Gamma}, \quad (3.4)$$

$$= H_1(t, X, Z) + 2G(t, x-z, y+\omega) \left(1 - \sqrt{\pi} A \sqrt{dt} \frac{\text{Erfc}\left(\frac{2A dt + y + \omega}{2\sqrt{dt}}\right)}{\Gamma}\right), \quad (3.5)$$

where

$$A := \frac{\theta}{d(1-\theta)},$$

$\Gamma(\ell) := e^{-\ell^2}$, and the definition of Erfc is recalled in Appendix A.1. Furthermore,

$$\|v(t, \cdot)\|_{L^{\infty}(\mathbb{R}_+^N)} \lesssim \frac{\int_{\mathbb{R}_+^N} (1+\omega) |v_0(Z)| dZ}{t^{\frac{N+1}{2}}}, \quad \forall t > 0, \quad (3.6)$$

and

$$\|v(t, \cdot)\|_{L^{\infty}(\mathbb{R}_+^N)} \lesssim \frac{C_{v_0}^{\#}}{(1+t)^{\frac{N+1}{2}}}, \quad \forall t > 0, \quad (3.7)$$

$$\text{where } C_{v_0}^{\#} := \left( \|v_0\|_{L^{\infty}(\mathbb{R}_+^N)}^{\frac{2}{N+1}} + \left( \int_{\mathbb{R}_+^N} (1+\omega) |v_0(Z)| dZ \right)^{\frac{2}{N+1}} \right)^{\frac{N+1}{2}}.$$

This theorem is rather classical, in particular points (i) and (ii). A possible way to prove it is to use some continuation arguments. For instance, to prove (i), one defines $\tilde{v}(t, X)$ as the solution of the Heat equation in the whole of $\mathbb{R}^N$ arising from the initial datum

$$\tilde{v}_0(x, y) := \begin{cases} v_0(x, y) & \text{if } y > 0, \\ v_0(x, -y) & \text{if } y < 0. \end{cases}$$

Then, by uniqueness of the (bounded) solutions of the Heat equation, $\tilde{v}$ remains even with respect to $y$, is $C^\infty$ in $\mathbb{R}^N \times (0, +\infty)$ by usual regularity results on parabolic equations. Thus
(i) When $\theta = 0$ (Neumann case), we have

$$H_0(t, X, Z) = G(t, x-z, y-\omega) + G(t, x-z, y+\omega),$$

and

$$\|v(t, \cdot)\|_{L^{\infty}(\mathbb{R}_+^N)} \lesssim \frac{\int_{\mathbb{R}_+^N} |v_0(Z)| dZ}{t^{\frac{N}{2}}}, \quad \forall t > 0.$$

(ii) When $\theta = 1$ (Dirichlet case), we have

$$H_1(t, X, Z) = G(t, x-z, y-\omega) - G(t, x-z, y+\omega),$$

and

$$\|v(t, \cdot)\|_{L^{\infty}(\mathbb{R}_+^N)} \lesssim \frac{\int_{\mathbb{R}_+^N} \omega |v_0(Z)| dZ}{t^{\frac{N+1}{2}}}, \quad \forall t > 0.$$

(iii) When $0 < \theta < 1$ (Robin case), we have

$$H_{\theta}(t, X, Z) = H_{0}(t, X, Z) - 2\sqrt{\pi} A \sqrt{dt} G(t, x-z, y+\omega) \frac{\text{Erfc}\left(\frac{2A dt + y + \omega}{2\sqrt{dt}}\right)}{\Gamma}, \quad (3.4)$$

$$= H_1(t, X, Z) + 2G(t, x-z, y+\omega) \left(1 - \sqrt{\pi} A \sqrt{dt} \frac{\text{Erfc}\left(\frac{2A dt + y + \omega}{2\sqrt{dt}}\right)}{\Gamma}\right), \quad (3.5)$$

where

$$A := \frac{\theta}{d(1-\theta)},$$

$\Gamma(\ell) := e^{-\ell^2}$, and the definition of Erfc is recalled in Appendix A.1. Furthermore,

$$\|v(t, \cdot)\|_{L^{\infty}(\mathbb{R}_+^N)} \lesssim \frac{\int_{\mathbb{R}_+^N} (1+\omega) |v_0(Z)| dZ}{t^{\frac{N+1}{2}}}, \quad \forall t > 0, \quad (3.6)$$

and

$$\|v(t, \cdot)\|_{L^{\infty}(\mathbb{R}_+^N)} \lesssim \frac{C_{v_0}^{\#}}{(1+t)^{\frac{N+1}{2}}}, \quad \forall t > 0, \quad (3.7)$$

$$\text{where } C_{v_0}^{\#} := \left( \|v_0\|_{L^{\infty}(\mathbb{R}_+^N)}^{\frac{2}{N+1}} + \left( \int_{\mathbb{R}_+^N} (1+\omega) |v_0(Z)| dZ \right)^{\frac{2}{N+1}} \right)^{\frac{N+1}{2}}.$$

This theorem is rather classical, in particular points (i) and (ii). A possible way to prove it is to use some continuation arguments. For instance, to prove (i), one defines $\tilde{v}(t, X)$ as the solution of the Heat equation in the whole of $\mathbb{R}^N$ arising from the initial datum

$$\tilde{v}_0(x, y) := \begin{cases} v_0(x, y) & \text{if } y > 0, \\ v_0(x, -y) & \text{if } y < 0. \end{cases}$$

Then, by uniqueness of the (bounded) solutions of the Heat equation, $\tilde{v}$ remains even with respect to $y$, is $C^\infty$ in $\mathbb{R}^N \times (0, +\infty)$ by usual regularity results on parabolic equations. Thus
it satisfies $\partial_y \tilde{v}(t, x, 0) = 0$ for all $t > 0$ and all $x \in \mathbb{R}^{N-1}$. Hence, $v(t, x, y)$ stands as the restriction of $\tilde{v}(t, x, y)$ to the upper-half-space, and we have

$$ 
\begin{aligned}
v(t, X) &= \int_{\mathbb{R}^N} G(t, X-Z) \tilde{v}_0(Z) dZ \\
&= \int_{\mathbb{R}_+^N} [G(t, x-z, y-\omega) + G(t, x-z, y+\omega)] v_0(z, \omega) dZ,
\end{aligned}
$$

from which we deduce the expression of $H_0$ and the decay rate. The point (ii) can be proven in the same way by considering this time, for $\tilde{v}_0$, the odd-continuation (with respect to $y$) of $v_0$. A more tricky (and less known) continuation argument can be used to reach the point (iii). However, in view of the sequel, let us present now the detailed proof of (iii) using rather a Fourier/Laplace transform strategy.

In the following, the $x$-Fourier transform of $u$ shall be denoted by $\hat{u}$, the $t$-Laplace one by $\hat{u}$ and the $x$-Fourier/$t$-Laplace one by $\hat{u}$. We refer to Appendix A.3 for conventions and notations related to these transforms.

*Proof of Theorem 3.1, (iii).* Applying the $t$-Laplace transform to the first equation in (3.2), namely $\partial_t v(t, x, y) = d(\Delta_x + \partial_{yy})v(t, x, y)$, yields the elliptic PDE

$$ d(\Delta_x \hat{v}(s, x, y) + \partial_{yy} \hat{v}(s, x, y)) - s\hat{v}(s, x, y) = -v_0(x, y), \quad (x, y) \in \mathbb{R}_+^N, $$

where $s > 0$ acts as a parameter. Then by applying the $x$-Fourier transform, we reach the linear second order ODE

$$ d\partial_{yy}\hat{v}(s, \xi, y) - \left(s + d\|\xi\|^2\right)\hat{v}(s, \xi, y) = -\hat{v}_0(\xi, y), \quad y > 0, \quad (3.8) $$

where $s > 0$ and $\xi \in \mathbb{R}^{N-1}$ act as parameters. Setting $\sigma := \sqrt{s + d\|\xi\|^2}$, solutions of (3.8) may be written

$$ \hat{v}(s, \xi, y) = e^{\frac{\sigma}{\sqrt{d}} y} \left( C_1 - \frac{1}{2\sqrt{d}\sigma} \int_0^y e^{-\frac{\sigma}{\sqrt{d}}\omega} \hat{v}_0(\xi, \omega) d\omega \right) + e^{-\frac{\sigma}{\sqrt{d}} y} \left( C_2 + \frac{1}{2\sqrt{d}\sigma} \int_0^y e^{\frac{\sigma}{\sqrt{d}}\omega} \hat{v}_0(\xi, \omega) d\omega \right), $$

where $C_1$ and $C_2$ are constants (with respect to $y$) to be determined. In order to insure $\lim_{y \to +\infty} \hat{v}(s, \xi, y) = 0$, one needs

$$ C_1 = \frac{1}{2\sqrt{d}\sigma} \int_0^{+\infty} e^{-\frac{\sigma}{\sqrt{d}}\omega} \hat{v}_0(\xi, \omega) d\omega. \quad (3.9) $$

Next, since

$$ \hat{v}(s, \xi, 0) = C_1 + C_2 \quad \text{and} \quad \partial_y \hat{v}(s, \xi, 0) = \frac{\sigma}{\sqrt{d}}(C_1 - C_2), $$

the Robin boundary conditions enforce

$$ C_2 = \left( 1 - \frac{2A\sqrt{d}}{\sigma + A\sqrt{d}} \right) C_1. \quad (3.10) $$
it satisfies $\partial_y \tilde{v}(t, x, 0) = 0$ for all $t > 0$ and all $x \in \mathbb{R}^{N-1}$. Hence, $v(t, x, y)$ stands as the restriction of $\tilde{v}(t, x, y)$ to the upper-half-space, and we have

$$ 
\begin{aligned}
v(t, X) &= \int_{\mathbb{R}^N} G(t, X - Z) \tilde{v}_0(Z) dZ \\
&= \int_{\mathbb{R}_+^N} [G(t, x - z, y - \omega) + G(t, x - z, y + \omega)] v_0(z, \omega) dZ,
\end{aligned}
$$

from which we deduce the expression of $H_0$ and the decay rate. The point (ii) can be proven in the same way by considering this time, for $\tilde{v}_0$, the odd-continuation (with respect to $y$) of $v_0$. A more tricky (and less known) continuation argument can be used to reach the point (iii). However, in view of the sequel, let us present now the detailed proof of (iii) using rather a Fourier/Laplace transform strategy.

In the following, the $x$-Fourier transform of $u$ shall be denoted by $\hat{u}$, the $t$-Laplace one by $\hat{u}$ and the $x$-Fourier/$t$-Laplace one by $\hat{u}$. We refer to Appendix A.3 for conventions and notations related to these transforms.

*Proof of Theorem 3.1, (iii).* Applying the $t$-Laplace transform to the first equation in (3.2), namely $\partial_t v(t, x, y) = d(\Delta_x + \partial_{yy})v(t, x, y)$, yields the elliptic PDE

$$ d(\Delta_x \hat{v}(s, x, y) + \partial_{yy} \hat{v}(s, x, y)) - s\hat{v}(s, x, y) = -v_0(x, y), \quad (x, y) \in \mathbb{R}_+^N, $$

where $s > 0$ acts as a parameter. Then by applying the $x$-Fourier transform, we reach the linear second order ODE

$$ d\partial_{yy}\hat{v}(s, \xi, y) - \left(s + d\|\xi\|^2\right)\hat{v}(s, \xi, y) = -\hat{v}_0(\xi, y), \quad y > 0, \quad (3.8) $$

where $s > 0$ and $\xi \in \mathbb{R}^{N-1}$ act as parameters. Setting $\sigma := \sqrt{s + d\|\xi\|^2}$, solutions of (3.8) may be written

$$ \hat{v}(s, \xi, y) = e^{\frac{\sigma}{\sqrt{d}} y} \left( C_1 - \frac{1}{2\sqrt{d}\sigma} \int_0^y e^{-\frac{\sigma}{\sqrt{d}}\omega} \hat{v}_0(\xi, \omega) d\omega \right) + e^{-\frac{\sigma}{\sqrt{d}} y} \left( C_2 + \frac{1}{2\sqrt{d}\sigma} \int_0^y e^{\frac{\sigma}{\sqrt{d}}\omega} \hat{v}_0(\xi, \omega) d\omega \right), $$

where $C_1$ and $C_2$ are constants (with respect to $y$) to be determined. In order to insure $\lim_{y \to +\infty} \hat{v}(s, \xi, y) = 0$, one needs

$$ C_1 = \frac{1}{2\sqrt{d}\sigma} \int_0^{+\infty} e^{-\frac{\sigma}{\sqrt{d}}\omega} \hat{v}_0(\xi, \omega) d\omega. \quad (3.9) $$

Next, since

$$ \hat{v}(s, \xi, 0) = C_1 + C_2 \quad \text{and} \quad \partial_y \hat{v}(s, \xi, 0) = \frac{\sigma}{\sqrt{d}}(C_1 - C_2), $$

the Robin boundary conditions enforce

$$ C_2 = \left( 1 - \frac{2A\sqrt{d}}{\sigma + A\sqrt{d}} \right) C_1. \quad (3.10) $$
This brings us to the following explicit expression for $\hat{v}$:

$$ \hat{v}(s, \xi, y) = \frac{1}{2\sqrt{d}} \int_{0}^{+\infty} \left( \frac{e^{-\frac{\sigma}{\sqrt{d}}|y-\omega|}}{\sigma} + \frac{e^{-\frac{\sigma}{\sqrt{d}}(y+\omega)}}{\sigma} - 2A\sqrt{d} \frac{e^{-\frac{\sigma}{\sqrt{d}}(y+\omega)}}{\sigma(\sigma + A\sqrt{d})} \right) \hat{v}_0(\xi, \omega) d\omega, $$

where we recognize known forms of t-Laplace transforms evaluated at $\sigma^2 = s + d\|\xi\|^2$ — see Lemma A.6, (iii) and (iv) — so that

$$ \hat{v}(s, \xi, y) = \frac{1}{2\sqrt{d}} \int_{0}^{+\infty} \mathcal{L} \left[ t \mapsto \frac{e^{-\frac{(y-\omega)^2}{4dt}}}{\sqrt{\pi t}} + \frac{e^{-\frac{(y+\omega)^2}{4dt}}}{\sqrt{\pi t}} - 2A\sqrt{d} \frac{\text{Erfc}}{\Gamma} \left( \frac{2Adt + y + \omega}{2\sqrt{dt}} \right) e^{-\frac{(y+\omega)^2}{4dt}} \right] (s + d\|\xi\|^2) \hat{v}_0(\xi, \omega) d\omega. $$

Using Lemma A.7 (iv), we can now get rid of the t-Laplace transform and get

$$ \hat{v}(t, \xi, y) = \frac{1}{2\sqrt{d}} \int_{0}^{+\infty} \left[ \frac{e^{-\frac{(y-\omega)^2}{4dt}}}{\sqrt{\pi t}} + \frac{e^{-\frac{(y+\omega)^2}{4dt}}}{\sqrt{\pi t}} \left( 1 - 2\sqrt{\pi} A \sqrt{dt} \frac{\text{Erfc}}{\Gamma} \left( \frac{2Adt + y + \omega}{2\sqrt{dt}} \right) \right) \right] e^{-dt\|\xi\|^2} \hat{v}_0(\xi, \omega) d\omega. $$

Finally, note that

$$ e^{-dt\|\xi\|^2} \hat{v}_0(\xi, \omega) = \mathcal{F} \left[ \frac{e^{-\frac{\|\cdot\|^2}{4dt}}}{(4\pi dt)^{\frac{N-1}{2}}} \right] (\xi) \times \mathcal{F} [v_0(\cdot, \omega)] (\xi) = \mathcal{F} \left[ \frac{e^{-\frac{\|\cdot\|^2}{4dt}}}{(4\pi dt)^{\frac{N-1}{2}}} * v_0(\cdot, \omega) \right] (\xi), $$

whence we can also get rid of the x-Fourier transform and reach

$$ v(t, x, y) = \int_{\mathbb{R}_+^N} \left[ \frac{e^{-\frac{\|x-z\|^2}{4dt}}}{(4\pi dt)^{\frac{N-1}{2}}} \left( \frac{e^{-\frac{(y-\omega)^2}{4dt}}}{\sqrt{4\pi dt}} + \frac{e^{-\frac{(y+\omega)^2}{4dt}}}{\sqrt{4\pi dt}} \left( 1 - 2\sqrt{\pi} A \sqrt{dt} \frac{\text{Erfc}}{\Gamma} \left( \frac{2Adt + y + \omega}{2\sqrt{dt}} \right) \right) \right) \right] v_0(Z) \, dZ. \quad (3.11) $$

Expressions (3.4) and (3.5) of $H_\theta$ are immediately deduced from (3.11).

We now turn to the decay rate. It is readily seen that

$$ -\sqrt{\pi} \frac{\text{Erfc}}{\Gamma}(\ell) \leq -\frac{1}{1+\ell}, \quad \forall \ell \geq 0. $$

As a consequence, for all $t > 0$, all $y \ge 0$, and all $\omega \ge 0$,

$$ -\sqrt{\pi} \frac{\text{Erfc}}{\Gamma} \left( \frac{2Adt + y + \omega}{2\sqrt{dt}} \right) \leq -\frac{2\sqrt{dt}}{2Adt + 2\sqrt{dt} + y + \omega}. $$

Using this into (3.5), we see that, for all $t > 0$, all $X \in \mathbb{R}_+^N$, and all $Z \in \mathbb{R}_+^N$,

$$ \begin{aligned} 0 < H_{\theta}(t, X, Z) - H_{1}(t, X, Z) & \le 2G(t, x-z, y+\omega) \frac{2\sqrt{dt} + y + \omega}{2Adt + 2\sqrt{dt} + y + \omega} \\ & \lesssim \frac{m(t)}{t^{\frac{N}{2}}}, \end{aligned} \quad (3.12) $$
This brings us to the following explicit expression for $\hat{v}$:

$$ \hat{v}(s, \xi, y) = \frac{1}{2\sqrt{d}} \int_{0}^{+\infty} \left( \frac{e^{-\frac{\sigma}{\sqrt{d}}|y-\omega|}}{\sigma} + \frac{e^{-\frac{\sigma}{\sqrt{d}}(y+\omega)}}{\sigma} - 2A\sqrt{d} \frac{e^{-\frac{\sigma}{\sqrt{d}}(y+\omega)}}{\sigma(\sigma + A\sqrt{d})} \right) \hat{v}_0(\xi, \omega) d\omega, $$

where we recognize known forms of t-Laplace transforms evaluated at $\sigma^2 = s + d\|\xi\|^2$ — see Lemma A.6, (iii) and (iv) — so that

$$ \hat{v}(s, \xi, y) = \frac{1}{2\sqrt{d}} \int_{0}^{+\infty} \mathcal{L} \left[ t \mapsto \frac{e^{-\frac{(y-\omega)^2}{4dt}}}{\sqrt{\pi t}} + \frac{e^{-\frac{(y+\omega)^2}{4dt}}}{\sqrt{\pi t}} - 2A\sqrt{d} \frac{\text{Erfc}}{\Gamma} \left( \frac{2Adt + y + \omega}{2\sqrt{dt}} \right) e^{-\frac{(y+\omega)^2}{4dt}} \right] (s + d\|\xi\|^2) \hat{v}_0(\xi, \omega) d\omega. $$

Using Lemma A.7 (iv), we can now get rid of the t-Laplace transform and get

$$ \hat{v}(t, \xi, y) = \frac{1}{2\sqrt{d}} \int_{0}^{+\infty} \left[ \frac{e^{-\frac{(y-\omega)^2}{4dt}}}{\sqrt{\pi t}} + \frac{e^{-\frac{(y+\omega)^2}{4dt}}}{\sqrt{\pi t}} \left( 1 - 2\sqrt{\pi} A \sqrt{dt} \frac{\text{Erfc}}{\Gamma} \left( \frac{2Adt + y + \omega}{2\sqrt{dt}} \right) \right) \right] e^{-dt\|\xi\|^2} \hat{v}_0(\xi, \omega) d\omega. $$

Finally, note that

$$ e^{-dt\|\xi\|^2} \hat{v}_0(\xi, \omega) = \mathcal{F} \left[ \frac{e^{-\frac{\|\cdot\|^2}{4dt}}}{(4\pi dt)^{\frac{N-1}{2}}} \right] (\xi) \times \mathcal{F} [v_0(\cdot, \omega)](\xi) = \mathcal{F} \left[ \frac{e^{-\frac{\|\cdot\|^2}{4dt}}}{(4\pi dt)^{\frac{N-1}{2}}} * v_0(\cdot, \omega) \right] (\xi), $$

whence we can also get rid of the x-Fourier transform and reach

$$ v(t, x, y) = \int_{\mathbb{R}_+^N} \left[ \frac{e^{-\frac{\|x-z\|^2}{4dt}}}{(4\pi dt)^{\frac{N-1}{2}}} \left( \frac{e^{-\frac{(y-\omega)^2}{4dt}}}{\sqrt{4\pi dt}} + \frac{e^{-\frac{(y+\omega)^2}{4dt}}}{\sqrt{4\pi dt}} \left( 1 - 2\sqrt{\pi} A \sqrt{dt} \frac{\text{Erfc}}{\Gamma} \left( \frac{2Adt + y + \omega}{2\sqrt{dt}} \right) \right) \right) \right] v_0(Z) \, dZ. \quad (3.11) $$

Expressions (3.4) and (3.5) of $H_\theta$ are immediately deduced from (3.11).

We now turn to the decay rate. It is readily seen that

$$ -\sqrt{\pi} \frac{\text{Erfc}}{\Gamma}(\ell) \le -\frac{1}{1+\ell}, \quad \forall \ell \ge 0. $$

As a consequence, for all $t > 0$, all $y \ge 0$, and all $\omega \ge 0$,

$$ -\sqrt{\pi} \frac{\text{Erfc}}{\Gamma} \left( \frac{2Adt + y + \omega}{2\sqrt{dt}} \right) \le -\frac{2\sqrt{dt}}{2Adt + 2\sqrt{dt} + y + \omega}. $$

Using this into (3.5), we see that, for all $t > 0$, all $X \in \mathbb{R}_+^N$, and all $Z \in \mathbb{R}_+^N$,

$$ \begin{aligned} 0 < H_{\theta}(t, X, Z) - H_{1}(t, X, Z) &\le 2G(t, x-z, y+\omega) \frac{2\sqrt{dt} + y + \omega}{2Adt + 2\sqrt{dt} + y + \omega} \\ &\lesssim \frac{m(t)}{t^{\frac{N}{2}}}, \end{aligned} \quad (3.12) $$
where

$$m(t) := \sup_{r \ge 0} \left( \frac{2\sqrt{dt} + r}{2Adt + 2\sqrt{dt} + r} e^{-\frac{r^2}{4dt}} \right).$$

An elementary analysis shows that

$$m(t) = \max_{0 \le r \le \sqrt{dt}} \left( \frac{2\sqrt{dt} + r}{2Adt + 2\sqrt{dt} + r} e^{-\frac{r^2}{4dt}} \right) \le \max_{0 \le r \le \sqrt{dt}} \left( \frac{2\sqrt{dt} + r}{2Adt + 2\sqrt{dt} + r} \right) \lesssim \frac{1}{\sqrt{t}}.$$

From this, (3.12), and the case (ii), we obtain the decay rate as stated in (3.6). Combining (3.6) and the fact that $\|v(t, \cdot)\|_{L^{\infty}(\mathbb{R}_+^N)} \le \|v_0\|_{L^{\infty}(\mathbb{R}_+^N)}$, one can check that (3.7) holds true. $\square$

Let us conclude this section with the following direct consequence of Theorem 3.1.

**Corollary 3.2** For any $0 \le \theta \le 1$, we have

$$H_{\theta}(t, X, Z) = \frac{e^{-\frac{\|x-z\|^2}{4dt}}}{(4\pi dt)^{\frac{N-1}{2}}} \times H_{\theta}^{(1)}(t, y, \omega),$$

where $H_{\theta}^{(1)}$ denotes $H_{\theta}$ for $N = 1$.

# 4 The solution to the field-road model Cauchy problem

In this section, we explicitly obtain the bounded solution to the Cauchy problem (1.1)–(2.1). Before starting, we recall that the $x$-Fourier transform of $u$ is denoted by $\hat{u}$, the $t$-Laplace one by $\hat{u}$ and the $x$-Fourier/$t$-Laplace one by $\hat{u}$. We refer to Appendix A.3 for conventions and notations related to these transforms.

**Proof of Theorem 2.1.** We start with establishing expression (2.3) for $u = u(t, x)$ which is the simplest to get. To do so, we apply the $t$-Laplace and then the $x$-Fourier transforms to the third equation in (1.1), namely $\partial_t u(t, x) = D\Delta_x u(t, x) + \nu v|_{y=0}(t, x) - \mu u(t, x)$. This yields the algebraic relation

$$-D\|\xi\|^2 \hat{u}(s, \xi) + \nu \hat{v}|_{y=0}(s, \xi) - \mu \hat{u}(s, \xi) = -\hat{u}_0(\xi) + s \hat{u}(s, \xi), \quad s > 0, \xi \in \mathbb{R}^{N-1},$$

which is recast

$$\hat{u}(s, \xi) = \frac{\hat{u}_0(\xi) + \nu \hat{v}|_{y=0}(s, \xi)}{\mu + s + D \|\xi\|^2}, \quad s > 0, \xi \in \mathbb{R}^{N-1}. \quad (4.1)$$

It is now possible to reach (2.3) from (4.1) by using properties and known forms of $t$-Laplace and $x$-Fourier transforms — see Lemmas A.6, (i) and (ii), and A.7, (ii) and (iii):

$$\begin{align*}
\hat{u}(s, \xi) &= \mathcal{L} \left[ t \mapsto e^{-(\mu+D\|\xi\|^2)t} \hat{u}_0(\xi) \right] (s) + \nu \mathcal{L} \left[ t \mapsto e^{-(\mu+D\|\xi\|^2)t} \right] (s) \times \mathcal{L} \left[ t \mapsto \hat{v}|_{y=0}(t, \xi) \right] (s) \\
&= \mathcal{L} \left[ t \mapsto e^{-(\mu+D\|\xi\|^2)t} \hat{u}_0(\xi) \right] (s) + \nu \mathcal{L} \left[ t \mapsto \int_0^t e^{-(\mu+D\|\xi\|^2)(t-\tau)} \hat{v}|_{y=0}(\tau, \xi) \, d\tau \right] (s) \\
&= \mathcal{LF} \left[ (t, x) \mapsto e^{-\mu t} [u_0 * G_R^D(t, \cdot)](x) \right] (s, \xi) \\
&\quad + \nu \mathcal{LF} \left[ (t, x) \mapsto \int_0^t e^{-\mu(t-\tau)} [v|_{y=0}(\tau, \cdot) * G_R^D(t-\tau, \cdot)](x) \, d\tau \right] (s, \xi) \\
&= \mathcal{LF} \left[ (t, x) \mapsto e^{-\mu t} \left( U(t, x) + \nu \int_0^t e^{\mu\tau} \int_{\mathbb{R}^{N-1}} G_R^D(t-\tau, x-z) v|_{y=0}(\tau, z) \, dz \, d\tau \right) \right] (s, \xi)
\end{align*}$$
where

$$m(t) := \sup_{r \ge 0} \left( \frac{2\sqrt{dt} + r}{2Adt + 2\sqrt{dt} + r} e^{-\frac{r^2}{4dt}} \right).$$

An elementary analysis shows that

$$m(t) = \max_{0 \le r \le \sqrt{dt}} \left( \frac{2\sqrt{dt} + r}{2Adt + 2\sqrt{dt} + r} e^{-\frac{r^2}{4dt}} \right) \le \max_{0 \le r \le \sqrt{dt}} \left( \frac{2\sqrt{dt} + r}{2Adt + 2\sqrt{dt} + r} \right) \lesssim \frac{1}{\sqrt{t}}.$$

From this, (3.12), and the case (ii), we obtain the decay rate as stated in (3.6). Combining (3.6) and the fact that $\|v(t, \cdot)\|_{L^{\infty}(\mathbb{R}_+^N)} \le \|v_0\|_{L^{\infty}(\mathbb{R}_+^N)}$, one can check that (3.7) holds true. $\square$

Let us conclude this section with the following direct consequence of Theorem 3.1.

**Corollary 3.2** For any $0 \le \theta \le 1$, we have

$$H_{\theta}(t, X, Z) = \frac{e^{-\frac{\|x-z\|^2}{4dt}}}{(4\pi dt)^{\frac{N-1}{2}}} \times H_{\theta}^{(1)}(t, y, \omega),$$

where $H_{\theta}^{(1)}$ denotes $H_{\theta}$ for $N = 1$.

# 4 The solution to the field-road model Cauchy problem

In this section, we explicitly obtain the bounded solution to the Cauchy problem (1.1)–(2.1). Before starting, we recall that the $x$-Fourier transform of $u$ is denoted by $\hat{u}$, the $t$-Laplace one by $\hat{u}$ and the $x$-Fourier/$t$-Laplace one by $\hat{u}$. We refer to Appendix A.3 for conventions and notations related to these transforms.

**Proof of Theorem 2.1.** We start with establishing expression (2.3) for $u = u(t, x)$ which is the simplest to get. To do so, we apply the $t$-Laplace and then the $x$-Fourier transforms to the third equation in (1.1), namely $\partial_t u(t, x) = D\Delta_x u(t, x) + \nu v|_{y=0}(t, x) - \mu u(t, x)$. This yields the algebraic relation

$$-D\|\xi\|^2 \hat{u}(s, \xi) + \nu \hat{v}|_{y=0}(s, \xi) - \mu \hat{u}(s, \xi) = -\hat{u}_0(\xi) + s \hat{u}(s, \xi), \quad s > 0, \xi \in \mathbb{R}^{N-1},$$

which is recast

$$\hat{u}(s, \xi) = \frac{\hat{u}_0(\xi) + \nu \hat{v}|_{y=0}(s, \xi)}{\mu + s + D \|\xi\|^2}, \quad s > 0, \xi \in \mathbb{R}^{N-1}. \quad (4.1)$$

It is now possible to reach (2.3) from (4.1) by using properties and known forms of $t$-Laplace and $x$-Fourier transforms — see Lemmas A.6, (i) and (ii), and A.7, (ii) and (iii):

$$\begin{align*}
\hat{u}(s, \xi) &= \mathcal{L} \left[ t \mapsto e^{-(\mu+D\|\xi\|^2)t} \hat{u}_0(\xi) \right] (s) + \nu \mathcal{L} \left[ t \mapsto e^{-(\mu+D\|\xi\|^2)t} \right] (s) \times \mathcal{L} \left[ t \mapsto \hat{v}|_{y=0}(t, \xi) \right] (s) \\
&= \mathcal{L} \left[ t \mapsto e^{-(\mu+D\|\xi\|^2)t} \hat{u}_0(\xi) \right] (s) + \nu \mathcal{L} \left[ t \mapsto \int_0^t e^{-(\mu+D\|\xi\|^2)(t-\tau)} \hat{v}|_{y=0}(\tau, \xi) d\tau \right] (s) \\
&= \mathcal{LF} \left[ (t, x) \mapsto e^{-\mu t} [u_0 * G_R^D(t, \cdot)](x) \right] (s, \xi) \\
&\quad + \nu \mathcal{LF} \left[ (t, x) \mapsto \int_0^t e^{-\mu(t-\tau)} [v|_{y=0}(\tau, \cdot) * G_R^D(t-\tau, \cdot)](x) d\tau \right] (s, \xi) \\
&= \mathcal{LF} \left[ (t, x) \mapsto e^{-\mu t} \left( U(t, x) + \nu \int_0^t e^{\mu\tau} \int_{\mathbb{R}^{N-1}} G_R^D(t-\tau, x-z) v|_{y=0}(\tau, z) dz d\tau \right) \right] (s, \xi)
\end{align*}$$
which provides (2.3).

We now turn to the expression (2.2) for $v = v(t, X)$. For convenience let us set

$$ \sigma := \sqrt{s + d\|\xi\|^2} \quad \text{and} \quad \Sigma := \sqrt{s + D\|\xi\|^2}. \qquad (4.2) $$

As in Section 3, applying t-Laplace then x-Fourier transforms on the first equation in (1.1), namely $\partial_t v(t, x, y) = d(\Delta_x + \partial_{yy})v(t, x, y)$, leads to the linear second order ODE (where $s > 0$ and $\xi \in \mathbb{R}^{N-1}$ serve as parameters)

$$ d\partial_{yy}\hat{v}(s, \xi, y) - (s + d\|\xi\|^2)\hat{v}(s, \xi, y) = -\hat{v}_0(\xi, y), \quad y > 0, \qquad (4.3) $$

whose solutions are

$$ \hat{v}(s, \xi, y) = e^{\frac{\sigma}{\sqrt{d}} y} \left( C_1 - \frac{1}{2\sqrt{d}\sigma} \int_0^y e^{-\frac{\sigma}{\sqrt{d}}\omega} \hat{v}_0(\xi, \omega) d\omega \right) + e^{-\frac{\sigma}{\sqrt{d}} y} \left( C_2 + \frac{1}{2\sqrt{d}\sigma} \int_0^y e^{\frac{\sigma}{\sqrt{d}}\omega} \hat{v}_0(\xi, \omega) d\omega \right), $$

where $C_1$ and $C_2$ are constants (with respect to $y$) to be determined. In order to insure

$$ \lim_{y \to +\infty} \hat{v}(s, \xi, y) = 0, \text{ one needs} $$

$$ C_1 = \frac{1}{2\sqrt{d}\sigma} \int_0^{+\infty} e^{-\frac{\sigma}{\sqrt{d}}\omega} \hat{v}_0(\xi, \omega) d\omega. \qquad (4.4) $$

Then, calling $A := \frac{\nu}{d}$ and $B := \frac{\mu}{d}$, the exchange condition — second equation in (1.1) — is recast $Av|_{y=0}(t, x) - \partial_y v|_{y=0}(t, x) = Bu(t, x)$; and since

$$ \hat{v}(s, \xi, 0) = C_1 + C_2 \quad \text{and} \quad \partial_y \hat{v}(s, \xi, 0) = \frac{\sigma}{\sqrt{d}}(C_1 - C_2), $$

the exchange condition enforces

$$ \begin{aligned} C_2 &= \left(1 - \frac{2A\sqrt{d}}{\sigma + A\sqrt{d}}\right) C_1 + \frac{B\sqrt{d}}{\sigma + A\sqrt{d}} \hat{u}(s, \xi) \\ &= \left(1 - \frac{2A\sqrt{d}}{\sigma + A\sqrt{d}}\right) C_1 + B\sqrt{d} \frac{\hat{u}_0(\xi) + \nu \hat{v}|_{y=0}(s, \xi)}{(\Sigma^2 + \mu)(\sigma + A\sqrt{d})}, \end{aligned} \qquad (4.5) $$

where we used the expression (4.1) for $\hat{u}(s, \xi)$ and the definition of $\Sigma$ in (4.2). Comparing (3.9), (3.10) with (4.4), (4.5), we see that the “deviation” of $\hat{v}(s, \xi, y)$ (of the field-road model) from the solution to the Heat equation in the half-space (see Section 3) stands in the second term in the right-hand side of (4.5). Hence, using the computations that have been performed in Section 3 and calling $V = V(t, X)$ the solution to (2.4), we reach an *implicit* expression for $\hat{v}(s, \xi, y)$, namely

$$ \hat{v}(s, \xi, y) = \hat{V}(s, \xi, y) + B\sqrt{d} e^{-\frac{\sigma}{\sqrt{d}}y} \frac{\hat{u}_0(\xi) + \nu \hat{v}|_{y=0}(s, \xi)}{(\Sigma^2 + \mu)(\sigma + A\sqrt{d})}. \qquad (4.6) $$
which provides (2.3).

We now turn to the expression (2.2) for $v = v(t, X)$. For convenience let us set

$$ \sigma := \sqrt{s + d\|\xi\|^2} \quad \text{and} \quad \Sigma := \sqrt{s + D\|\xi\|^2}. \qquad (4.2) $$

As in Section 3, applying t-Laplace then x-Fourier transforms on the first equation in (1.1), namely $\partial_t v(t, x, y) = d(\Delta_x + \partial_{yy})v(t, x, y)$, leads to the linear second order ODE (where $s > 0$ and $\xi \in \mathbb{R}^{N-1}$ serve as parameters)

$$ d\partial_{yy}\hat{v}(s, \xi, y) - (s + d\|\xi\|^2)\hat{v}(s, \xi, y) = -\hat{v}_0(\xi, y), \quad y > 0, \qquad (4.3) $$

whose solutions are

$$ \hat{v}(s, \xi, y) = e^{\frac{\sigma}{\sqrt{d}} y} \left( C_1 - \frac{1}{2\sqrt{d}\sigma} \int_0^y e^{-\frac{\sigma}{\sqrt{d}}\omega} \hat{v}_0(\xi, \omega) d\omega \right) + e^{-\frac{\sigma}{\sqrt{d}} y} \left( C_2 + \frac{1}{2\sqrt{d}\sigma} \int_0^y e^{\frac{\sigma}{\sqrt{d}}\omega} \hat{v}_0(\xi, \omega) d\omega \right), $$

where $C_1$ and $C_2$ are constants (with respect to $y$) to be determined. In order to insure $\lim_{y \to +\infty} \hat{v}(s, \xi, y) = 0$, one needs

$$ C_1 = \frac{1}{2\sqrt{d}\sigma} \int_0^{+\infty} e^{-\frac{\sigma}{\sqrt{d}}\omega} \hat{v}_0(\xi, \omega) d\omega. \qquad (4.4) $$

Then, calling $A := \frac{\nu}{d}$ and $B := \frac{\mu}{d}$, the exchange condition — second equation in (1.1) — is recast $Av|_{y=0}(t, x) - \partial_y v|_{y=0}(t, x) = Bu(t, x)$; and since

$$ \hat{v}(s, \xi, 0) = C_1 + C_2 \quad \text{and} \quad \partial_y \hat{v}(s, \xi, 0) = \frac{\sigma}{\sqrt{d}}(C_1 - C_2), $$

the exchange condition enforces

$$ \begin{aligned} C_2 &= \left(1 - \frac{2A\sqrt{d}}{\sigma + A\sqrt{d}}\right) C_1 + \frac{B\sqrt{d}}{\sigma + A\sqrt{d}} \hat{u}(s, \xi) \\ &= \left(1 - \frac{2A\sqrt{d}}{\sigma + A\sqrt{d}}\right) C_1 + B\sqrt{d} \frac{\hat{u}_0(\xi) + \nu \hat{v}|_{y=0}(s, \xi)}{(\Sigma^2 + \mu)(\sigma + A\sqrt{d})}, \end{aligned} \qquad (4.5) $$

where we used the expression (4.1) for $\hat{u}(s, \xi)$ and the definition of $\Sigma$ in (4.2). Comparing (3.9), (3.10) with (4.4), (4.5), we see that the “deviation” of $\hat{v}(s, \xi, y)$ (of the field-road model) from the solution to the Heat equation in the half-space (see Section 3) stands in the second term in the right-hand side of (4.5). Hence, using the computations that have been performed in Section 3 and calling $V = V(t, X)$ the solution to (2.4), we reach an *implicit* expression for $\hat{v}(s, \xi, y)$, namely

$$ \hat{v}(s, \xi, y) = \hat{V}(s, \xi, y) + B\sqrt{d} e^{-\frac{\sigma}{\sqrt{d}}y} \frac{\hat{u}_0(\xi) + \nu \hat{v}|_{y=0}(s, \xi)}{(\Sigma^2 + \mu)(\sigma + A\sqrt{d})}. \qquad (4.6) $$
Evaluating now (4.6) at $y = 0$ yields

$$ \hat{v}|_{y=0}(s, \xi) = \hat{V}|_{y=0}(s, \xi) + B\sqrt{d} \frac{\hat{u}_0(\xi) + \nu \hat{v}|_{y=0}(s, \xi)}{(\Sigma^2 + \mu)(\sigma + A\sqrt{d})}, $$

whence

$$ \hat{v}|_{y=0}(s, \xi) = \frac{B\sqrt{d}}{(\Sigma^2 + \mu)(\sigma + A\sqrt{d}) - \nu B\sqrt{d}} \hat{u}_0(\xi) + \frac{(\Sigma^2 + \mu)(\sigma + A\sqrt{d})}{(\Sigma^2 + \mu)(\sigma + A\sqrt{d}) - \nu B\sqrt{d}} \hat{V}|_{y=0}(s, \xi). \quad (4.7) $$

Plugging (4.7) into (4.6) now yields the explicit expression

$$ \hat{v}(s, \xi, y) = \hat{V}(s, \xi, y) + B\sqrt{d} \frac{e^{-\frac{\sigma}{\sqrt{d}}y}}{(\Sigma^2 + \mu)(\sigma + A\sqrt{d}) - \nu B\sqrt{d}} \left[ \hat{u}_0(\xi) + \nu \hat{V}|_{y=0}(s, \xi) \right]. \quad (4.8) $$

Our task is now to take the inverse Fourier/Laplace transform of the second term in the above right-hand-side. Letting

$$ \delta := (D - d)\|\xi\|^2, $$

one has $\Sigma^2 = \sigma^2 + \delta$. Then define

$$ \hat{\Lambda}(s, \xi, y) := \frac{e^{-\frac{\sigma}{\sqrt{d}}y}}{(\Sigma^2 + \mu)(\sigma + A\sqrt{d}) - \nu B\sqrt{d}} = \frac{e^{-\frac{\sigma}{\sqrt{d}}y}}{\sigma^3 + A\sqrt{d}\sigma^2 + (\mu + \delta)\sigma + A\sqrt{d}\delta}. \quad (4.9) $$

To match up $\hat{\Lambda}(s, \xi, y)$ with known t-Laplace transforms, we need to expand (4.9) into partial fractions, which requires a study of the polynomials

$$ P_{\delta}(\sigma) := \sigma^3 + A\sqrt{d}\sigma^2 + (\mu + \delta)\sigma + A\sqrt{d}\delta, \quad \delta \ge 0 \text{ (since } D > d). $$

We show in Appendix A.2 that the three complex roots $\alpha(\delta)$, $\beta(\delta)$, $\gamma(\delta)$ of $P_\delta$ are simple for almost all $\xi \in \mathbb{R}^{N-1}$. Hence, for these $\xi$,

$$ \hat{\Lambda}(s, \xi, y) = a \frac{e^{-\frac{\sigma}{\sqrt{d}}y}}{\sigma - \alpha} + b \frac{e^{-\frac{\sigma}{\sqrt{d}}y}}{\sigma - \beta} + c \frac{e^{-\frac{\sigma}{\sqrt{d}}y}}{\sigma - \gamma}, \quad (4.10) $$

where $a$, $b$ and $c$ have been defined in (2.9), and do not depend on the $s$ variable. Expression (4.10) of $\hat{\Lambda}$ shows up a known form of $t$-Laplace transform — see Lemma A.6, (v) — evaluated at $\sigma^2 = s + d\|\xi\|^2$, namely

$$ \hat{\Lambda}(s, \xi, y) = \mathcal{L} \left[ t \mapsto e^{-\frac{y^2}{4dt}} \left( \frac{a+b+c}{\sqrt{\pi t}} + a\alpha \frac{\text{Erfc}}{\Gamma} \left( \frac{-2\alpha\sqrt{d}t+y}{2\sqrt{dt}} \right) + b\beta \frac{\text{Erfc}}{\Gamma} \left( \frac{-2\beta\sqrt{d}t+y}{2\sqrt{dt}} \right) \right. \right. \\ \left. \left. + c\gamma \frac{\text{Erfc}}{\Gamma} \left( \frac{-2\gamma\sqrt{d}t+y}{2\sqrt{dt}} \right) \right) \right] (s+d\|\xi\|^2). $$
Evaluating now (4.6) at $y = 0$ yields

$$ \hat{v}_{|y=0}(s, \xi) = \hat{V}_{|y=0}(s, \xi) + B\sqrt{d} \frac{\hat{u}_0(\xi) + \nu \hat{v}_{|y=0}(s, \xi)}{(\Sigma^2 + \mu)(\sigma + A\sqrt{d})}, $$

whence

$$ \hat{v}_{|y=0}(s, \xi) = \frac{B\sqrt{d}}{(\Sigma^2 + \mu)(\sigma + A\sqrt{d}) - \nu B\sqrt{d}} \hat{u}_0(\xi) + \frac{(\Sigma^2 + \mu)(\sigma + A\sqrt{d})}{(\Sigma^2 + \mu)(\sigma + A\sqrt{d}) - \nu B\sqrt{d}} \hat{V}_{|y=0}(s, \xi). \quad (4.7) $$

Plugging (4.7) into (4.6) now yields the explicit expression

$$ \hat{v}(s, \xi, y) = \hat{V}(s, \xi, y) + B\sqrt{d} \frac{e^{-\frac{\sigma}{\sqrt{d}}y}}{(\Sigma^2 + \mu)(\sigma + A\sqrt{d}) - \nu B\sqrt{d}} \left[ \hat{u}_0(\xi) + \nu \hat{V}_{|y=0}(s, \xi) \right]. \quad (4.8) $$

Our task is now to take the inverse Fourier/Laplace transform of the second term in the above right-hand-side. Letting

$$ \delta := (D - d)\|\xi\|^2, $$

one has $\Sigma^2 = \sigma^2 + \delta$. Then define

$$ \hat{\Lambda}(s, \xi, y) := \frac{e^{-\frac{\sigma}{\sqrt{d}}y}}{(\Sigma^2 + \mu)(\sigma + A\sqrt{d}) - \nu B\sqrt{d}} = \frac{e^{-\frac{\sigma}{\sqrt{d}}y}}{\sigma^3 + A\sqrt{d}\sigma^2 + (\mu + \delta)\sigma + A\sqrt{d}\delta}. \quad (4.9) $$

To match up $\hat{\Lambda}(s, \xi, y)$ with known t-Laplace transforms, we need to expand (4.9) into partial fractions, which requires a study of the polynomials

$$ P_{\delta}(\sigma) := \sigma^3 + A\sqrt{d}\sigma^2 + (\mu + \delta)\sigma + A\sqrt{d}\delta, \quad \delta \ge 0 \text{ (since } D > d). $$

We show in Appendix A.2 that the three complex roots $\alpha(\delta)$, $\beta(\delta)$, $\gamma(\delta)$ of $P_{\delta}$ are simple for almost all $\xi \in \mathbb{R}^{N-1}$. Hence, for these $\xi$,

$$ \hat{\Lambda}(s, \xi, y) = a \frac{e^{-\frac{\sigma}{\sqrt{d}}y}}{\sigma - \alpha} + b \frac{e^{-\frac{\sigma}{\sqrt{d}}y}}{\sigma - \beta} + c \frac{e^{-\frac{\sigma}{\sqrt{d}}y}}{\sigma - \gamma}, \quad (4.10) $$

where $a$, $b$ and $c$ have been defined in (2.9), and do not depend on the $s$ variable. Expression (4.10) of $\hat{\Lambda}$ shows up a known form of $t$-Laplace transform — see Lemma A.6, (v) — evaluated at $\sigma^2 = s + d\|\xi\|^2$, namely

$$ \hat{\Lambda}(s, \xi, y) = \mathcal{L} \left[ t \mapsto e^{-\frac{y^2}{4dt}} \left( \frac{a+b+c}{\sqrt{\pi t}} + a\alpha \frac{\text{Erfc}}{\Gamma} \left( \frac{-2\alpha\sqrt{d}t+y}{2\sqrt{dt}} \right) + b\beta \frac{\text{Erfc}}{\Gamma} \left( \frac{-2\beta\sqrt{d}t+y}{2\sqrt{dt}} \right) \right. \right. \\ \left. \left. + c\gamma \frac{\text{Erfc}}{\Gamma} \left( \frac{-2\gamma\sqrt{d}t+y}{2\sqrt{dt}} \right) \right) \right] (s+d\|\xi\|^2). $$
Observing that $a + b + c = 0$, using Lemma A.7, (iv), and recalling the notation $\Phi_{\bullet}$ defined in statement of Theorem 2.1, see (2.10), we thus reach

$$ \hat{\Lambda}(s, \xi, y) = \mathcal{L} \left[ t \mapsto e^{-\frac{y^2}{4dt}} [a\alpha\Phi_{\alpha} + b\beta\Phi_{\beta} + c\gamma\Phi_{\gamma}] (t, \xi, y) e^{-dt\|\xi\|^2} \right] (s). \quad (4.11) $$

We know from Lemma A.6, (i), that $e^{-dt\|\xi\|^2} = \mathcal{F}[G_R^d(t, \cdot)](\xi)$ — notice that here we do mean $G_R^d$, and not $G_R^D$, see (2.6). Denoting

$$ \Phi := a\alpha\Phi_{\alpha} + b\beta\Phi_{\beta} + c\gamma\Phi_{\gamma}, $$

we also know — see Lemma A.7, (i) — that $\Phi(t, \xi, y) = \frac{1}{(2\pi)^{N-1}} \mathcal{F}[\mathcal{F}[\Phi(t, -\cdot, y)]](\xi)$. Hence, from Lemma A.7, (ii), we end up with

$$ \begin{aligned} \hat{\Lambda}(s, \xi, y) &= \mathcal{L}\mathcal{F} \left[ (t, x) \mapsto \frac{e^{-\frac{y^2}{4dt}}}{(2\pi)^{N-1}} \int_{\eta \in \mathbb{R}^{N-1}} \int_{\chi \in \mathbb{R}^{N-1}} \Phi(t, \chi, y) G_R^d(t, x-\eta) e^{i\chi \cdot \eta} d\chi d\eta \right] (s, \xi) \\ &= \mathcal{L}\mathcal{F} \left[ (t, x) \mapsto \frac{e^{-\frac{y^2}{4dt}}}{(2\pi)^{N-1}} \int_{\chi \in \mathbb{R}^{N-1}} \Phi(t, \chi, y) e^{i\chi \cdot x} \int_{\eta \in \mathbb{R}^{N-1}} G_R^d(t, x-\eta) e^{-i\chi \cdot (x-\eta)} d\eta d\chi \right] (s, \xi) \\ &= \mathcal{L}\mathcal{F} \left[ (t, x) \mapsto \frac{e^{-\frac{y^2}{4dt}}}{(2\pi)^{N-1}} \int_{\chi \in \mathbb{R}^{N-1}} \Phi(t, \chi, y) e^{-dt\|\chi\|^2+i\chi \cdot x} d\chi \right] (s, \xi). \end{aligned} \quad (4.12) $$

In other words, and as expected, we do have

$$ \hat{\Lambda}(s, \xi, y) = \mathcal{L}\mathcal{F} \left[ (t, x) \mapsto \Lambda(t, x, y) \right] (s, \xi), $$

where $\Lambda = \Lambda(t, x, y)$ has been introduced in the statement of Theorem 2.1, see (2.7). Finally, our last move consists in using Lemma A.7, (ii) and (iii), in (4.8), leading us, knowing $\Lambda$, to

$$ \hat{v}(s, \xi, y) = \mathcal{L}\mathcal{F} \left[ (t, x) \mapsto V(t, x, y) + B\sqrt{d} \int_{\mathbb{R}^{N-1}} \Lambda(t, z, y) u_0(x-z) dz \right. \\ \left. \qquad + \nu B\sqrt{d} \int_0^t \int_{\mathbb{R}^{N-1}} \Lambda(\tau, z, y) V|_{y=0}(t-\tau, x-z) dz d\tau \right] (s, \xi) $$

which provides (2.2).

<figure><img src="image_1.png" /></figure>

**Remark 4.1** (When $D \le d$) Let us first discuss the case $D < d$. Then, when $\xi$ browses $\mathbb{R}^{N-1}$, $\delta = (D-d)\|\xi\|^2$ browses $(-\infty, 0]$. As in Appendix A.2, it can be shown that, for almost all $\delta \le 0$ (and thus for almost all $\xi \in \mathbb{R}^{N-1}$), the roots of $P_{\delta}$ are simple. We thus can write again (4.10) and the above proof readily applies.

We now turn on the critical case $D = d$. When $\xi$ browses $\mathbb{R}^{N-1}$, $\delta = (D-d)\|\xi\|^2$ remains stuck at zero — so that $\Phi(t, \xi, y)$ is actually independent on $\xi$ — and only the polynomial $P_0$ plays a role. If $\mu \ne \frac{\nu^2}{4d}$, then $P_0$ has three simple roots (see Figure X and Appendix A.2) and the above proof again readily applies. On the other hand, if $\mu = \frac{\nu^2}{4d}$, then $P_0$ has a simple root and a double root. In this case, one should go back to (4.9) and, rather than (4.10), use the adequate expansion into partial fractions. Details are omitted.
Observing that $a + b + c = 0$, using Lemma A.7, (iv), and recalling the notation $\Phi_{\bullet}$ defined in statement of Theorem 2.1, see (2.10), we thus reach

$$ \hat{\Lambda}(s, \xi, y) = \mathcal{L} \left[ t \mapsto e^{-\frac{y^2}{4dt}} \left[ a\alpha\Phi_{\alpha} + b\beta\Phi_{\beta} + c\gamma\Phi_{\gamma} \right] (t, \xi, y) e^{-dt\|\xi\|^2} \right] (s). \quad (4.11) $$

We know from Lemma A.6, (i), that $e^{-dt\|\xi\|^2} = \mathcal{F}[G_R^d(t, \cdot)](\xi)$ — notice that here we do mean $G_R^d$, and not $G_R^D$, see (2.6). Denoting

$$ \Phi := a\alpha\Phi_{\alpha} + b\beta\Phi_{\beta} + c\gamma\Phi_{\gamma}, $$

we also know — see Lemma A.7, (i) — that $\Phi(t, \xi, y) = \frac{1}{(2\pi)^{N-1}} \mathcal{F}[\mathcal{F}[\Phi(t, -\cdot, y)]](\xi)$. Hence, from Lemma A.7, (ii), we end up with

$$ \begin{aligned} \hat{\Lambda}(s, \xi, y) &= \mathcal{L}\mathcal{F} \left[ (t, x) \mapsto \frac{e^{-\frac{y^2}{4dt}}}{(2\pi)^{N-1}} \int_{\eta \in \mathbb{R}^{N-1}} \int_{\chi \in \mathbb{R}^{N-1}} \Phi(t, \chi, y) G_R^d(t, x-\eta) e^{i\chi \cdot \eta} d\chi d\eta \right] (s, \xi) \\ &= \mathcal{L}\mathcal{F} \left[ (t, x) \mapsto \frac{e^{-\frac{y^2}{4dt}}}{(2\pi)^{N-1}} \int_{\chi \in \mathbb{R}^{N-1}} \Phi(t, \chi, y) e^{i\chi \cdot x} \int_{\eta \in \mathbb{R}^{N-1}} G_R^d(t, x-\eta) e^{-i\chi \cdot (x-\eta)} d\eta d\chi \right] (s, \xi) \\ &= \mathcal{L}\mathcal{F} \left[ (t, x) \mapsto \frac{e^{-\frac{y^2}{4dt}}}{(2\pi)^{N-1}} \int_{\chi \in \mathbb{R}^{N-1}} \Phi(t, \chi, y) e^{-dt\|\chi\|^2+i\chi \cdot x} d\chi \right] (s, \xi). \end{aligned} \quad (4.12) $$

In other words, and as expected, we do have

$$ \hat{\Lambda}(s, \xi, y) = \mathcal{L}\mathcal{F} \left[ (t, x) \mapsto \Lambda(t, x, y) \right] (s, \xi), $$

where $\Lambda = \Lambda(t, x, y)$ has been introduced in the statement of Theorem 2.1, see (2.7). Finally, our last move consists in using Lemma A.7, (ii) and (iii), in (4.8), leading us, knowing $\Lambda$, to

$$ \hat{v}(s, \xi, y) = \mathcal{L}\mathcal{F} \left[ (t, x) \mapsto V(t, x, y) + B\sqrt{d} \int_{\mathbb{R}^{N-1}} \Lambda(t, z, y) u_0(x-z) dz \right. \\ \left. \qquad + \nu B\sqrt{d} \int_0^t \int_{\mathbb{R}^{N-1}} \Lambda(\tau, z, y) V|_{y=0}(t-\tau, x-z) dz d\tau \right] (s, \xi) $$

which provides (2.2).

<figure><img src="image_3.png" /></figure>

**Remark 4.1** (When $D \le d$) Let us first discuss the case $D < d$. Then, when $\xi$ browses $\mathbb{R}^{N-1}$, $\delta = (D-d)\|\xi\|^2$ browses $(-\infty, 0]$. As in Appendix A.2, it can be shown that, for almost all $\delta \le 0$ (and thus for almost all $\xi \in \mathbb{R}^{N-1}$), the roots of $P_{\delta}$ are simple. We thus can write again (4.10) and the above proof readily applies.

We now turn on the critical case $D = d$. When $\xi$ browses $\mathbb{R}^{N-1}$, $\delta = (D-d)\|\xi\|^2$ remains stuck at zero — so that $\Phi(t, \xi, y)$ is actually independent on $\xi$ — and only the polynomial $P_0$ plays a role. If $\mu \ne \frac{\nu^2}{4d}$, then $P_0$ has three simple roots (see Figure X and Appendix A.2) and the above proof again readily applies. On the other hand, if $\mu = \frac{\nu^2}{4d}$, then $P_0$ has a simple root and a double root. In this case, one should go back to (4.9) and, rather than (4.10), use the adequate expansion into partial fractions. Details are omitted.
# 5 The decay rate of the field-road model

In this section, we estimate the decay rate of the $L^\infty$ norm of the solution to the Cauchy problem (1.1)–(2.11). Owing to the parabolic comparison principle, we can assume without loss of generality that $v_0$ and $u_0$ are smooth.

*Proof of Theorem 2.2.* We first handle the case of $v$ — from which the one of $u$ shall easily ensue. Reminding the expression (2.2) of $v$,

$$v(t, X) = V(t, X) + \frac{\mu}{\sqrt{d}} \underbrace{\int_{\mathbb{R}^{N-1}} \Lambda(t, z, y) u_0(x-z) dz}_{\text{Let us call this } P(t, x, y),} \\ + \frac{\mu \nu}{\sqrt{d}} \underbrace{\int_0^t \int_{\mathbb{R}^{N-1}} \Lambda(s, z, y) V|_{y=0}(t-s, x-z) dz ds}_{\text{and that } Q(t, x, y)}. \quad (5.1)$$

The first term, namely $V(t, X)$, is the solution to the Cauchy problem (2.4) whose $L^\infty$ control is given in statement of Theorem 3.1 (iii):

$$\|V(t, \cdot)\|_{L^{\infty}(\mathbb{R}_{+}^{N})} \lesssim \frac{\int_{\mathbb{R}_{+}^{N}} (1+\omega)|v_0(Z)| dZ}{t^{\frac{N+1}{2}}}, \quad \forall t > 0. \quad (5.2)$$

To deal with terms $P$ and $Q$, we will need the following lemma, whose proof is postponed at the end of this section.

**Lemma 5.1** ($L^\infty$ control of $\Phi$) *Recall the notation introduced in the proof of Theorem 2.1:*

$$\Phi(t, \xi, y) := [a\alpha\Phi_{\alpha} + b\beta\Phi_{\beta} + c\gamma\Phi_{\gamma}](t, \xi, y). \quad (5.3)$$

Then

$$S(t) := \sup_{y \ge 0} \sup_{\xi \in \mathbb{R}^{N-1}} |\Phi(t, \xi, y)| \lesssim \frac{1}{\sqrt{1+t}}, \quad \forall t > 0. \quad (5.4)$$

From here, the estimate of $P$ is straight. Indeed, using expression (2.7) of $\Lambda$ as written in Theorem 2.1, we get

$$\begin{aligned} |P(t, x, y)| &= \left| \int_{z \in \mathbb{R}^{N-1}} \frac{e^{-\frac{y^2}{4dt}}}{(2\pi)^{N-1}} \int_{\xi \in \mathbb{R}^{N-1}} \Phi(t, \xi, y) e^{-dt\|\xi\|^2 + i\xi \cdot z} u_0(x-z) d\xi dz \right| \\ &\le \frac{S(t)}{(2\pi)^{N-1}} \int_{\xi \in \mathbb{R}^{N-1}} e^{-dt\|\xi\|^2} d\xi \int_{z \in \mathbb{R}^{N-1}} |u_0(z)| dz \\ &= \frac{S(t)}{(4\pi dt)^{\frac{N-1}{2}}} \int_{z \in \mathbb{R}^{N-1}} |u_0(z)| dz. \end{aligned}$$

Whence, thanks to Lemma 5.1,

$$\|P(t, \cdot)\|_{L^{\infty}(\mathbb{R}_{+}^{N})} \lesssim \frac{\int_{\mathbb{R}^{N-1}} |u_0(z)| dz}{t^{\frac{N}{2}}}, \quad \forall t > 0. \quad (5.5)$$
**Remark 4.2** (Adding linear growth) One may also want to consider the case with linear growth both in the field and on the road. If, in the purely diffusive system (1.1), $+pv$ and $+qu$ are added to the $v$-equation and the $u$-equation respectively, then by considering $(\tilde{v}(t, x), \tilde{u}(t, x)) := e^{-qt}(v(t, x), u(t, x))$, the system is recast (after dropping the tildes and defining $r := p - q$)

$$ egin{cases} 
\partial_t v = d\Delta v + rv, & t > 0, x \in \mathbb{R}^{N-1}, y > 0, \\ 
-d\partial_y v|_{y=0} = \mu u - \nu v|_{y=0}, & t > 0, x \in \mathbb{R}^{N-1}, \\ 
\partial_t u = D\Delta u + \nu v|_{y=0} - \mu u, & t > 0, x \in \mathbb{R}^{N-1}, 
\end{cases} $$

which is the linearized system around the null steady state of the Fisher-KPP system (1.2) originally introduced in [10, 9]. Then the first change in the proof is that (4.3) is transferred into

$$ d\partial_{yy}\hat{v}(s, \xi, y) - \left(s + d\|\xi\|^2 - r\right)\hat{v}(s, \xi, y) = -\hat{v}_0(\xi, y), \quad y > 0. \quad (4.13) $$

Since individuals are expected to better grow in the field than on the road, we would typically have $r > 0$ and, because of that, the structure of the solutions to the linear second order ODE (4.13) does depend upon parameters $s$ and $r$, and different cases should be considered. As a result our method readily applies but, to reach a rather explicit expression of the solution possibly yielding new insights, further heavy computations would be necessary.

# 5 The decay rate of the field-road model

In this section, we estimate the decay rate of the $L^\infty$ norm of the solution to the Cauchy problem (1.1)–(2.11). Owing to the parabolic comparison principle, we can assume without loss of generality that $v_0$ and $u_0$ are smooth.

*Proof of Theorem 2.2.* We first handle the case of $v$ — from which the one of $u$ shall easily ensue. Reminding the expression (2.2) of $v$,

$$ v(t, X) = V(t, X) + \frac{\mu}{\sqrt{d}} \underbrace{\int_{\mathbb{R}^{N-1}} \Lambda(t, z, y) u_0(x-z) dz}_{\text{Let us call this } P(t, x, y),} + \frac{\mu \nu}{\sqrt{d}} \underbrace{\int_0^t \int_{\mathbb{R}^{N-1}} \Lambda(s, z, y) V|_{y=0}(t-s, x-z) dz ds}_{\text{and that } Q(t, x, y)}. \quad (5.1) $$

The first term, namely $V(t, X)$, is the solution to the Cauchy problem (2.4) whose $L^\infty$ control is given in statement of Theorem 3.1 (iii):

$$ \|V(t, \cdot)\|_{L^\infty(\mathbb{R}_+^N)} \lesssim \frac{\int_{\mathbb{R}_+^N} (1+\omega)|v_0(Z)| dZ}{t^{\frac{N+1}{2}}}, \quad \forall t > 0. \quad (5.2) $$

To deal with terms $P$ and $Q$, we will need the following lemma, whose proof is postponed at the end of this section.
Controlling $Q$ is more refined. Observe from Theorem 3.1 and Corollary 3.2 that

$$V(t, x, y) = \int_{\eta \in \mathbb{R}^{N-1}} \int_{\omega=0}^{\infty} G_{R}^{d}(t, x - \eta) H_{\theta}^{(1)}(t, y, \omega) v_{0}(\eta, \omega) \, d\omega \, d\eta,$$

where $\theta := \frac{\nu}{1+\nu} \in (0, 1)$, so that

$$ 
\begin{aligned}
Q(t, x, y) &= \int_{s=0}^{t} \int_{z \in \mathbb{R}^{N-1}} \frac{e^{-\frac{y^2}{4ds}}}{(2\pi)^{N-1}} \int_{\xi \in \mathbb{R}^{N-1}} \Phi(s, \xi, y) e^{-ds\|\xi\|^2 + i\xi \cdot z} \\
& \qquad \int_{\eta \in \mathbb{R}^{N-1}} \int_{\omega=0}^{\infty} G_{R}^{d}(t-s, x-z-\eta) H_{\theta}^{(1)}(t-s, 0, \omega) v_{0}(\eta, \omega) \, d\omega \, d\eta \, d\xi \, dz \, ds \\
&= \int_{s=0}^{t} \frac{e^{-\frac{y^2}{4ds}}}{(2\pi)^{N-1}} \int_{\xi \in \mathbb{R}^{N-1}} \Phi(s, \xi, y) e^{-ds\|\xi\|^2} \int_{\eta \in \mathbb{R}^{N-1}} \int_{\omega=0}^{\infty} H_{\theta}^{(1)}(t-s, 0, \omega) v_{0}(\eta, \omega) \\
& \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \int_{z \in \mathbb{R}^{N-1}} G_{R}^{d}(t-s, x-z-\eta) e^{i\xi \cdot z} dz \, d\omega \, d\eta \, d\xi \, ds.
\end{aligned}
$$

The integral in $z$ is nothing else than $\mathcal{F}[G_R^d(t-s, \cdot + x - \eta)](\xi) = e^{-d(t-s)\|\xi\|^2 + i\xi \cdot (x-\eta)}$ so that

$$ 
\begin{aligned}
Q(t, x, y) &= \int_{s=0}^{t} \frac{e^{-\frac{y^2}{4ds}}}{(2\pi)^{N-1}} \int_{\xi \in \mathbb{R}^{N-1}} \Phi(s, \xi, y) e^{-dt\|\xi\|^2} \int_{\eta \in \mathbb{R}^{N-1}} \int_{\omega=0}^{\infty} H_{\theta}^{(1)}(t-s, 0, \omega) v_{0}(\eta, \omega) \\
& \qquad e^{i\xi \cdot (x-\eta)} \, d\omega \, d\eta \, d\xi \, ds \\
&= \int_{s=0}^{t} \frac{e^{-\frac{y^2}{4ds}}}{(2\pi)^{N-1}} \int_{\xi \in \mathbb{R}^{N-1}} \Phi(s, \xi, y) e^{-dt\|\xi\|^2} e^{i\xi \cdot x} \int_{\omega=0}^{\infty} H_{\theta}^{(1)}(t-s, 0, \omega) \hat{v}_{0}(\xi, \omega) \, d\omega \, d\xi \, ds.
\end{aligned}
$$

As a result,

$$|Q(t, x, y)| \le \frac{1}{(2\pi)^{N-1}} \int_{s=0}^{t} \mathcal{S}(s) \int_{\xi \in \mathbb{R}^{N-1}} e^{-dt\|\xi\|^2} \int_{\omega=0}^{\infty} H_{\theta}^{(1)}(t-s, 0, \omega) |\hat{v}_{0}(\xi, \omega)| \, d\omega \, d\xi \, ds.$$

From the control on $\mathcal{S}$ in Lemma 5.1 and estimate (3.7) in Theorem 3.1 (with $N=1$), we reach

$$|Q(t, x, y)| \lesssim \tilde{C}_{v_0} \int_{s=0}^{t} \int_{\xi \in \mathbb{R}^{N-1}} e^{-dt\|\xi\|^2} \frac{1}{\sqrt{1+s}(1+t-s)} d\xi \, ds \lesssim \frac{\tilde{C}_{v_0}}{t^{\frac{N-1}{2}}} \int_{0}^{t} \frac{1}{\sqrt{1+s}(1+t-s)} ds,$$

for some $\tilde{C}_{v_0} \ge 0$. Finally, an elementary computation shows

$$\int_{0}^{t} \frac{1}{\sqrt{1+s}(1+t-s)} ds = \frac{1}{\sqrt{2+t}} \left[ \ln \left( \frac{1+\sqrt{\frac{1+t}{2+t}}}{1-\sqrt{\frac{1+t}{2+t}}} \right) - \ln \left( \frac{1+\sqrt{\frac{1}{2+t}}}{1-\sqrt{\frac{1}{2+t}}} \right) \right] \le \frac{\ln(4t+6)}{\sqrt{2+t}},$$

from which there comes

$$\|Q(t, \cdot)\|_{L^{\infty}(\mathbb{R}_{+}^{N})} \lesssim \tilde{C}_{v_0} \frac{\ln(1+t)}{t^{\frac{N}{2}}}, \quad \forall t > 0. \qquad (5.6)$$

Gathering (5.2), (5.5) and (5.6) brings us to the announced control (2.12) on $v$.
**Lemma 5.1** ($L^\infty$ control of $\Phi$) Recall the notation introduced in the proof of Theorem 2.1:

$$ \Phi(t, \xi, y) := [a\alpha\Phi_{\alpha} + b\beta\Phi_{\beta} + c\gamma\Phi_{\gamma}](t, \xi, y). \quad (5.3) $$

Then

$$ S(t) := \sup_{y \ge 0} \sup_{\xi \in \mathbb{R}^{N-1}} |\Phi(t, \xi, y)| \lesssim \frac{1}{\sqrt{1+t}}, \quad \forall t > 0. \quad (5.4) $$

From here, the estimate of $P$ is straight. Indeed, using expression (2.7) of $\Lambda$ as written in Theorem 2.1, we get

$$ \begin{aligned} |P(t, x, y)| &= \left| \int_{z \in \mathbb{R}^{N-1}} \frac{e^{-\frac{y^2}{4dt}}}{(2\pi)^{N-1}} \int_{\xi \in \mathbb{R}^{N-1}} \Phi(t, \xi, y) e^{-dt\|\xi\|^2 + i\xi \cdot z} u_0(x-z) \, d\xi \, dz \right| \\ &\le \frac{S(t)}{(2\pi)^{N-1}} \int_{\xi \in \mathbb{R}^{N-1}} e^{-dt\|\xi\|^2} \, d\xi \int_{z \in \mathbb{R}^{N-1}} |u_0(z)| \, dz \\ &= \frac{S(t)}{(4\pi dt)^{\frac{N-1}{2}}} \int_{z \in \mathbb{R}^{N-1}} |u_0(z)| \, dz. \end{aligned} $$

Whence, thanks to Lemma 5.1,

$$ \|P(t, \cdot)\|_{L^{\infty}(\mathbb{R}_{+}^{N})} \lesssim \frac{\int_{\mathbb{R}^{N-1}} |u_0(z)| \, dz}{t^{\frac{N}{2}}}, \quad \forall t > 0. \quad (5.5) $$

Controlling $Q$ is more refined. Observe from Theorem 3.1 and Corollary 3.2 that

$$ V(t, x, y) = \int_{\eta \in \mathbb{R}^{N-1}} \int_{\omega=0}^{\infty} G_{R}^{d}(t, x - \eta) H_{\theta}^{(1)}(t, y, \omega) v_{0}(\eta, \omega) \, d\omega \, d\eta, $$

where $\theta := \frac{\nu}{1+\nu} \in (0, 1)$, so that

$$ \begin{aligned} Q(t, x, y) &= \int_{s=0}^{t} \int_{z \in \mathbb{R}^{N-1}} \frac{e^{-\frac{y^2}{4ds}}}{(2\pi)^{N-1}} \int_{\xi \in \mathbb{R}^{N-1}} \Phi(s, \xi, y) e^{-ds\|\xi\|^2 + i\xi \cdot z} \\ & \qquad \int_{\eta \in \mathbb{R}^{N-1}} \int_{\omega=0}^{\infty} G_{R}^{d}(t-s, x-z-\eta) H_{\theta}^{(1)}(t-s, 0, \omega) v_{0}(\eta, \omega) \, d\omega \, d\eta \, d\xi \, dz \, ds \\ &= \int_{s=0}^{t} \frac{e^{-\frac{y^2}{4ds}}}{(2\pi)^{N-1}} \int_{\xi \in \mathbb{R}^{N-1}} \Phi(s, \xi, y) e^{-ds\|\xi\|^2} \int_{\eta \in \mathbb{R}^{N-1}} \int_{\omega=0}^{\infty} H_{\theta}^{(1)}(t-s, 0, \omega) v_{0}(\eta, \omega) \\ & \qquad \int_{z \in \mathbb{R}^{N-1}} G_{R}^{d}(t-s, x-z-\eta) e^{i\xi \cdot z} dz \, d\omega \, d\eta \, d\xi \, ds. \end{aligned} $$

The integral in $z$ is nothing else than $\mathcal{F}[G_R^d(t-s, \cdot + x - \eta)](\xi) = e^{-d(t-s)\|\xi\|^2 + i\xi \cdot (x-\eta)}$ so that

$$ \begin{aligned} Q(t, x, y) &= \int_{s=0}^{t} \frac{e^{-\frac{y^2}{4ds}}}{(2\pi)^{N-1}} \int_{\xi \in \mathbb{R}^{N-1}} \Phi(s, \xi, y) e^{-dt\|\xi\|^2} \int_{\eta \in \mathbb{R}^{N-1}} \int_{\omega=0}^{\infty} H_{\theta}^{(1)}(t-s, 0, \omega) v_{0}(\eta, \omega) \\ & \qquad e^{i\xi \cdot (x-\eta)} \, d\omega \, d\eta \, d\xi \, ds \\ &= \int_{s=0}^{t} \frac{e^{-\frac{y^2}{4ds}}}{(2\pi)^{N-1}} \int_{\xi \in \mathbb{R}^{N-1}} \Phi(s, \xi, y) e^{-dt\|\xi\|^2} e^{i\xi \cdot x} \int_{\omega=0}^{\infty} H_{\theta}^{(1)}(t-s, 0, \omega) \hat{v}_{0}(\xi, \omega) \, d\omega \, d\xi \, ds. \end{aligned} $$
Let us now turn to $u$. We recall its expression as written in (2.3):

$$u(t, x) = e^{-\mu t} U(t, x) + \nu \int_{0}^{t} e^{-\mu(t-s)} \int_{\mathbb{R}^{N-1}} G_{R}^{D}(t-s, x-z) v|_{y=0}(s, z) dz ds,$$

where $U = U(t, x)$ denotes the well-known solution to the Cauchy problem (2.5). Hence, thanks the control (2.12) we just established on $v$,

$$|u(t, x)| \lesssim \frac{\|u_0\|_{L^1(\mathbb{R}^{N-1})} e^{-\mu t}}{t^{\frac{N-1}{2}}} + C_{v_0} e^{-\mu t} \int_0^t e^{\mu s} \frac{\ln(1+s)}{(1+s)^{\frac{N}{2}}} ds + C_{v_0}^{u_0} e^{-\mu t} \int_0^t e^{\mu s} \frac{1}{(1+s)^{\frac{N}{2}}} ds. \quad (5.7)$$

We now treat the two temporal convolutions in (5.7). To do so, remark that

$$\frac{d}{ds} \left( \frac{e^{\mu s} \ln(1+s)}{\mu (1+s)^{\frac{N}{2}}} \right) = e^{\mu s} \frac{\ln(1+s)}{(1+s)^{\frac{N}{2}}} + \frac{e^{\mu s} (1 - \frac{N}{2} \ln(1+s))}{\mu (1+s)^{\frac{N+2}{2}}} \stackrel{s \to \infty}{\sim} e^{\mu s} \frac{\ln(1+s)}{(1+s)^{\frac{N}{2}}}$$

yields $\int_0^t e^{\mu s} \frac{\ln(1+s)}{(1+s)^{\frac{N}{2}}} ds \stackrel{t \to \infty}{\sim} \frac{e^{\mu t}}{\mu} \frac{\ln(1+t)}{(1+t)^{\frac{N}{2}}}$. Similarly,

$$\frac{d}{ds} \left( \frac{e^{\mu s}}{\mu} \frac{1}{(1+s)^{\frac{N}{2}}} \right) = e^{\mu s} \frac{1}{(1+s)^{\frac{N}{2}}} - \frac{e^{\mu s}}{\mu} \frac{\frac{N}{2}}{(1+s)^{\frac{N+2}{2}}} \stackrel{s \to \infty}{\sim} e^{\mu s} \frac{1}{(1+s)^{\frac{N}{2}}}$$

yields $\int_0^t e^{\mu s} \frac{1}{(1+s)^{\frac{N}{2}}} ds \stackrel{t \to \infty}{\sim} \frac{e^{\mu t}}{\mu} \frac{1}{(1+t)^{\frac{N}{2}}}$. Therefore,

$$|u(t, x)| \lesssim \frac{\|u_0\|_{L^1(\mathbb{R}^{N-1})} e^{-\mu t}}{t^{\frac{N-1}{2}}} + \frac{C_{v_0} \ln(1+t)}{(1+t)^{\frac{N}{2}}} + \frac{C_{v_0}^{u_0}}{(1+t)^{\frac{N}{2}}},$$

which provides the control (2.13) on $u$ up to changing the value of $C_{v_0}^{u_0}$ if necessary. $\square$

It now remains to prove Lemma 5.1.

**Proof of Lemma 5.1.** The proof is related to the study of the $\delta$-indexed polynomials $P_\delta$ in Appendix A.2. As $\delta = (D-d)\|\xi\|^2$ browses $\mathbb{R}_+$, the polynomials $P_\delta$ may provide double or triple roots for, at most, two values of $\delta$ — see Lemmas A.2 to A.5 (i). In Lemmas A.2 to A.5 we build a closed set $\mathcal{E}$ of $\mathbb{R}_+$ so that, for $\delta \in \mathcal{E}$, the roots $\alpha, \beta, \gamma$ remain away from each other and so that each connected component of $\mathcal{E}^c$ contains one and only one $\delta_i$ ($i = 0, 1$ or 2) which provides a multiple root $\lambda_i$. This proof is divided into three steps:

(i) We treat the cases where $\delta \in \mathcal{E}$ for which, thanks to the boundedness of $a, b$ and $c$, we can consider $a\alpha\Phi_\alpha, b\beta\Phi_\beta$ and $c\gamma\Phi_\gamma$ independently.

We deal then with the situations where $\delta \in \mathcal{E}^c$. In those, because $a, b$ and $c$ may be unbounded we must consider $a\alpha\Phi_\alpha, b\beta\Phi_\beta$ and $c\gamma\Phi_\gamma$ together. Let $\delta \in \mathcal{E}^c$, then there is $i \in \{0, 1, 2\}$ such that $\delta \in (\delta_i - \eta, \delta_i + \eta)$ — recall that $\delta = \delta_i$ provides the multiple root $\lambda_i$ and that $\eta$ is defined in Lemmas A.3 to A.5.

(ii) We treat the cases where the root $\lambda_i$ is double.

(iii) We treat the cases where the root $\lambda_i$ is triple.
As a result,

$$|Q(t, x, y)| \le \frac{1}{(2\pi)^{N-1}} \int_{s=0}^{t} \mathcal{S}(s) \int_{\xi \in \mathbb{R}^{N-1}} e^{-dt\|\xi\|^2} \int_{\omega=0}^{\infty} H_{\theta}^{(1)}(t-s, 0, \omega) |\hat{v}_0(\xi, \omega)| d\omega d\xi ds.$$

From the control on $\mathcal{S}$ in Lemma 5.1 and estimate (3.7) in Theorem 3.1 (with $N = 1$), we reach

$$|Q(t, x, y)| \lesssim \tilde{C}_{v_0} \int_{s=0}^{t} \int_{\xi \in \mathbb{R}^{N-1}} e^{-dt\|\xi\|^2} \frac{1}{\sqrt{1+s}(1+t-s)} d\xi ds \lesssim \frac{\tilde{C}_{v_0}}{t^{\frac{N-1}{2}}} \int_{0}^{t} \frac{1}{\sqrt{1+s}(1+t-s)} ds,$$

for some $\tilde{C}_{v_0} \ge 0$. Finally, an elementary computation shows

$$\int_{0}^{t} \frac{1}{\sqrt{1+s}(1+t-s)} ds = \frac{1}{\sqrt{2}+t} \left[ \ln \left( \frac{1+\sqrt{\frac{1+t}{2+t}}}{1-\sqrt{\frac{1+t}{2+t}}} \right) - \ln \left( \frac{1+\sqrt{\frac{1}{2+t}}}{1-\sqrt{\frac{1}{2+t}}} \right) \right] \le \frac{\ln(4t+6)}{\sqrt{2}+t},$$

from which there comes

$$\|Q(t, \cdot)\|_{L^{\infty}(\mathbb{R}_{+}^{N})} \lesssim \tilde{C}_{v_0} \frac{\ln(1+t)}{t^{\frac{N}{2}}}, \quad \forall t > 0. \qquad (5.6)$$

Gathering (5.2), (5.5) and (5.6) brings us to the announced control (2.12) on $v$.

Let us now turn to $u$. We recall its expression as written in (2.3):

$$u(t, x) = e^{-\mu t} U(t, x) + \nu \int_{0}^{t} e^{-\mu(t-s)} \int_{\mathbb{R}^{N-1}} G_{R}^{D}(t-s, x-z) v|_{y=0}(s, z) dz ds,$$

where $U = U(t, x)$ denotes the well-known solution to the Cauchy problem (2.5). Hence, thanks the control (2.12) we just established on $v$,

$$|u(t, x)| \lesssim \frac{\|u_0\|_{L^1(\mathbb{R}^{N-1})} e^{-\mu t}}{t^{\frac{N-1}{2}}} + C_{v_0} e^{-\mu t} \int_0^t e^{\mu s} \frac{\ln(1+s)}{(1+s)^{\frac{N}{2}}} ds + C_{v_0} e^{-\mu t} \int_0^t e^{\mu s} \frac{1}{(1+s)^{\frac{N}{2}}} ds. \qquad (5.7)$$

We now treat the two temporal convolutions in (5.7). To do so, remark that

$$\frac{d}{ds} \left( \frac{e^{\mu s} \ln(1+s)}{\mu (1+s)^{\frac{N}{2}}} \right) = e^{\mu s} \frac{\ln(1+s)}{(1+s)^{\frac{N}{2}}} + \frac{e^{\mu s}}{\mu} \frac{1 - \frac{N}{2} \ln(1+s)}{(1+s)^{\frac{N+2}{2}}} \stackrel{s \to \infty}{\sim} e^{\mu s} \frac{\ln(1+s)}{(1+s)^{\frac{N}{2}}}$$

yields $\int_0^t e^{\mu s} \frac{\ln(1+s)}{(1+s)^{\frac{N}{2}}} ds \stackrel{t \to \infty}{\sim} \frac{e^{\mu t}}{\mu} \frac{\ln(1+t)}{(1+t)^{\frac{N}{2}}}$. Similarly,

$$\frac{d}{ds} \left( \frac{e^{\mu s}}{\mu (1+s)^{\frac{N}{2}}} \right) = e^{\mu s} \frac{1}{(1+s)^{\frac{N}{2}}} - \frac{e^{\mu s}}{\mu} \frac{\frac{N}{2}}{(1+s)^{\frac{N+2}{2}}} \stackrel{s \to \infty}{\sim} e^{\mu s} \frac{1}{(1+s)^{\frac{N}{2}}}$$

yields $\int_0^t e^{\mu s} \frac{1}{(1+s)^{\frac{N}{2}}} ds \stackrel{t \to \infty}{\sim} \frac{e^{\mu t}}{\mu} \frac{1}{(1+t)^{\frac{N}{2}}}$. Therefore,

$$|u(t, x)| \lesssim \frac{\|u_0\|_{L^1(\mathbb{R}^{N-1})} e^{-\mu t}}{t^{\frac{N-1}{2}}} + \frac{C_{v_0} \ln(1+t)}{(1+t)^{\frac{N}{2}}} + \frac{C_{v_0}^{\mu_0}}{(1+t)^{\frac{N}{2}}},$$

which provides the control (2.13) on $u$ up to changing the value of $C_{v_0}^{\mu_0}$ if necessary. $\square$
Before starting, we emphasize that (A.5) insures $\mathfrak{Re}(\alpha)$, $\mathfrak{Re}(\beta)$, $\mathfrak{Re}(\gamma) \le 0$ which allows us to use the estimates (A.4) at point $z = \frac{-2\cdot\sqrt{dt+y}}{2\sqrt{dt}}$ for $\cdot$ in the convex hull of $\alpha$, $\beta$ and $\gamma$. Notice that, in the sequel, we shall often write $\Phi_{\cdot}$ for $\Phi_{\cdot}(t, \xi, y)$.

Step (i). Because $\inf_{\delta \in \mathcal{E}} (|\alpha - \beta|, |\alpha - \gamma|, |\beta - \gamma|) > \varepsilon$, we have

$$ \sup_{\delta \in \mathcal{E}} (|a|, |b|, |c|) < \frac{1}{\varepsilon^2}. \qquad (5.8) $$

Next, using (A.4), we have, for $\lambda \in \{\alpha, \beta, \gamma\}$,

$$ |\lambda\Phi_{\lambda}| = \left| \lambda \frac{\text{Erfc}}{\Gamma} \left( \frac{-2\lambda\sqrt{dt} + y}{2\sqrt{dt}} \right) \right| \lesssim \frac{2|\lambda|\sqrt{dt}}{-2\lambda\sqrt{dt} + y} \le \frac{1}{\sqrt{t}}, $$

the last inequality coming from $\mathfrak{Re}(\lambda) \le 0$ and $y \ge 0$. Gathering the latter inequality with (5.8) provides, in view of (5.3),

$$ \sup_{y \ge 0} \sup_{\delta \in \mathcal{E}} |\Phi(t, \xi, y)| \lesssim \frac{1}{\sqrt{1+t}} \quad \text{for all } t \ge 1. \qquad (5.9) $$

We now turn to the case where $0 < t \le 1$ (still when $\delta \in \mathcal{E}$). Let $\delta_{\infty} > 0$ be defined as in Lemma A.1. It is readily seen that

$$ \sup_{t \in (0,1]} \sup_{y \ge 0} \sup_{\delta \in \mathcal{E} \cap [0, \delta_{\infty}]} |\Phi(t, \xi, y)| < +\infty. \qquad (5.10) $$

Next let us consider $\delta \ge \delta_{\infty}$. Since $\gamma = \bar{\beta}$, $c = \bar{b}$, $\Phi_{\gamma} = \bar{\Phi}_{\beta}$, we deduce from (5.3) that

$$ |\Phi| \le |a\alpha\Phi_{\alpha}| + 2|b\beta\Phi_{\beta}|. \qquad (5.11) $$

The first term in the above right-hand-side is clearly bounded while, for the second,

$$ |b\beta\Phi_{\beta}| = \frac{|\beta| |\Phi_{\beta}|}{|\beta - \alpha| |\beta - \gamma|} \le \frac{|\Phi_{\beta}|}{|\beta - \alpha|} \le \frac{|\Phi_{\beta}|}{\varepsilon} < +\infty. $$

Hence

$$ \sup_{t \in (0,1]} \sup_{y \ge 0} \sup_{\delta \in \mathcal{E} \cap [\delta_{\infty}, +\infty)} |\Phi(t, \xi, y)| < +\infty. \qquad (5.12) $$

From (5.10) and (5.12), we deduce that (5.9) actually holds true for all $t > 0$, and we are done with this case.

Step (ii). Without loss of generality, we may, to fix ideas, suppose that $\delta \to \delta_i$ ($i = 1$ or 2) provides the merging of the simple roots $\alpha$ and $\beta$ into the double root $\lambda_i$ — see Lemmas A.4 and A.5. Recall in that case that

$$ \inf_{\delta \in (\delta_i - \eta, \delta_i + \eta)} \inf_{\lambda \in [\alpha, \beta]} (|\lambda|, |\lambda - \gamma|) > \varepsilon, $$

for some $\eta > 0$. Then, for $\delta \in (\delta_i - \eta, \delta_i + \eta)$ and $\lambda \in [\alpha, \beta]$, one sets

$$ \psi_{\lambda} = \psi_{\lambda}(t, \xi, y) := \frac{\lambda(\Phi_{\lambda} - \Phi_{\gamma})}{\lambda - \gamma}, $$
It now remains to prove Lemma 5.1.

**Proof of Lemma 5.1.** The proof is related to the study of the $\delta$-indexed polynomials $P_\delta$ in Appendix A.2. As $\delta = (D-d)\|\xi\|^2$ browses $\mathbb{R}_+$, the polynomials $P_\delta$ may provide double or triple roots for, at most, two values of $\delta$ — see Lemmas A.2 to A.5 (i). In Lemmas A.2 to A.5 we build a closed set $\mathcal{E}$ of $\mathbb{R}_+$ so that, for $\delta \in \mathcal{E}$, the roots $\alpha, \beta, \gamma$ remain away from each other and so that each connected component of $\mathcal{E}^c$ contains one and only one $\delta_i$ ($i = 0, 1$ or 2) which provides a multiple root $\lambda_i$. This proof is divided into three steps:

(i) We treat the cases where $\delta \in \mathcal{E}$ for which, thanks to the boundedness of $a, b$ and $c$, we can consider $a\alpha\Phi_\alpha, b\beta\Phi_\beta$ and $c\gamma\Phi_\gamma$ independently.

We deal then with the situations where $\delta \in \mathcal{E}^c$. In those, because $a, b$ and $c$ may be unbounded we must consider $a\alpha\Phi_\alpha, b\beta\Phi_\beta$ and $c\gamma\Phi_\gamma$ together. Let $\delta \in \mathcal{E}^c$, then there is $i \in \{0, 1, 2\}$ such that $\delta \in (\delta_i - \eta, \delta_i + \eta)$ — recall that $\delta = \delta_i$ provides the multiple root $\lambda_i$ and that $\eta$ is defined in Lemmas A.3 to A.5.

(ii) We treat the cases where the root $\lambda_i$ is double.

(iii) We treat the cases where the root $\lambda_i$ is triple.

Before starting, we emphasize that (A.5) insures $\mathfrak{Re}(\alpha), \mathfrak{Re}(\beta), \mathfrak{Re}(\gamma) \le 0$ which allows us to use the estimates (A.4) at point $z = \frac{-2\cdot\sqrt{dt+y}}{2\sqrt{dt}}$ for $\cdot$ in the convex hull of $\alpha, \beta$ and $\gamma$. Notice that, in the sequel, we shall often write $\Phi_\cdot$ for $\Phi_\cdot(t, \xi, y)$.

Step (i). Because $\inf_{\delta \in \mathcal{E}} (|\alpha - \beta|, |\alpha - \gamma|, |\beta - \gamma|) > \varepsilon$, we have

$$ \sup_{\delta \in \mathcal{E}} (|a|, |b|, |c|) < \frac{1}{\varepsilon^2}. \qquad (5.8) $$

Next, using (A.4), we have, for $\lambda \in \{\alpha, \beta, \gamma\}$,

$$ |\lambda\Phi_\lambda| = \left| \lambda \frac{\text{Erfc}}{\Gamma} \left( \frac{-2\lambda\sqrt{dt} + y}{2\sqrt{dt}} \right) \right| \lesssim \frac{2|\lambda|\sqrt{dt}}{-2\lambda\sqrt{dt} + y} \le \frac{1}{\sqrt{t}}, $$

the last inequality coming from $\mathfrak{Re}(\lambda) \le 0$ and $y \ge 0$. Gathering the latter inequality with (5.8) provides, in view of (5.3),

$$ \sup_{y \ge 0} \sup_{\delta \in \mathcal{E}} |\Phi(t, \xi, y)| \lesssim \frac{1}{\sqrt{1+t}} \quad \text{for all } t \ge 1. \qquad (5.9) $$

We now turn to the case where $0 < t \le 1$ (still when $\delta \in \mathcal{E}$). Let $\delta_\infty > 0$ be defined as in Lemma A.1. It is readily seen that

$$ \sup_{t \in (0,1]} \sup_{y \ge 0} \sup_{\delta \in \mathcal{E} \cap [0, \delta_\infty]} |\Phi(t, \xi, y)| < +\infty. \qquad (5.10) $$

Next let us consider $\delta \ge \delta_\infty$. Since $\gamma = \bar{\beta}, c = \bar{b}, \Phi_\gamma = \bar{\Phi}_\beta$, we deduce from (5.3) that

$$ |\Phi| \le |a\alpha\Phi_\alpha| + 2|b\beta\Phi_\beta|. \qquad (5.11) $$
and we claim that

$$ \sup_{y \ge 0} \sup_{\delta \in (\delta_i - \eta, \delta_i + \eta)} \sup_{\lambda \in [\alpha, \beta]} |\psi'_\lambda| \lesssim \frac{1}{\sqrt{1+t}}, \quad \forall t > 0, \quad (5.13) $$

where $\psi'_\lambda$ denotes $(\partial_\bullet \psi_\bullet)_\lambda(t, \xi, y)$. Using $a\alpha + b\beta + c\gamma = 0$ we get

$$ \begin{aligned} |\Phi| &= |a\alpha(\Phi_\alpha - \Phi_\gamma) + b\beta(\Phi_\beta - \Phi_\gamma)| \\ &= \frac{|\psi_\alpha - \psi_\beta|}{|\alpha - \beta|} \\ &\le \sup_{y \ge 0} \sup_{\delta \in (\delta_i - \eta, \delta_i + \eta)} \sup_{\lambda \in [\alpha, \beta]} |\psi'_\lambda| \\ &\lesssim \frac{1}{\sqrt{1+t}}, \end{aligned} $$

where we used the mean value inequality to provide third line.

It thus only remains to prove the claim (5.13). We have

$$ \psi'_{\lambda} = \frac{\lambda \Phi'_{\lambda}}{\lambda - \gamma} - \frac{\gamma \Phi_{\lambda}}{(\lambda - \gamma)^2} + \frac{\gamma \Phi_{\gamma}}{(\lambda - \gamma)^2}, $$

where $\Phi'_\lambda$ denotes $(\partial_\bullet \Phi_\bullet)_\lambda(t, \xi, y)$. It is at first clear that

$$ \sup_{t \in (0,1]} \sup_{y \ge 0} \sup_{\delta \in (\delta_i - \eta, \delta_i + \eta)} \sup_{\lambda \in [\alpha, \beta]} |\psi'_\lambda| < +\infty. \quad (5.14) $$

Next, using control (A.4), we have, for $\lambda \in [\alpha, \beta]$,

$$ \begin{aligned} |\Phi'_{\lambda}| &= \sqrt{t} \left| \left( \frac{\mathrm{Erfc}}{\Gamma} \right)' \left( \frac{-2\lambda\sqrt{d}t + y}{2\sqrt{dt}} \right) \right| \lesssim \sqrt{t} \frac{4dt}{|-2\lambda\sqrt{d}t + y|^2} \le \frac{1}{|\lambda|^2\sqrt{t}}, \\ |\Phi_{\lambda}| &= \left| \frac{\mathrm{Erfc}}{\Gamma} \left( \frac{-2\lambda\sqrt{d}t + y}{2\sqrt{dt}} \right) \right| \lesssim \frac{2\sqrt{dt}}{-2\lambda\sqrt{d}t + y} \le \frac{1}{|\lambda|\sqrt{t}}, \\ |\gamma\Phi_{\gamma}| &= \left| \gamma \frac{\mathrm{Erfc}}{\Gamma} \left( \frac{-2\gamma\sqrt{d}t + y}{2\sqrt{dt}} \right) \right| \lesssim \frac{|\gamma|2\sqrt{dt}}{-2\gamma\sqrt{d}t + y} \le \frac{1}{\sqrt{t}}, \end{aligned} $$

and so, because $|\lambda| > \varepsilon$ — see Lemmas A.4 and A.5, (iii) —, one gets, for all $t \ge 1$,

$$ \sup_{y \ge 0} \sup_{\delta \in (\delta_i - \eta, \delta_i + \eta)} \sup_{\lambda \in [\alpha, \beta]} |\psi'_\lambda| \lesssim \frac{1}{\sqrt{1+t}}. \quad (5.15) $$

Gathering (5.14) and (5.15) finally gives (5.13) for all $t > 0$ and we are done with this case.

Step (iii). We recall in this last case that $\delta \to \delta_0$ provokes the merging of the three simple roots $\alpha, \beta, \gamma$ into the triple root $\lambda_0 = -\frac{A\sqrt{d}}{3}$. Without loss of generality it may be assumed that $\alpha$ is the real root, and taking $\varepsilon_R$ and $\varepsilon_I$ as defined in Lemma A.3, (ii), we have

$$ \alpha = \lambda_0 - 2\varepsilon_R, \quad \beta = \lambda_0 + \varepsilon_R + i\varepsilon_I, \quad \gamma = \lambda_0 + \varepsilon_R - i\varepsilon_I. $$
The first term in the above right-hand-side is clearly bounded while, for the second,

$$|b\beta\Phi_{\beta}| = \frac{|\beta| |\Phi_{\beta}|}{|\beta - \alpha| |\beta - \gamma|} \le \frac{|\Phi_{\beta}|}{|\beta - \alpha|} \le \frac{|\Phi_{\beta}|}{\varepsilon} < +\infty.$$

Hence

$$\sup_{t \in (0,1]} \sup_{y \ge 0} \sup_{\delta \in \mathcal{E} \cap [\delta_{\infty}, +\infty)} |\Phi(t, \xi, y)| < +\infty. \quad (5.12)$$

From (5.10) and (5.12), we deduce that (5.9) actually holds true for all $t > 0$, and we are done with this case.

Step (ii). Without loss of generality, we may, to fix ideas, suppose that $\delta \to \delta_i$ ($i = 1$ or 2) provides the merging of the simple roots $\alpha$ and $\beta$ into the double root $\lambda_i$ — see Lemmas A.4 and A.5. Recall in that case that

$$\inf_{\delta \in (\delta_i - \eta, \delta_i + \eta)} \inf_{\lambda \in [\alpha, \beta]} (|\lambda|, |\lambda - \gamma|) > \varepsilon,$$

for some $\eta > 0$. Then, for $\delta \in (\delta_i - \eta, \delta_i + \eta)$ and $\lambda \in [\alpha, \beta]$, one sets

$$\psi_{\lambda} = \psi_{\lambda}(t, \xi, y) := \frac{\lambda(\Phi_{\lambda} - \Phi_{\gamma})}{\lambda - \gamma},$$

and we claim that

$$\sup_{y \ge 0} \sup_{\delta \in (\delta_i - \eta, \delta_i + \eta)} \sup_{\lambda \in [\alpha, \beta]} |\psi'_{\lambda}| \lesssim \frac{1}{\sqrt{1+t}}, \quad \forall t > 0, \quad (5.13)$$

where $\psi'_{\lambda}$ denotes $(\partial_{\bullet}\psi_{\bullet})_{\lambda}(t, \xi, y)$. Using $a\alpha + b\beta + c\gamma = 0$ we get

$$\begin{aligned} |\Phi| &= |a\alpha(\Phi_{\alpha} - \Phi_{\gamma}) + b\beta(\Phi_{\beta} - \Phi_{\gamma})| \\ &= \frac{|\psi_{\alpha} - \psi_{\beta}|}{|\alpha - \beta|} \\ &\le \sup_{y \ge 0} \sup_{\delta \in (\delta_i - \eta, \delta_i + \eta)} \sup_{\lambda \in [\alpha, \beta]} |\psi'_{\lambda}| \\ &\lesssim \frac{1}{\sqrt{1+t}}, \end{aligned}$$

where we used the mean value inequality to provide third line.

It thus only remains to prove the claim (5.13). We have

$$\psi'_{\lambda} = \frac{\lambda\Phi'_{\lambda}}{\lambda - \gamma} - \frac{\gamma\Phi_{\lambda}}{(\lambda - \gamma)^2} + \frac{\gamma\Phi_{\gamma}}{(\lambda - \gamma)^2},$$

where $\Phi'_{\lambda}$ denotes $(\partial_{\bullet}\Phi_{\bullet})_{\lambda}(t, \xi, y)$. It is at first clear that

$$\sup_{t \in (0,1]} \sup_{y \ge 0} \sup_{\delta \in (\delta_i - \eta, \delta_i + \eta)} \sup_{\lambda \in [\alpha, \beta]} |\psi'_{\lambda}| < +\infty. \quad (5.14)$$

Next, using control (A.4), we have, for $\lambda \in [\alpha, \beta]$,

$$|\Phi'_{\lambda}| = \sqrt{t} \left| \left( \frac{\mathrm{Erfc}}{\Gamma} \right)' \left( \frac{-2\lambda\sqrt{d}t + y}{2\sqrt{dt}} \right) \right| \lesssim \sqrt{t} \frac{4dt}{|-2\lambda\sqrt{d}t + y|^2} \le \frac{1}{|\lambda|^2 \sqrt{t}},$$
Using again $a\alpha + b\beta + c\gamma = 0$ yields

$$
\begin{align*}
\Phi &= b\beta(\Phi_{\beta} - \Phi_{\alpha}) + c\gamma(\Phi_{\gamma} - \Phi_{\alpha}) \\
&= \frac{\beta}{\beta - \gamma} \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} - \frac{\gamma}{\beta - \gamma} \frac{\Phi_{\gamma} - \Phi_{\alpha}}{\gamma - \alpha} \\
&= \frac{\beta}{\beta - \gamma} \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} - \frac{\gamma}{\beta - \gamma} \left( \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} \right) \\
&= \mathfrak{Re} \left( \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} \right) + i \frac{\beta + \gamma}{\beta - \gamma} \mathfrak{Im} \left( \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} \right) \\
&= \mathfrak{Re} \left( \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} \right) + \frac{\lambda_0 + \varepsilon_R}{\varepsilon_I} \mathfrak{Im} \left( \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} \right). \tag{5.16}
\end{align*}
$$

The main issue is now to control second term in (5.16) because one has to compensate the term $\varepsilon_I$ which vanishes as $\delta \to \delta_0$. To do so, we use a Taylor-Lagrange expansion of $\Phi_\bullet$ at $\bullet = \alpha$ which provides, for some $\lambda \in [\alpha, \beta]$,

$$
\Phi_{\beta} = \Phi_{\alpha} + (\beta - \alpha)\Phi'_{\alpha} + \frac{(\beta - \alpha)^2}{2}\Phi''_{\lambda},
$$

where $\Phi'_\alpha$ denotes $(\partial_\bullet \Phi_\bullet)_\alpha$, and $\Phi''_\lambda$ denotes $(\partial_\bullet \Phi_\bullet)_\lambda$. Hence, one gets

$$
\frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} = \Phi'_{\alpha} + \frac{3\varepsilon_{R} + i\varepsilon_{I}}{2} \Phi''_{\lambda}.
$$

Now, because $\alpha$ is real,

$$
\mathfrak{Im} \left( \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} \right) = \frac{3\varepsilon_{R}}{2} \mathfrak{Im}(\Phi_{\lambda''}) + \frac{\varepsilon_{I}}{2} \mathfrak{Re}(\Phi_{\lambda''}),
$$

and therefore, using that $\varepsilon_R \stackrel{\delta \to \delta_0}{=} \mathcal{O}(\varepsilon_I)$,

$$
\left| \mathfrak{Im} \left( \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} \right) \right| \le \left( \frac{3}{2} |\varepsilon_{R}| + \frac{1}{2} |\varepsilon_{I}| \right) |\Phi_{\lambda}''| \lesssim |\varepsilon_{I} \Phi_{\lambda}''|.
$$

As a consequence, by using the mean value inequality to control the first term in (5.16), one gets

$$
|\Phi| \lesssim \sup_{y \ge 0} \sup_{\delta \in (\delta_0 - \eta, \delta_0 + \eta)} \sup_{\ell \in [\alpha, \beta]} |\Phi'_\ell| + \sup_{y \ge 0} \sup_{\delta \in (\delta_0 - \eta, \delta_0 + \eta)} \sup_{\lambda \in [\alpha, \beta]} |\lambda_0 + \varepsilon_R| |\Phi_\lambda''|.
$$

It is first clear that

$$
\sup_{t \in (0,1]} \sup_{y \ge 0} \sup_{\delta \in (\delta_0 - \eta, \delta_0 + \eta)} \sup_{\ell \in [\alpha, \beta]} |\Phi'_{\ell}| < +\infty, \quad (5.17)
$$

and

$$
\sup_{t \in (0,1]} \sup_{y \ge 0} \sup_{\delta \in (\delta_0 - \eta, \delta_0 + \eta)} \sup_{\lambda \in [\alpha, \beta]} |\Phi_{\lambda}''| < +\infty. \quad (5.18)
$$

Now, using control (A.4), we have, for $\ell$ and $\lambda \in [\alpha, \beta]$,

$$
|\Phi'_{\ell}| = \sqrt{t} \left| \left( \frac{\mathrm{Erfc}}{\Gamma} \right)' \left( \frac{-2\ell\sqrt{d}t + y}{2\sqrt{dt}} \right) \right| \lesssim \sqrt{t} \frac{4dt}{|-2\ell\sqrt{d}t + y|^2} \le \frac{1}{|\ell|^2\sqrt{t}},
$$
$$
\begin{align*}
|\Phi_{\lambda}| &= \left| \frac{\mathrm{Erfc}}{\Gamma} \left( \frac{-2\lambda\sqrt{d}t + y}{2\sqrt{dt}} \right) \right| \lesssim \frac{2\sqrt{dt}}{|-2\lambda\sqrt{d}t + y|} \le \frac{1}{|\lambda|\sqrt{t}}, \\
|\gamma\Phi_{\gamma}| &= \left| \gamma \frac{\mathrm{Erfc}}{\Gamma} \left( \frac{-2\gamma\sqrt{d}t + y}{2\sqrt{dt}} \right) \right| \lesssim \frac{|\gamma|2\sqrt{dt}}{|-2\gamma\sqrt{d}t + y|} \le \frac{1}{\sqrt{t}},
\end{align*}
$$

and so, because $|\lambda| > \varepsilon$ — see Lemmas A.4 and A.5, (iii) —, one gets, for all $t \ge 1$,

$$
\sup_{y \ge 0} \sup_{\delta \in (\delta_i - \eta, \delta_i + \eta)} \sup_{\lambda \in [\alpha, \beta]} |\psi'_{\lambda}| \lesssim \frac{1}{\sqrt{1+t}}. \qquad (5.15)
$$

Gathering (5.14) and (5.15) finally gives (5.13) for all $t > 0$ and we are done with this case.

Step (iii). We recall in this last case that $\delta \to \delta_0$ provokes the merging of the three simple roots $\alpha, \beta, \gamma$ into the triple root $\lambda_0 = -\frac{A\sqrt{d}}{3}$. Without loss of generality it may be assumed that $\alpha$ is the real root, and taking $\varepsilon_R$ and $\varepsilon_I$ as defined in Lemma A.3, (ii), we have

$$
\alpha = \lambda_0 - 2\varepsilon_R, \quad \beta = \lambda_0 + \varepsilon_R + i\varepsilon_I, \quad \gamma = \lambda_0 + \varepsilon_R - i\varepsilon_I.
$$

Using again $a\alpha + b\beta + c\gamma = 0$ yields

$$
\begin{align*}
\Phi &= b\beta(\Phi_{\beta} - \Phi_{\alpha}) + c\gamma(\Phi_{\gamma} - \Phi_{\alpha}) \\
&= \frac{\beta}{\beta - \gamma} \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} - \frac{\gamma}{\beta - \gamma} \frac{\Phi_{\gamma} - \Phi_{\alpha}}{\gamma - \alpha} \\
&= \frac{\beta}{\beta - \gamma} \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} - \frac{\gamma}{\beta - \gamma} \overline{\left( \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} \right)} \\
&= \mathfrak{Re} \left( \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} \right) + i \frac{\beta + \gamma}{\beta - \gamma} \mathfrak{Im} \left( \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} \right) \\
&= \mathfrak{Re} \left( \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} \right) + \frac{\lambda_0 + \varepsilon_R}{\varepsilon_I} \mathfrak{Im} \left( \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} \right). \tag{5.16}
\end{align*}
$$

The main issue is now to control second term in (5.16) because one has to compensate the
term $\varepsilon_I$ which vanishes as $\delta \to \delta_0$. To do so, we use a Taylor-Lagrange expansion of $\Phi_\bullet$ at
$\bullet = \alpha$ which provides, for some $\lambda \in [\alpha, \beta]$,

$$
\Phi_{\beta} = \Phi_{\alpha} + (\beta - \alpha)\Phi'_{\alpha} + \frac{(\beta - \alpha)^2}{2} \Phi''_{\lambda},
$$

where $\Phi'_\alpha$ denotes $(\partial_\bullet \Phi_\bullet)_\alpha$, and $\Phi''_\lambda$ denotes $(\partial_{\bullet\bullet} \Phi_\bullet)_\lambda$. Hence, one gets

$$
\frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} = \Phi'_{\alpha} + \frac{3\varepsilon_{R} + i\varepsilon_{I}}{2} \Phi''_{\lambda}.
$$

Now, because $\alpha$ is real,

$$
\mathfrak{Im} \left( \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} \right) = \frac{3\varepsilon_{R}}{2} \mathfrak{Im}(\Phi_{\lambda}^{\prime\prime}) + \frac{\varepsilon_{I}}{2} \mathfrak{Re}(\Phi_{\lambda}^{\prime\prime}),
$$

and therefore, using that $\varepsilon_R \stackrel{\delta \to \delta_0}{=} \mathcal{O}(\varepsilon_I)$,

$$
\left| \mathfrak{Im} \left( \frac{\Phi_{\beta} - \Phi_{\alpha}}{\beta - \alpha} \right) \right| \le \left( \frac{3}{2} |\varepsilon_{R}| + \frac{1}{2} |\varepsilon_{I}| \right) |\Phi_{\lambda}''| \lesssim |\varepsilon_{I} \Phi_{\lambda}''|.
$$
$$|\Phi_{\lambda}''| = t \left| \left( \frac{\mathrm{Erfc}}{\Gamma} \right)'' \left( \frac{-2\lambda\sqrt{d}t + y}{2\sqrt{dt}} \right) \right| \lesssim t \frac{8(dt)^{\frac{3}{2}}}{|-2\lambda\sqrt{d}t + y|^3} \le \frac{1}{|\lambda|^3\sqrt{t}},$$

and thus, because $|\ell|$ and $|\lambda| > \varepsilon$ — see Lemma A.3, (iv) —, one gets, for all $t \ge 1$,

$$\sup_{y \ge 0} \sup_{\delta \in (\delta_0 - \eta, \delta_0 + \eta)} \sup_{\ell \in [\alpha, \beta]} |\Phi'_{\ell}| \lesssim \frac{1}{\sqrt{1+t}}, \qquad (5.19)$$

$$\sup_{y \ge 0} \sup_{\delta \in (\delta_0 - \eta, \delta_0 + \eta)} \sup_{\lambda \in [\alpha, \beta]} |\Phi_{\lambda}''| \lesssim \frac{1}{\sqrt{1+t}}. \qquad (5.20)$$

Gathering (5.17), (5.18), (5.19) and (5.20) finally provides $|\Phi| \lesssim \frac{1}{\sqrt{1+t}}$ for all $t > 0$, and we are done with this case which concludes the proof of the lemma. $\square$
As a consequence, by using the mean value inequality to control the first term in (5.16), one gets

$$|\Phi| \lesssim \sup_{y \ge 0} \sup_{\delta \in (\delta_0 - \eta, \delta_0 + \eta)} \sup_{\ell \in [\alpha, \beta]} |\Phi'_\ell| + \sup_{y \ge 0} \sup_{\delta \in (\delta_0 - \eta, \delta_0 + \eta)} \sup_{\lambda \in [\alpha, \beta]} |\lambda_0 + \varepsilon_R| |\Phi''_\lambda|.$$

It is first clear that

$$\sup_{t \in (0,1]} \sup_{y \ge 0} \sup_{\delta \in (\delta_0 - \eta, \delta_0 + \eta)} \sup_{\ell \in [\alpha, \beta]} |\Phi'_\ell| < +\infty, \quad (5.17)$$

and

$$\sup_{t \in (0,1]} \sup_{y \ge 0} \sup_{\delta \in (\delta_0 - \eta, \delta_0 + \eta)} \sup_{\lambda \in [\alpha, \beta]} |\Phi''_\lambda| < +\infty. \quad (5.18)$$

Now, using control (A.4), we have, for $\ell$ and $\lambda \in [\alpha, \beta]$,

$$|\Phi'_{\ell}| = \sqrt{t} \left| \left( \frac{\mathrm{Erfc}}{\Gamma} \right)' \left( \frac{-2\ell\sqrt{d}t + y}{2\sqrt{dt}} \right) \right| \lesssim \sqrt{t} \frac{4dt}{|-2\ell\sqrt{d}t + y|^2} \le \frac{1}{|\ell|^2\sqrt{t}},$$

$$|\Phi''_{\lambda}| = t \left| \left( \frac{\mathrm{Erfc}}{\Gamma} \right)'' \left( \frac{-2\lambda\sqrt{d}t + y}{2\sqrt{dt}} \right) \right| \lesssim t \frac{8(dt)^{\frac{3}{2}}}{|-2\lambda\sqrt{d}t + y|^3} \le \frac{1}{|\lambda|^3\sqrt{t}},$$

and thus, because $|\ell|$ and $|\lambda| > \varepsilon$ — see Lemma A.3, (iv) —, one gets, for all $t \ge 1$,

$$\sup_{y \ge 0} \sup_{\delta \in (\delta_0 - \eta, \delta_0 + \eta)} \sup_{\ell \in [\alpha, \beta]} |\Phi'_{\ell}| \lesssim \frac{1}{\sqrt{1+t}}, \quad (5.19)$$

$$\sup_{y \ge 0} \sup_{\delta \in (\delta_0 - \eta, \delta_0 + \eta)} \sup_{\lambda \in [\alpha, \beta]} |\Phi''_{\lambda}| \lesssim \frac{1}{\sqrt{1+t}}. \quad (5.20)$$

Gathering (5.17), (5.18), (5.19) and (5.20) finally provides $|\Phi| \lesssim \frac{1}{\sqrt{1+t}}$ for all $t > 0$, and we are done with this case which concludes the proof of the lemma. $\square$
# 6 Numerical explorations

Through this section, we assume $N = 2$, meaning that the road is a line and we perform some numerical simulations of the Cauchy problem (1.1)–(2.11). To do so, we truncate the unbounded domain $\mathbb{R} \times (0, +\infty)$ and work on a bounded box $(-2M, 2M) \times (0, M)$, where $M > 0$. To preserve the quantity of individuals, we impose the no-flux boundary conditions through the artificial frontiers. The resolution method is based on a classical finite difference scheme. Since we consider $M \gg 1$ and an initial datum $(v_0, u_0)$ supported in a “relatively small” compact compared to the size of the box, we are confident that the numerical solution should be close to the real solution.

**On the decay rate.** When $D > d$, the asymptotic decay rate of the $L^\infty$ norm of the solution $(v, u)$ is expected to be of the magnitude $O((1+t)^{-1})$, see Theorem 2.2 and the discussion on the logarithmic term right after. This is numerically confirmed by the left panel of Figure I.

On the other hand, when $D \le d$ (which is not the essence of the field-road model), one should, based on Remark 4.1, redo the lengthy computations of Section 5 to capture the asymptotic decay rate. We are rather confident that this should not alter the result, and this is sustained by numerical simulations, see the right panel of Figure I.

<figure><img src="image_2.png" /><figcaption>**Figure I** — Decay in log-scale of the solution to the Cauchy problem (1.1)–(2.11) with $N = 2$, $M = 800$, $\mu = 10$ and $\nu = 1$, arising from the datum $(v_0, u_0) = (\mathbb{1}_{[-10,10]} \times [10,30], \mathbb{1}_{[-10,10]})$.</figcaption></figure>

**On the level sets.** We now turn to numerically explore the form of the level sets of the solution $(v, u)$. The results are presented in Figure II.

In particular on the bottom panel of Figure II, corresponding to a situation where $D = 100$ and $d = 1$, we observe the following. There holds

$$-\partial_y v|_{y=0}(t=1000, \cdot) < 0 \quad \text{in the "middle of the road"},$$
# 6 Numerical explorations

Through this section, we assume $N = 2$, meaning that the road is a line and we perform some numerical simulations of the Cauchy problem (1.1)–(2.11). To do so, we truncate the unbounded domain $\mathbb{R} \times (0, +\infty)$ and work on a bounded box $(-2M, 2M) \times (0, M)$, where $M > 0$. To preserve the quantity of individuals, we impose the no-flux boundary conditions through the artificial frontiers. The resolution method is based on a classical finite difference scheme. Since we consider $M \gg 1$ and an initial datum $(v_0, u_0)$ supported in a “relatively small” compact compared to the size of the box, we are confident that the numerical solution should be close to the real solution.

**On the decay rate.** When $D > d$, the asymptotic decay rate of the $L^\infty$ norm of the solution $(v, u)$ is expected to be of the magnitude $O((1+t)^{-1})$, see Theorem 2.2 and the discussion on the logarithmic term right after. This is numerically confirmed by the left panel of Figure I.

On the other hand, when $D \le d$ (which is not the essence of the field-road model), one should, based on Remark 4.1, redo the lengthy computations of Section 5 to capture the asymptotic decay rate. We are rather confident that this should not alter the result, and this is sustained by numerical simulations, see the right panel of Figure I.

<figure><img src="image_4.png" /><figcaption>**Figure I** — Decay in log-scale of the solution to the Cauchy problem (1.1)–(2.11) with $N = 2$, $M = 800$, $\mu = 10$ and $\nu = 1$, arising from the datum $(v_0, u_0) = (\mathbb{1}_{[-10,10]} \times [10,30], \mathbb{1}_{[-10,10]}).$</figcaption></figure>

**On the level sets.** We now turn to numerically explore the form of the level sets of the solution $(v, u)$. The results are presented in Figure II.

In particular on the bottom panel of Figure II, corresponding to a situation where $D = 100$ and $d = 1$, we observe the following. There holds

$$-\partial_y v|_{y=0}(t = 1000, \cdot) < 0 \quad \text{in the "middle of the road"},$$
meaning that the individuals mainly switch from the field to the road. On the other hand, there holds

$$-\partial_y v|_{y=0}(t=1000, \cdot) > 0 \text{ "far away" on the road,}$$

meaning that the individuals mainly switch from the road to the field. Roughly speaking, the road sucks up individuals in the central region (corresponding to the bulk of the population) and spits them out in the far away region (corresponding to the tails of the population).

On the other hand, on the top panel of Figure II, corresponding to a situation where $D = 0.1$ and $d = 1$, there holds $-\partial_y v|_{y=0}(t=1000, \cdot) \approx 0$, meaning that there are very little exchanges between the field and the road.

<figure><img src="image_3.png" /><figcaption>**Figure II** — Two snapshots at time $t = 1000$ of a few level sets of the solution to the Cauchy problem (1.1)–(2.11) with $N = 2$, $M = 200$, $d = \mu = \nu = 1$ and two different values for $D$, arising from the datum $(v_0, u_0) = (1_{[-5,5]} \times [0,5], 0)$.</figcaption></figure>

**On the flux.** In view of the above considerations, we are now interested in

$$F(t, x) := \mu u(t, x) - \nu v|_{y=0}(t, x) = -d \partial_y v|_{y=0}(t, x)$$

which is the flux entering the field.

When $D = 100$ and $d = 1$, Figure III confirms that the flux is negative in the middle of the road, and positive far away. Furthermore if we denote $x_0(t)$ the (rightmost) position where the flux changes sign, that is $F(t, x_0(t)) = 0$, the left panel of Figure IV suggests that $x_0(t)$ asymptotically behaves like $O(t^{1/2})$. Last, the right panel of Figure IV indicates that the asymptotic decay rate of the flux at $x = 0$ is of the magnitude $O((1 + t)^{-3/2})$.
meaning that the individuals mainly switch from the field to the road. On the other hand, there holds

$$ -\partial_y v|_{y=0}(t = 1000, \cdot) > 0 \text{ "far away" on the road,} $$

meaning that the individuals mainly switch from the road to the field. Roughly speaking, the road sucks up individuals in the central region (corresponding to the bulk of the population) and spits them out in the far away region (corresponding to the tails of the population).

On the other hand, on the top panel of Figure II, corresponding to a situation where $D = 0.1$ and $d = 1$, there holds $-\partial_y v|_{y=0}(t = 1000, \cdot) \approx 0$ (since the level sets of $v$ are almost perpendicular to the road), meaning that there are very few exchanges between the field and the road. This can be observed with further accuracy in Figure V.

<figure><img src="image_5.png" /><figcaption>**Figure II** — Two snapshots at time $t = 1000$ of a few level sets of the solution to the Cauchy problem (1.1)–(2.11) with $N = 2$, $M = 200$, $d = \mu = \nu = 1$ and two different values for $D$, arising from the datum $(v_0, u_0) = (\mathbb{1}_{[-5,5]} \times [0,5], 0)$.</figcaption></figure>

**On the flux.** In view of the above considerations, we are now interested in

$$ F(t, x) := \mu u(t, x) - \nu v|_{y=0}(t, x) = -d \partial_y v|_{y=0}(t, x) $$

which is the flux entering the field.

When $D = 100$ and $d = 1$, Figure III confirms that the flux is negative in the middle of the road, and positive far away. Furthermore if we denote $x_0(t)$ the (rightmost) position where the flux changes sign, that is $F(t, x_0(t)) = 0$, the left panel of Figure IV suggests that $x_0(t)$ asymptotically behaves like $O(t^{1/2})$. Last, the right panel of Figure IV indicates that the asymptotic decay rate of the flux at $x = 0$ is of the magnitude $O((1 + t)^{-3/2})$.
<figure><img src="image_4.png" /><figcaption>**Figure III** — Plot of the flux entering the field, $F(t, x) = \mu u(t, x) - \nu v|_{y=0}(t, x)$, of the solution to the Cauchy problem (1.1)–(2.11) with $N=2$, $M=200$, $d=\mu=\nu=1$ and $D=100$, arising from the datum $(v_0, u_0) \equiv (\mathbb{1}_{[-5,5]} \times [0,5], 0)$, for eight different values of time.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure IV** — Plots in log-scale of $x_0(t)$ (on left) and $|F(t, 0)|$ (on right) for the solution to the Cauchy problem (1.1)–(2.11) with $N=2$, $M=200$, $d=\mu=\nu=1$ and $D=100$, arising from the datum $(v_0, u_0) \equiv (\mathbb{1}_{[-5,5]} \times [0,5], 0)$.</figcaption></figure>

When $D = 0.1$, $d = 1$, Figure V indicates that the flux is positive in the middle of the road, and negative far away, that is the opposite situation than the previous case. Furthermore if we denote $x_0(t)$ the (rightmost) position where the flux changes sign, that is $F(t, x_0(t)) = 0$, the left panel of Figure VI suggest that $x_0(t)$, again, asymptotically behaves like $O(t^{1/2})$. Last, the right panel of Figure VI indicates that the asymptotic decay rate of the flux at $x=0$ is of the magnitude $O((1+t)^{-2})$, which is in contrast with the previous case. In particular, the decay of the flux is faster, which is consistent with the observation $-\partial_y v|_{y=0}(t=1000, \cdot) \approx 0$ on the top panel of Figure II.
<figure><img src="image_6.png" /><figcaption>**Figure III** — Plot of the flux entering the field, $F(t, x) = \mu u(t, x) - \nu v|_{y=0}(t, x)$, of the solution to the Cauchy problem (1.1)–(2.11) with $N = 2$, $M = 200$, $d = \mu = \nu = 1$ and $D = 100$, arising from the datum $(v_0, u_0) \equiv (\mathbb{1}_{[-5,5]} \times [0,5], 0)$, for eight different values of time.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure IV** — Plots in log-scale of $x_0(t)$ (on left) and $|F(t, 0)|$ (on right) for the solution to the Cauchy problem (1.1)–(2.11) with $N = 2$, $M = 200$, $d = \mu = \nu = 1$ and $D = 100$, arising from the datum $(v_0, u_0) \equiv (\mathbb{1}_{[-5,5]} \times [0,5], 0)$.</figcaption></figure>

When $D = 0.1$, $d = 1$, Figure V indicates that the flux is positive in the middle of the road, and negative far away, that is the opposite situation than the previous case. Furthermore if we denote $x_0(t)$ the (rightmost) position where the flux changes sign, that is $F(t, x_0(t)) = 0$, the left panel of Figure VI suggest that $x_0(t)$, again, asymptotically behaves like $O(t^{1/2})$. Last, the right panel of Figure VI indicates that the asymptotic decay rate of the flux at $x=0$ is of the magnitude $O((1+t)^{-2})$, which is in contrast with the previous case.
<figure><img src="image_6.png" /><figcaption>**Figure V** — Plot of the flux entering the field, $F(t, x) = \mu u(t, x) - \nu v|_{y=0}(t, x)$, of the solution to the Cauchy problem (1.1)–(2.11) with $N = 2$, $M = 200$, $d = \mu = \nu = 1$ and $D = 0.1$, arising from the datum $(v_0, u_0) \equiv (\mathbb{1}_{[-5,5] \times [0,5]}, 0)$, for eight different values of time.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure VI** — Plots in log-scale of $x_0(t)$ (on left) and $|F(t, 0)|$ (on right) for the solution to the Cauchy problem (1.1)–(2.11) with $N = 2$, $M = 200$, $d = \mu = \nu = 1$ and $D = 0.1$, arising from the datum $(v_0, u_0) \equiv (\mathbb{1}_{[-5,5] \times [0,5]}, 0)$.</figcaption></figure>

**On the threshold value for the shape of the flux.** From the above, the shape of the flux function $F$ is very dependent on the value of $D$ ($d = 1$ being fixed). As seen in Figure VII, the threshold value should be $D = 2$. We suspect that the reason for this is that the road is one dimensional while the field is two dimensional. Indeed, when $D > 2$, the diffusion “normalized by dimension” is larger on the road, which induces a lack of individuals in the middle of the road and, thus, a negative flux entering the field — see the left panel of Figure VII.
<figure><img src="image_8.png" /><figcaption>**Figure V** — Plot of the flux entering the field, $F(t, x) = \mu u(t, x) - \nu v|_{y=0}(t, x)$, of the solution to the Cauchy problem (1.1)–(2.11) with $N = 2$, $M = 200$, $d = \mu = \nu = 1$ and $D = 0.1$, arising from the datum $(v_0, u_0) \equiv (\mathbb{1}_{[-5,5]} \times [0,5], 0)$, for eight different values of time.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure VI** — Plots in log-scale of $x_0(t)$ (on left) and $|F(t, 0)|$ (on right) for the solution to the Cauchy problem (1.1)–(2.11) with $N = 2$, $M = 200$, $d = \mu = \nu = 1$ and $D = 0.1$, arising from the datum $(v_0, u_0) \equiv (\mathbb{1}_{[-5,5]} \times [0,5], 0)$.</figcaption></figure>

**On the threshold value for the shape of the flux.** From the above, the shape of the flux function $F$ is very dependent on the value of $D$ ($d = 1$ being fixed). As we may see in Figure VII, it seems there is a threshold value at $D = 2$. This could be related to the enhanced spreading speed for the Fisher-KPP system (1.2), noticed in the original paper [10], when $D > 2d$. It would be interesting to have a result that highlights some connections between the purely diffusive system (1.1) from one side, and the nonlinear system (1.2) from the other side, especially to identify more clearly the reason for this threshold.
<figure><img src="image_8.png" /><figcaption>**Figure VII** — Plot of $F(3000, 0)$ (on left) and the decay rate of $|F(t, 0)|$ (on right) with different values of $D$ for the solution to the Cauchy problem (1.1)–(2.11) with $N = 2$, $M = 300$, $d = \mu = \nu = 1$, arising from the datum $(v_0, u_0) \equiv (1_{[-20,20]} \times [0,20], 0)$.</figcaption></figure>

Let us make a concluding remark on the asymptotic decay rate of the flux. Notice that we have focused on the flux at $x = 0$. However, the $L^\infty$ norm of the flux is not always attained at this point. In particular, we have numerically observed that, in the vicinity of $x = 0$, the flux is “dromedary shaped” (one bump) when $D$ is much smaller than 2, but “camel shaped” (two bumps) when $D$ is slightly smaller than 2. This oscillation phenomenon would deserve further investigations.

# A Appendix

## A.1 The complementary error function $Erfc$

Recalling that $\Gamma(\ell) = e^{-\ell^2}$, the Erfc function is defined by

$$ \mathrm{Erfc}(x) := \frac{2}{\sqrt{\pi}} \int_{x}^{+\infty} \Gamma(y) \, dy, \quad x \in \mathbb{R}, \qquad (\mathrm{A.1}) $$

and has an holomorphic continuation, namely

$$ \mathrm{Erfc}(z) := 1 - \frac{2}{\sqrt{\pi}} \sum_{k=0}^{+\infty} \frac{(-1)^k}{(2k+1)k!} z^{2k+1}, \quad z \in \mathbb{C}. \qquad (\mathrm{A.2}) $$

For any small $\delta > 0$, we have the asymptotic expansion, see [19, page 393] or [14, page 262],

$$ \sqrt{\pi} \frac{\mathrm{Erfc}}{\Gamma}(z) = \frac{1}{z} - \frac{1}{2z^3} + \frac{3}{4z^5} + o\left(\frac{1}{|z|^6}\right), \quad \text{as } |z| \to +\infty, \qquad (\mathrm{A.3}) $$

uniformly in $\{z \in \mathbb{C} : |\arg(z)| \le \frac{3\pi}{4} - \delta\}$ — see Figure IX. By computing the derivatives of $\frac{\mathrm{Erfc}}{\Gamma}$ and using (A.3) one can check that, for all non-zero $z \in \mathbb{R}_+ + i\mathbb{R}$,

$$ \left| \frac{\mathrm{Erfc}}{\Gamma}(z) \right| \lesssim \frac{1}{|z|}, \qquad \left| \frac{d}{dz} \frac{\mathrm{Erfc}}{\Gamma}(z) \right| \lesssim \frac{1}{|z|^2}, \qquad \left| \frac{d^2}{dz^2} \frac{\mathrm{Erfc}}{\Gamma}(z) \right| \lesssim \frac{1}{|z|^3}. \qquad (\mathrm{A.4}) $$
<figure><img src="image_10.png" /><figcaption>**Figure VII** — Plot of $F(3000, 0)$ (on left) and the decay rate of $|F(t, 0)|$ (on right) with different values of $D$ for the solution to the Cauchy problem (1.1)–(2.11) with $N = 2$, $M = 300$, $d = \mu = \nu = 1$, arising from the datum $(v_0, u_0) \equiv (1_{[-20,20]} \times [0,20], 0)$.</figcaption></figure>

Let us make a concluding remark on the asymptotic decay rate of the flux. Notice that we have focused on the flux at $x = 0$. However, the $L^\infty$ norm of the flux is not always attained at this point. In particular, we have numerically observed that, in the vicinity of $x = 0$, the flux is “dromedary shaped” (one bump) when $D$ is much smaller than 2, but “camel shaped” (two bumps) when $D$ is slightly smaller than 2. This oscillation phenomenon would deserve further investigations.

# A Appendix

## A.1 The complementary error function $Erfc$

Recalling that $\Gamma(\ell) = e^{-\ell^2}$, the Erfc function is defined by

$$ \mathrm{Erfc}(x) := \frac{2}{\sqrt{\pi}} \int_{x}^{+\infty} \Gamma(y) \, dy, \quad x \in \mathbb{R}, \qquad (\mathrm{A.1}) $$

and has an holomorphic continuation, namely

$$ \mathrm{Erfc}(z) := 1 - \frac{2}{\sqrt{\pi}} \sum_{k=0}^{+\infty} \frac{(-1)^k}{(2k+1)k!} z^{2k+1}, \quad z \in \mathbb{C}. \qquad (\mathrm{A.2}) $$

For any small $\delta > 0$, we have the asymptotic expansion, see [19, page 393] or [14, page 262],

$$ \sqrt{\pi} \frac{\mathrm{Erfc}}{\Gamma}(z) = \frac{1}{z} - \frac{1}{2z^3} + \frac{3}{4z^5} + o\left(\frac{1}{|z|^6}\right), \quad \text{as } |z| \to +\infty, \qquad (\mathrm{A.3}) $$

uniformly in $\{z \in \mathbb{C} : |\arg(z)| \le \frac{3\pi}{4} - \delta\}$ — see Figure IX. By computing the derivatives of $\frac{\mathrm{Erfc}}{\Gamma}$ and using (A.3) one can check that, for all non-zero $z \in \mathbb{R}_+ + i\mathbb{R}$,

$$ \left| \frac{\mathrm{Erfc}}{\Gamma}(z) \right| \lesssim \frac{1}{|z|}, \qquad \left| \frac{d}{dz} \frac{\mathrm{Erfc}}{\Gamma}(z) \right| \lesssim \frac{1}{|z|^2}, \qquad \left| \frac{d^2}{dz^2} \frac{\mathrm{Erfc}}{\Gamma}(z) \right| \lesssim \frac{1}{|z|^3}. \qquad (\mathrm{A.4}) $$
<figure><img src="image_9.png" /><figcaption>Figure VIII — Plot of $|\mathrm{Erfc}(z)|$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure IX — Plot of $|\frac{\mathrm{Erfc}}{\Gamma}(z)|$.</figcaption></figure>

## A.2 The polynomials $P_\delta$

Let $\mu, \nu, d > 0$ be given. For $\delta \ge 0$, we consider the polynomial

$$P_{\delta}(\sigma) := \sigma^3 + A\sqrt{d}\sigma^2 + (\mu + \delta)\sigma + A\sqrt{d}\delta,$$

where $A := \frac{\nu}{d}$. We denote $\alpha(\delta)$, $\beta(\delta)$ and $\gamma(\delta)$ its three complex roots. Notice that $0$ is a root if and only if $\delta = 0$. Recall that $P_\delta$ has multiple roots if and only if the discriminant

$$\Delta_{\delta} := 18A^{2}d(\mu + \delta)\delta - 4A^{4}d^{2}\delta + A^{2}d(\mu + \delta)^{2} - 4(\mu + \delta)^{3} - 27A^{2}d\delta^{2}$$

vanishes. As a consequence, $\alpha, \beta$ and $\gamma$ are distinct for almost all $\delta \ge 0$.

As easily checked, the real roots belong to $(-A\sqrt{d}, 0]$ and, since $\alpha + \beta + \gamma = -A\sqrt{d}$, we have

$$-A\sqrt{d} < \mathfrak{Re}(\alpha(\delta)), \mathfrak{Re}(\beta(\delta)), \mathfrak{Re}(\gamma(\delta)) \le 0, \quad \text{for all } \delta \ge 0. \quad (A.5)$$

Then, a more detailed analysis provides further information about the nature and the behaviour of the roots, as stated in the following lemmas.

**Lemma A.1 (Large values of $\delta$)** There is $\delta_\infty > 0$ large enough such that, for all $\delta \ge \delta_\infty$, the three roots are simple, one is real (say $\alpha$) and the two others are complex conjugated (say $\mathfrak{Im}(\beta) > 0$ and $\mathfrak{Im}(\gamma) < 0$). Moreover, we have

$$\lim_{\delta \to +\infty} \alpha(\delta) = -A\sqrt{d}, \quad \lim_{\delta \to +\infty} \mathfrak{Re}(\beta(\delta)) = 0, \quad \lim_{\delta \to +\infty} \mathfrak{Im}(\beta(\delta)) = +\infty,$$

so that, up to increasing $\delta_\infty$ if necessary, $|\beta - \gamma| > |\beta|$ for all $\delta \ge \delta_\infty$.

Next we fix $\nu > 0$ and $d > 0$ (and thus $A$) and describe the situations while decreasing the parameter $\mu$ from $+\infty$ to $0$. Most of the results can be seen at a glance on Figure X.
<figure><img src="image_11.png" /><figcaption>Figure VIII — Plot of $|\mathrm{Erfc}(z)|$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure IX — Plot of $|\frac{\mathrm{Erfc}}{\Gamma}(z)|$.</figcaption></figure>

## A.2 The polynomials $P_\delta$

Let $\mu, \nu, d > 0$ be given. For $\delta \ge 0$, we consider the polynomial

$$P_{\delta}(\sigma) := \sigma^3 + A\sqrt{d}\sigma^2 + (\mu + \delta)\sigma + A\sqrt{d}\delta,$$

where $A := \frac{\nu}{d}$. We denote $\alpha(\delta)$, $\beta(\delta)$ and $\gamma(\delta)$ its three complex roots. Notice that $0$ is a root if and only if $\delta = 0$. Recall that $P_\delta$ has multiple roots if and only if the discriminant

$$\Delta_{\delta} := 18A^{2}d(\mu + \delta)\delta - 4A^{4}d^{2}\delta + A^{2}d(\mu + \delta)^{2} - 4(\mu + \delta)^{3} - 27A^{2}d\delta^{2}$$

vanishes. As a consequence, $\alpha, \beta$ and $\gamma$ are distinct for almost all $\delta \ge 0$.

As easily checked, the real roots belong to $(-A\sqrt{d}, 0]$ and, since $\alpha + \beta + \gamma = -A\sqrt{d}$, we have

$$-A\sqrt{d} < \mathfrak{Re}(\alpha(\delta)), \mathfrak{Re}(\beta(\delta)), \mathfrak{Re}(\gamma(\delta)) \le 0, \quad \text{for all } \delta \ge 0. \quad (A.5)$$

Then, a more detailed analysis provides further information about the nature and the behavior of the roots, as stated in the following lemmas.

**Lemma A.1 (Large values of $\delta$)** There is $\delta_\infty > 0$ large enough such that, for all $\delta \ge \delta_\infty$, the three roots are simple, one is real (say $\alpha$) and the two others are complex conjugated (say $\mathfrak{Im}(\beta) > 0$ and $\mathfrak{Im}(\gamma) < 0$). Moreover, we have

$$\lim_{\delta \to +\infty} \alpha(\delta) = -A\sqrt{d}, \quad \lim_{\delta \to +\infty} \mathfrak{Re}(\beta(\delta)) = 0, \quad \lim_{\delta \to +\infty} \mathfrak{Im}(\beta(\delta)) = +\infty,$$

so that, up to increasing $\delta_\infty$ if necessary, $|\beta - \gamma| > |\beta|$ for all $\delta \ge \delta_\infty$.

Next we fix $\nu > 0$ and $d > 0$ (and thus $A$) and describe the situations while decreasing the parameter $\mu$ from $+\infty$ to $0$. Most of the results can be seen at a glance on Figure X.
<figure><img src="image_11.png" /><figcaption>**Figure X** — Multiplicity of the roots of $P_\delta$ for fixed $\nu$, $d$. Fix the last parameter $\mu > 0$, then $\delta$ browses $[0, +\infty)$.</figcaption></figure>

**Lemma A.2** (Only simple roots) Assume $\mu > \frac{8A^2d}{27}$ is fixed and set $\mathcal{E} := \mathbb{R}_+$.

(i) For all $\delta \in \mathcal{E}$, the three roots are simple, one is real and the two others are complex conjugated.

(ii) There is $\varepsilon > 0$ such that $\inf_{\delta \in \mathcal{E}} (|\alpha - \beta|, |\alpha - \gamma|, |\beta - \gamma|) > \varepsilon$.

**Lemma A.3** (Appearance of a triple root) Assume $\mu = \frac{8A^2d}{27}$ and set $\delta_0 := \frac{A^2d}{27}$.

(i) If $\delta \neq \delta_0$, the three roots are simple, one is real and the two others are complex conjugated, and as $\delta$ encounters $\delta_0$, the three simple roots merge in a triple root $\lambda_0 = -\frac{A\sqrt{d}}{3}$.

Without loss of generality, we assume that $\alpha$ is the real root for all $\delta \ge 0$.

(ii) Letting $\varepsilon_R := \Re(\beta) - \lambda_0$ and $\varepsilon_I := \Im(\beta)$, one has

$$ \alpha = \lambda_0 - 2\varepsilon_R, \quad \beta = \lambda_0 + \varepsilon_R + i\varepsilon_I, \quad \gamma = \lambda_0 + \varepsilon_R - i\varepsilon_I, $$

and $\varepsilon_R \stackrel{\delta \to \delta_0}{=} \mathcal{O}(\varepsilon_I)$.

Define $\eta := \frac{\delta_0}{2}$ and $\mathcal{E} := \mathbb{R}_+ \setminus (\delta_0 - \eta, \delta_0 + \eta)$, then there is $\varepsilon > 0$ such that

(iii) $\inf_{\delta \in \mathcal{E}} (|\alpha - \beta|, |\alpha - \gamma|, |\beta - \gamma|) > \varepsilon$,

(iv) $\inf_{\delta \in \mathcal{E}^c} \inf_{\lambda \in [\alpha, \beta]} |\lambda| > \varepsilon$.

**Lemma A.4** (Appearance of a double root, two times) Assume $\frac{A^2d}{4} \le \mu < \frac{8A^2d}{27}$, then there are $0 \le \delta_1 < \delta_2$ such that the following holds.

(i) If $\delta \neq \delta_1$ and $\delta \neq \delta_2$, the three roots are simple and as $\delta$ encounters $\delta_i$ ($i = 1$ or $2$), two of the three simple roots merge in a double root $\lambda_i$.
<figure><img src="image_13.png" /><figcaption>**Figure X** — Multiplicity of the roots of $P_\delta$ for fixed $\nu$, $d$. Fix the last parameter $\mu > 0$, then $\delta$ browses $[0, +\infty)$.</figcaption></figure>

**Lemma A.2** (Only simple roots) Assume $\mu > \frac{8A^2d}{27}$ is fixed and set $\mathcal{E} := \mathbb{R}_+$.

(i) For all $\delta \in \mathcal{E}$, the three roots are simple, one is real and the two others are complex conjugated.

(ii) There is $\varepsilon > 0$ such that $\inf_{\delta \in \mathcal{E}} (|\alpha - \beta|, |\alpha - \gamma|, |\beta - \gamma|) > \varepsilon$.

**Lemma A.3** (Appearance of a triple root) Assume $\mu = \frac{8A^2d}{27}$ and set $\delta_0 := \frac{A^2d}{27}$.

(i) If $\delta \neq \delta_0$, the three roots are simple, one is real and the two others are complex conjugated, and as $\delta$ encounters $\delta_0$, the three simple roots merge in a triple root $\lambda_0 = -\frac{A\sqrt{d}}{3}$.

Without loss of generality, we assume that $\alpha$ is the real root for all $\delta \ge 0$.

(ii) Letting $\varepsilon_R := \mathfrak{Re}(\beta) - \lambda_0$ and $\varepsilon_I := \mathfrak{Im}(\beta)$, one has

$$ \alpha = \lambda_0 - 2\varepsilon_R, \quad \beta = \lambda_0 + \varepsilon_R + i\varepsilon_I, \quad \gamma = \lambda_0 + \varepsilon_R - i\varepsilon_I, $$

and $\varepsilon_R \stackrel{\delta \to \delta_0}{=} \mathcal{O}(\varepsilon_I)$.

Define $\eta := \frac{\delta_0}{2}$ and $\mathcal{E} := \mathbb{R}_+ \setminus (\delta_0 - \eta, \delta_0 + \eta)$, then there is $\varepsilon > 0$ such that

(iii) $\inf_{\delta \in \mathcal{E}} (|\alpha - \beta|, |\alpha - \gamma|, |\beta - \gamma|) > \varepsilon$,

(iv) $\inf_{\delta \in \mathcal{E}^c} \inf_{\lambda \in [\alpha, \beta]} |\lambda| > \varepsilon$.

**Lemma A.4** (Appearance of a double root, two times) Assume $\frac{A^2d}{4} \le \mu < \frac{8A^2d}{27}$, then there are $0 \le \delta_1 < \delta_2$ such that the following holds.

(i) If $\delta \neq \delta_1$ and $\delta \neq \delta_2$, the three roots are simple and as $\delta$ encounters $\delta_i$ ($i = 1$ or $2$), two of the three simple roots merge in a double root $\lambda_i$.
Without loss of generality, we assume that, for $i = 1$ or $2$, the double root $\lambda_i$ derives from the merging of $\alpha$ and $\beta$. Define $\eta := \frac{\delta_2 - \delta_1}{2}$ and $\mathcal{E} := \mathbb{R}_+ \setminus [(\delta_1 - \eta, \delta_1 + \eta) \cup (\delta_2 - \eta, \delta_2 + \eta)]$, then there is $\varepsilon > 0$ such that

$$ (ii) \inf_{\delta \in \mathcal{E}} (|\alpha - \beta|, |\alpha - \gamma|, |\beta - \gamma|) > \varepsilon, $$

$$ (iii) \inf_{\delta \in \mathcal{E}^c} \inf_{\lambda \in [\alpha, \beta]} (|\lambda|, |\lambda - \gamma|) > \varepsilon. $$

**Lemma A.5** (Appearance of a double root, one time) Assume $0 < \mu < \frac{A^2d}{4}$, then there is $\delta_2 > 0$ such that the following holds.

(i) If $\delta \neq \delta_2$, the three roots are simple and as $\delta$ encounters $\delta_2$, two of the three simple roots merge in a double root $\lambda_2$.

Without loss of generality, we assume that the double root $\lambda_2$ derives from the merging of $\alpha$ and $\beta$. Define $\eta := \frac{\delta_2}{2}$ and $\mathcal{E} := \mathbb{R}_+ \setminus (\delta_2 - \eta, \delta_2 + \eta)$, then there is $\varepsilon > 0$ such that

$$ (ii) \inf_{\delta \in \mathcal{E}} (|\alpha - \beta|, |\alpha - \gamma|, |\beta - \gamma|) > \varepsilon, $$

$$ (iii) \inf_{\delta \in \mathcal{E}^c} \inf_{\lambda \in [\alpha, \beta]} (|\lambda|, |\lambda - \gamma|) > \varepsilon. $$

## A.3 Fourier, Laplace and Fourier/Laplace transforms

The transforms of a function $w$ depending on $t > 0$ and $x \in \mathbb{R}^{N-1}$ we use are

* the $x$-Fourier one (said “hat $w$”):

$$ \hat{w}(t, \xi) = \mathcal{F}[w(t, \cdot)](\xi) := \int_{\mathbb{R}^{N-1}} w(t, z)e^{-i\xi \cdot z} dz; $$

* the $t$-Laplace one (said “frown $w$”):

$$ \hat{w}(s, x) = \mathcal{L}[w(\cdot, x)](s) := \int_{0}^{+\infty} w(\tau, x)e^{-s\tau} d\tau; $$

* the $x$-Fourier/$t$-Laplace one (said “frown-hat $w$”):

$$ \hat{w}(s, \xi) = \mathcal{FL}[w(\cdot, \cdot)](s, \xi) := \int_{\mathbb{R}^{N-1}} \int_{0}^{+\infty} w(\tau, z)e^{-(s\tau + i\xi \cdot z)} d\tau dz. $$

The following Fourier and Laplace transforms are well-known — see, e.g., [30, page 250, lines 84 and 88] for (iii) and (iv).

**Lemma A.6** (Useful $x$-Fourier and $t$-Laplace transforms) For all $a \in \mathbb{R}$ and all $b \in \mathbb{C}$,

$$ (i) \mathcal{F} \left[ x \mapsto \frac{e^{-\frac{\|x\|^2}{4a}}}{(4\pi a)^{\frac{N-1}{2}}} \right] (\xi) = e^{-a\|\xi\|^2}, \quad a > 0. $$

$$ (ii) \mathcal{L}[t \mapsto e^{-at}] (s) = \frac{1}{s+a}, \quad s > -a. $$
Without loss of generality, we assume that, for $i = 1$ or $2$, the double root $\lambda_i$ derives from the merging of $\alpha$ and $\beta$. Define $\eta := \frac{\delta_2 - \delta_1}{2}$ and $\mathcal{E} := \mathbb{R}_+ \setminus [(\delta_1 - \eta, \delta_1 + \eta) \cup (\delta_2 - \eta, \delta_2 + \eta)]$, then there is $\varepsilon > 0$ such that

$$ (ii) \inf_{\delta \in \mathcal{E}} (|\alpha - \beta|, |\alpha - \gamma|, |\beta - \gamma|) > \varepsilon, $$

$$ (iii) \inf_{\delta \in \mathcal{E}^c} \inf_{\lambda \in [\alpha, \beta]} (|\lambda|, |\lambda - \gamma|) > \varepsilon. $$

**Lemma A.5 (Appearance of a double root, one time)** Assume $0 < \mu < \frac{A^2d}{4}$, then there is $\delta_2 > 0$ such that the following holds.

(i) If $\delta \neq \delta_2$, the three roots are simple and as $\delta$ encounters $\delta_2$, two of the three simple roots merge in a double root $\lambda_2$.

Without loss of generality, we assume that the double root $\lambda_2$ derives from the merging of $\alpha$ and $\beta$. Define $\eta := \frac{\delta_2}{2}$ and $\mathcal{E} := \mathbb{R}_+ \setminus (\delta_2 - \eta, \delta_2 + \eta)$, then there is $\varepsilon > 0$ such that

$$ (ii) \inf_{\delta \in \mathcal{E}} (|\alpha - \beta|, |\alpha - \gamma|, |\beta - \gamma|) > \varepsilon, $$

$$ (iii) \inf_{\delta \in \mathcal{E}^c} \inf_{\lambda \in [\alpha, \beta]} (|\lambda|, |\lambda - \gamma|) > \varepsilon. $$

## A.3 Fourier, Laplace and Fourier/Laplace transforms

The transforms of a function $w$ depending on $t > 0$ and $x \in \mathbb{R}^{N-1}$ we use are

* the $x$-Fourier one (said “hat $w$”):

$$ \hat{w}(t, \xi) = \mathcal{F}[w(t, \cdot)](\xi) := \int_{\mathbb{R}^{N-1}} w(t, z)e^{-i\xi \cdot z} dz; $$

* the $t$-Laplace one (said “frown $w$”):

$$ \hat{w}(s, x) = \mathcal{L}[w(\cdot, x)](s) := \int_{0}^{+\infty} w(\tau, x)e^{-s\tau} d\tau; $$

* the $x$-Fourier/$t$-Laplace one (said “frown-hat $w$”):

$$ \hat{w}(s, \xi) = \mathcal{FL}[w(\cdot, \cdot)](s, \xi) := \int_{\mathbb{R}^{N-1}} \int_{0}^{+\infty} w(\tau, z)e^{-(s\tau + i\xi \cdot z)} d\tau dz. $$

The following Fourier and Laplace transforms are well-known — see, e.g., [30, page 250, lines 84 and 88] for (iii) and (iv).

**Lemma A.6 (Useful $x$-Fourier and $t$-Laplace transforms)** For all $a \in \mathbb{R}$ and all $b \in \mathbb{C}$,

$$ (i) \mathcal{F} \left[ x \mapsto \frac{e^{-\frac{\|x\|^2}{4a}}}{(4\pi a)^{\frac{N-1}{2}}} \right] (\xi) = e^{-a\|\xi\|^2}, \quad a > 0. $$

$$ (ii) \mathcal{L}[t \mapsto e^{-at}] (s) = \frac{1}{a+s}, \quad s > -a. $$
$$ (iii) \mathcal{L} \left[ t \mapsto \frac{e^{-\frac{a^2}{4t}}}{\sqrt{\pi t}} \right] (s) = \frac{e^{-a\sqrt{s}}}{\sqrt{s}}, \quad s > 0. $$

$$ (iv) \mathcal{L} \left[ t \mapsto \frac{\text{Erfc}}{\Gamma} \left( \frac{a+2bt}{2\sqrt{t}} \right) e^{-\frac{a^2}{4t}} \right] (s) = \frac{e^{-a\sqrt{s}}}{\sqrt{s}(\sqrt{s}+b)}, \quad s > \mathfrak{Re}(b)^2. $$

$$ (v) \mathcal{L} \left[ t \mapsto e^{-\frac{a^2}{4t}} \left( \frac{1}{\sqrt{\pi t}} + b \frac{\text{Erfc}}{\Gamma} \left( \frac{a-2bt}{2\sqrt{t}} \right) \right) \right] (s) = \frac{e^{-a\sqrt{s}}}{\sqrt{s}-b}, \quad s > \mathfrak{Re}(b)^2. $$

where $\Gamma(\ell) := e^{-\ell^2}$ and the definition of Erfc being recalled in Appendix A.1.

**Lemma A.7** (Properties on $x$-Fourier and $t$-Laplace transforms) For $f,g: \mathbb{R} \to \mathbb{C}$, all $a \in \mathbb{R}$ and $s > 0$,

$$ (i) f(x) = \frac{1}{(2\pi)^{N-1}} \int_{\mathbb{R}^{N-1}} \mathcal{F}[f](\xi) e^{i\xi \cdot x} d\xi, $$

$$ (ii) \mathcal{F}[f] \times \mathcal{F}[g] = \mathcal{F}[f * g] = \mathcal{F}[x \mapsto \int_{\mathbb{R}^{N-1}} f(x-z)g(z)dz], $$

$$ (iii) \mathcal{L}[f] \times \mathcal{L}[g] = \mathcal{L}\left[t \mapsto \int_{0}^{t} f(t-\tau)g(\tau)d\tau\right], $$

$$ (iv) \mathcal{L}[t \mapsto e^{-at}f(t)](s) = \mathcal{L}[f](s+a). $$

**Acknowledgement.** Matthieu Alfaro is supported by the *région Normandie* project BIOMANORMAN 21E04343 and the ANR project DEEV ANR-20-CE40-0011-01. Samuel Tréton would like to acknowledge the *région Normandie* for the financial support of his PhD.

## References

[1] E. AFFILI, *A Fisher-KPP model with a fast diffusion line in periodic media*, arXiv:2009.14760 [math], (2020).

[2] M. ALFARO, *Fujita blow up phenomena and hair trigger effect: The role of dispersal tails*, Annales de l'Institut Henri Poincaré C, Analyse non linéaire, 34 (2017), pp. 1309-1327.

[3] M. ALFARO AND O. KAVIAN, *Blow-up phenomena for positive solutions of semilinear diffusion equations in a half-space: The influence of the dispersion kernel*, arXiv:2004.09102 [math], (2020).

[4] W. C. ALLEE, *Animal Aggregations: A Study in General Sociology*, 1931.

[5] H. BERESTYCKI, A.-C. COULON, J.-M. ROQUEJOFFRE, AND L. ROSSI, *Speed-up of reaction-diffusion fronts by a line of fast diffusion*, Séminaire Laurent Schwartz — EDP et applications, (2013), pp. 1-25.

[6] ——, *The effect of a line with non-local diffusion on Fisher-KPP propagation*, arXiv:1506.09046 [math], (2015).

[7] H. BERESTYCKI, R. DUCASSE, AND L. ROSSI, *Generalized principal eigenvalues for heterogeneous road-field systems*, arXiv:1810.13180 [math], (2018).

[8] ——, *Influence of a road on a population in an ecological niche facing climate change*, arXiv:1903.02221 [math], (2019).
$$ (iii) \mathcal{L} \left[ t \mapsto \frac{e^{-\frac{a^2}{4t}}}{\sqrt{\pi t}} \right] (s) = \frac{e^{-a\sqrt{s}}}{\sqrt{s}}, \quad s > 0. $$

$$ (iv) \mathcal{L} \left[ t \mapsto \frac{\text{Erfc}}{\Gamma} \left( \frac{a+2bt}{2\sqrt{t}} \right) e^{-\frac{a^2}{4t}} \right] (s) = \frac{e^{-a\sqrt{s}}}{\sqrt{s}(\sqrt{s}+b)}, \quad s > \text{Re}(b)^2. $$

$$ (v) \mathcal{L} \left[ t \mapsto e^{-\frac{a^2}{4t}} \left( \frac{1}{\sqrt{\pi t}} + b \frac{\text{Erfc}}{\Gamma} \left( \frac{a-2bt}{2\sqrt{t}} \right) \right) \right] (s) = \frac{e^{-a\sqrt{s}}}{\sqrt{s}-b}, \quad s > \text{Re}(b)^2. $$

where $\Gamma(\ell) := e^{-\ell^2}$ and the definition of Erfc being recalled in Appendix A.1.

**Lemma A.7** (Properties on $x$-Fourier and $t$-Laplace transforms) For $f, g: \mathbb{R} \to \mathbb{C}$, all $a \in \mathbb{R}$ and $s > 0$,

$$ (i) f(x) = \frac{1}{(2\pi)^{N-1}} \int_{\mathbb{R}^{N-1}} \mathcal{F}[f](\xi) e^{i\xi \cdot x} d\xi, $$

$$ (ii) \mathcal{F}[f] \times \mathcal{F}[g] = \mathcal{F}[f * g] = \mathcal{F}[x \mapsto \int_{\mathbb{R}^{N-1}} f(x-z)g(z)dz], $$

$$ (iii) \mathcal{L}[f] \times \mathcal{L}[g] = \mathcal{L}\left[t \mapsto \int_{0}^{t} f(t-\tau)g(\tau)d\tau\right], $$

$$ (iv) \mathcal{L}\left[t \mapsto e^{-at}f(t)\right](s) = \mathcal{L}[f](s+a). $$

**Acknowledgement.** Matthieu Alfaro is supported by the **région Normandie** project BIOMANORMAN 21E04343 and the ANR project DEEV ANR-20-CE40-0011-01. Samuel Tréton would like to acknowledge the **région Normandie** for the financial support of his PhD. The authors are grateful to the anonymous referee whose precise comments have improved the presentation of the results.

## References

[1] E. AFFILI, *A Fisher-KPP model with a fast diffusion line in periodic media*, arXiv:2009.14760 [math], (2020).

[2] M. ALFARO, *Fujita blow up phenomena and hair trigger effect: The role of dispersal tails*, Annales de l'Institut Henri Poincaré C, Analyse non linéaire, 34 (2017), pp. 1309-1327.

[3] M. ALFARO AND O. KAVIAN, *Blow-up phenomena for positive solutions of semilinear diffusion equations in a half-space: the influence of the dispersion kernel*, Ann. Fac. Sci. Toulouse Math. (6), 31 (2022), pp. 1259-1286.

[4] W. C. ALLEE, *Animal Aggregations: A Study in General Sociology*, 1931.

[5] H. BERESTYCKI, A.-C. COULON, J.-M. ROQUEJOFFRE, AND L. ROSSI, *Speed-up of reaction-diffusion fronts by a line of fast diffusion*, Séminaire Laurent Schwartz — EDP et applications, (2013), pp. 1-25.

[6] ——, *The effect of a line with nonlocal diffusion on Fisher-KPP propagation*, Math. Models Methods Appl. Sci., 25 (2015), pp. 2519-2562.

[7] H. BERESTYCKI, R. DUCASSE, AND L. ROSSI, *Generalized principal eigenvalues for heterogeneous road-field systems*, Commun. Contemp. Math., 22 (2020), pp. 1950013, 35.
[9] H. BERESTYCKI, J.-M. ROQUEJOFFRE, AND L. ROSSI, *Fisher-KPP propagation in the presence of a line: Further effects*, Nonlinearity, 26 (2013), pp. 2623-2640.

[10] ——, *The influence of a line with fast diffusion on Fisher-KPP propagation*, Journal of Mathematical Biology, 66 (2013), pp. 743-766.

[11] ——, *The shape of expansion induced by a line with fast diffusion in Fisher-KPP equations*, Communications in Mathematical Physics, 343 (2016), pp. 207-232.

[12] ——, *Travelling waves, spreading and extinction for Fisher-KPP propagation driven by a line with fast diffusion*, Nonlinear Analysis, 137 (2016), pp. 171-189.

[13] B. BOGOSEL, T. GILETTI, AND A. TELLINI, *Propagation for KPP bulk-surface systems in a general cylindrical domain*, Nonlinear Analysis, 213 (2021), p. 112528.

[14] A. A. M. CUYT, V. PETERSEN, B. VERDONK, H. WAADELAND, AND W. B. JONES, *Handbook of Continued Fractions for Special Functions*, Springer Science & Business Media, Apr. 2008.

[15] R. DUCASSE, *Influence of the geometry on a field-road model: The case of a conical field*, Journal of the London Mathematical Society, 97 (2018), pp. 441-469.

[16] H. FUJITA, *On the blowing up of solutions of the Cauchy problem for $u_t = \Delta u + u^{1+\alpha}$*, (1966).

[17] M. GATTO, E. BERTUZZO, L. MARI, S. MICCOLI, L. CARRARO, R. CASAGRANDI, AND A. RINALDO, *Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures*, Proceedings of the National Academy of Sciences, 117 (2020), pp. 10484-10491.

[18] T. GILETTI, L. MONSAINGEON, AND M. ZHOU, *A KPP road-field system with spatially periodic exchange terms*, Nonlinear Analysis, 128 (2015), pp. 273-302.

[19] P. HENRICI, *Applied and computational complex analysis. vol.2: Special functions, integral transforms, asymptotics, continued fractions*, (1977).

[20] H. W. MCKENZIE, E. H. MERRILL, R. J. SPITERI, AND M. A. LEWIS, *How linear features alter predator movement and the functional response*, Interface Focus, 2 (2012), pp. 205-216.

[21] J. D. MURRAY, *Mathematical Biology. Vol. 1: An Introduction*, vol. 17 of Interdisciplinary Applied Mathematics, Springer, 3rd ed., 2002.

[22] ——, *Mathematical Biology. Vol. 2: Spatial Models and Biomedical Applications*, vol. 18 of Interdisciplinary Applied Mathematics, Springer, 3rd ed., 2003.

[23] A. PAUTHIER, *The influence of a line with fast diffusion and nonlocal exchange terms on Fisher-KPP propagation*, arXiv:1408.1584 [math], (2015).

[24] ——, *Road-field reaction-diffusion system: A new threshold for long range exchanges*, arXiv:1504.05437 [math], (2015).
[8] ——, *Influence of a road on a population in an ecological niche facing climate change*, J. Math. Biol., 81 (2020), pp. 1059–1097.

[9] H. BERESTYCKI, J.-M. ROQUEJOFFRE, AND L. ROSSI, *Fisher-KPP propagation in the presence of a line: Further effects*, Nonlinearity, 26 (2013), pp. 2623–2640.

[10] ——, *The influence of a line with fast diffusion on Fisher-KPP propagation*, Journal of Mathematical Biology, 66 (2013), pp. 743–766.

[11] ——, *The shape of expansion induced by a line with fast diffusion in Fisher-KPP equations*, Communications in Mathematical Physics, 343 (2016), pp. 207–232.

[12] ——, *Travelling waves, spreading and extinction for Fisher-KPP propagation driven by a line with fast diffusion*, Nonlinear Analysis, 137 (2016), pp. 171–189.

[13] B. BOGOSEL, T. GILETTI, AND A. TELLINI, *Propagation for KPP bulk-surface systems in a general cylindrical domain*, Nonlinear Analysis, 213 (2021), p. 112528.

[14] A. A. M. CUYT, V. PETERSEN, B. VERDONK, H. WAADELAND, AND W. B. JONES, *Handbook of Continued Fractions for Special Functions*, Springer Science & Business Media, Apr. 2008.

[15] R. DUCASSE, *Influence of the geometry on a field-road model: The case of a conical field*, Journal of the London Mathematical Society, 97 (2018), pp. 441–469.

[16] H. FUJITA, *On the blowing up of solutions of the Cauchy problem for $u_t = \Delta u + u^{1+\alpha}$*, J. Fac. Sci. Univ. Tokyo Sect. I, 13 (1966), pp. 109–124 (1966).

[17] M. GATTO, E. BERTUZZO, L. MARI, S. MICCOLI, L. CARRARO, R. CASAGRANDI, AND A. RINALDO, *Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures*, Proceedings of the National Academy of Sciences, 117 (2020), pp. 10484–10491.

[18] T. GILETTI, L. MONSAINGEON, AND M. ZHOU, *A KPP road-field system with spatially periodic exchange terms*, Nonlinear Analysis, 128 (2015), pp. 273–302.

[19] P. HENRICI, *Applied and computational complex analysis. vol.2: Special functions, integral transforms, asymptotics, continued fractions*, (1977).

[20] H. W. MCKENZIE, E. H. MERRILL, R. J. SPITERI, AND M. A. LEWIS, *How linear features alter predator movement and the functional response*, Interface Focus, 2 (2012), pp. 205–216.

[21] J. D. MURRAY, *Mathematical Biology. Vol. 1: An Introduction*, vol. 17 of Interdisciplinary Applied Mathematics, Springer, 3rd ed., 2002.

[22] ——, *Mathematical Biology. Vol. 2: Spatial Models and Biomedical Applications*, vol. 18 of Interdisciplinary Applied Mathematics, Springer, 3rd ed., 2003.

[23] A. PAUTHIER, *The influence of a line with fast diffusion and nonlocal exchange terms on Fisher-KPP propagation*, arXiv:1408.1584 [math], (2015).
[25] ——, *Uniform dynamics for Fisher-KPP propagation driven by a line of fast diffusion under a singular limit*, Nonlinearity, 28 (2015), pp. 3891–3920.

[26] C. ROBINET, C.-E. IMBERT, J. ROUSSELET, D. SAUVARD, J. GARCIA, F. GOUSSARD, AND A. ROQUES, *Human-mediated long-distance jumps of the pine processionary moth in Europe*, Biological Invasions, 14 (2012), pp. 1557–1569.

[27] L. ROSSI, A. TELLINI, AND E. VALDINOCI, *The effect on Fisher-KPP propagation in a cylinder with fast diffusion on the boundary*, SIAM Journal on Mathematical Analysis, 49 (2017), pp. 4595–4624.

[28] B. V. SCHMID, U. BÜNTGEN, W. R. EASTERDAY, C. GINZLER, L. WALLØE, B. BRAMANTI, AND N. C. STENSETH, *Climate-driven introduction of the Black Death and successive plague reintroductions into Europe*, Proceedings of the National Academy of Sciences, 112 (2015), pp. 3020–3025.

[29] J. SMOLLER, *Shock Waves and Reaction-Diffusion Equations*, Springer Science & Business Media, 1994.

[30] M. R. SPIEGEL, *Theory and Problems of Laplace Transforms*, Schaum Publishing Co., New York, 1965.

[31] A. TELLINI, *Propagation speed in a strip bounded by a line with different diffusion*, Journal of Differential Equations, 260 (2016), pp. 5956–5986.

[32] A. TYCHONOFF, *Théoremes d'unicité pour l'équation de la chaleur*, 42 (1935), pp. 199–215.

[33] M. ZHANG, *Spreading speeds and pulsating fronts for a field-road model in a spatially periodic habitat*, Journal of Differential Equations, 304 (2021), pp. 191–228.
[24] ——, *Road-field reaction-diffusion system: A new threshold for long range exchanges*, arXiv:1504.05437 [math], (2015).

[25] ——, *Uniform dynamics for Fisher-KPP propagation driven by a line of fast diffusion under a singular limit*, Nonlinearity, 28 (2015), pp. 3891–3920.

[26] C. ROBINET, C.-E. IMBERT, J. ROUSSELET, D. SAUVARD, J. GARCIA, F. GOUSSARD, AND A. ROQUES, *Human-mediated long-distance jumps of the pine processionary moth in Europe*, Biological Invasions, 14 (2012), pp. 1557–1569.

[27] L. ROSSI, A. TELLINI, AND E. VALDINOCI, *The effect on Fisher-KPP propagation in a cylinder with fast diffusion on the boundary*, SIAM Journal on Mathematical Analysis, 49 (2017), pp. 4595–4624.

[28] B. V. SCHMID, U. BÜNTGEN, W. R. EASTERDAY, C. GINZLER, L. WALLØE, B. BRAMANTI, AND N. C. STENSETH, *Climate-driven introduction of the Black Death and successive plague reintroductions into Europe*, Proceedings of the National Academy of Sciences, 112 (2015), pp. 3020–3025.

[29] J. SMOLLER, *Shock Waves and Reaction-Diffusion Equations*, Springer Science & Business Media, 1994.

[30] M. R. SPIEGEL, *Theory and Problems of Laplace Transforms*, Schaum Publishing Co., New York, 1965.

[31] A. TELLINI, *Propagation speed in a strip bounded by a line with different diffusion*, Journal of Differential Equations, 260 (2016), pp. 5956–5986.

[32] A. TYCHONOFF, *Théoremes d'unicité pour l'équation de la chaleur*, 42 (1935), pp. 199–215.

[33] M. ZHANG, *Spreading speeds and pulsating fronts for a field-road model in a spatially periodic habitat*, Journal of Differential Equations, 304 (2021), pp. 191–228.