# Gödel-type universes in $f(T)$ gravity

Di Liu, Puxun Wu and Hongwei Yu  Center of Nonlinear Science and Department of
Physics, Ningbo University, Ningbo, Zhejiang, 315211 China

###### Abstract

The issue of causality in $f(T)$ gravity is investigated by examining the
possibility of existence of the closed timelike curves in the Gödel-type
metric. By assuming a perfect fluid as the matter source, we find that the
fluid must have an equation of state parameter greater than minus one in order
to allow the Gödel solutions to exist, and furthermore the critical radius
$r_{c}$, beyond which the causality is broken down, is finite and it depends
on both matter and gravity. Remarkably, for certain $f(T)$ models, the perfect
fluid that allows the Gödel-type solutions can even be normal matter, such as
pressureless matter or radiation. However, if the matter source is a special
scalar field rather than a perfect fluid, then $r_{c}\rightarrow\infty$ and
the causality violation is thus avoided.

###### pacs:

04.50.Kd, 04.20.Jb, 98.80.Jk

##  I Introduction

General relativity (GR) is established in the framework of the Levi-Civita
connection, therefore there is only curvature rather than torsion in the
spacetime. On the other hand, one can also introduce other connections, such
as the Weitzenböck connection, into the same spacetime where only torsion is
reserved. Thus, there is no such a thing as curvature or torsion of spacetime,
but only curvature or torsion of connection. Basing on the Weitzenböck
connection, Einstein Einstein  introduced firstly the Teleparallel Gravity
(TG) in his endeavor to unify gravity and electromagnetism with the
introduction of a tetrad field. TG can, as is well known, show up as a theory
completely equivalent to GR since the difference between their actions (the
actions of TG and GR are the torsion scalar $T$ and Ricci scalar $R$,
respectively) is just a derivative term FNGtnb ; FNGtn1 ; FNGtn2 ; FNGtn3 ;
FNGtn4 ; FNGtne .

Recently, a modification of TG, called $f(T)$ theory Bengochea2009 ;
Ferraro2007 ; Linder ; Zheng2011 ; Ferraro2008 ; Ferraro2011 ; pwhy2011 ;
Wu2011 ; Bamba2011 ; Wu2010a ; Ben2011ab ; Wu2010b ; Zhang2011bb ; FTbe ; FT1
; FT2 ; FT3 ; FT4 ; FT5 ; FT6 ; FT7 ; FT8 ; FT9 ; FT10 ; FT11 ; FT12 ; FT13 ;
FT14 ; FT15 ; FT16 ; FT17 ; FT18 ; FT19 ; FT20 ; FT21 ; FT22 ; FT23 ; FT24 ;
FT25 ; FT26 ; FT27 ; FT28 ; FT29 ; FT30 ; FT31 ; FT32 ; FT33 ; FT34 ; FT35 ;
FT36 ; FT37 ; FT38 ; FTed ; Ferraro2011a ; Li2011aa ; Li2011bb ; Li2011b ;
LiM2011 ; Miao2011 , has spurred an increasing deal of attention, as it can
explain the present accelerated cosmic expansion discovered from observations
(the Type Ia supernova R98 ; P99 , the cosmic microwave background radiation
Spa ; Spb , and the large scale structure T2004 ; E2005 , etc.) without the
need of dark energy. $f(T)$ theory is obtained by generalizing the action $T$
of TG to an arbitrary function $f$ of $T$, which is very analogous to $f(R)$
theory (see FeNoj08 ; FeSot10 ; Felice2010 ; FeNoj11 ; FebCli  for recent
review) where the action $R$ of GR is generalized to be $f(R)$. An advantage
of $f(T)$ theory is that its field equation is only second order, while in
$f(R)$ gravity it is forth order.

It has been found that $f(T)$ theory can give an inflation without an inflaton
Ferraro2007 ; Ferraro2008 , avoid the big bang singularity problem in the
standard cosmological model Ferraro2011 , realize the crossing of phantom
divide line for the effective equation of state Wu2011 ; Bamba2011 , and yield
an usual early cosmic evolution Wu2010b ; Zhang2011bb . But, at the same, this
theory lacks the local Lorentz invariance Li2011aa ; Li2011bb , and this
results in the appearance of extra degrees of freedom LiM2011 , the broken
down of the first law of black hole thermodynamic Miao2011 , and the problem
in cosmic large scale structure Li2011b .

In this paper, we plan to study the causality issue of $f(T)$ theory by
examining the possibility of existence of the closed timelike curves in the
Gödel spacetime Godel . The Gödel metric is the first cosmological solution
with rotating matter to the Einstein equation in GR. Since the Gödel solution
is very convenient for studying whether the closed timelike curves exist, it
has been used widely to test the causality issue. For example, Gödel found
that the closed timelike solution cannot be excluded in GR, assuming a
cosmological constant or a perfect fluid with its pressure equal to the energy
density. Gödel’s work has been generalized to include other matter sources,
such as, the vector field Sombe ; SomRe79 ; SomRay80 , scalar field Hiscockbe
; HisCha ; HisPan , spinor field Villalbabe ; VilPim ; VilKre ; VilLea ;
VilHered ; Reboucas1983  and tachyon field Reboucas . In addition, the Gödel-
type universes mjrteibe ; mjrtei1 ; mjrteied ; Reboucas1983  have also been
studied in the framework of other theories of gravitation, such as TG TGGod ,
$f(R)$ gravity RebClif05e ; Reb09b1 ; Reb10ed  and string-inspired
gravitational theory stri ; barrow1998 .

Here, assuming that the matter source is the perfect fluid or a scalar field,
we aim to find out the condition for non-violation of causality in $f(T)$
gravity. The paper is organized as follows. We give, in Sec. II, a brief
review of $f(T)$ theory and the vierbein of a general cylindrical symmetry
metric in Sec.III. The Gödel-type universe in $f(T)$ theory is discussed in
Sec. IV. With an assumption of different matter sources, we investigate the
issue of causality in Sec. V. Finally, we present our conclusions in Sec. VI.

##  II  $f(T)$ gravity

In this section, we give a brief view of $f(T)$ gravity. We use the Greek
alphabet ($\mu$, $\nu$, $\cdots$= 0, 1, 2, 3) to denote tensor indices, that
is, indices related to spacetime, and middle part of the Latin alphabet ($i$,
$j$, $\cdots$= 0, 1, 2, 3) to denote tangent space (local Lorentzian) indices.
TG, instead of using the metric tensor, uses tetrad, $e_{\mu}^{i}$ or
$e_{i}^{\mu}$ (frame or coframe), as the dynamical object. The relation
between frame and coframe is

$e_{i}^{\mu}e_{\mu}^{j}=\delta_{i}^{j}\;,\qquad
e_{i}^{\mu}e_{\nu}^{i}=\delta_{\nu}^{\mu},$ (1)

and the relation between tetrad and metric tensor is

$g_{\mu\nu}=e_{\mu}^{i}e_{\nu}^{j}\eta_{ij}\;,\qquad\eta_{ij}=e_{i}^{\mu}e_{j}^{\nu}g_{\mu\nu}\;,$
(2)

where $\eta_{ij}=diag(1,-1,-1,-1)$ is the Minkowski metric.

Different from GR, the Weitzenböck connection is used in TG

${\Gamma}_{\mu\nu}^{\lambda}=e_{i}^{\lambda}\partial_{\nu}e_{\mu}^{i}=-e_{\mu}^{i}\partial_{\nu}e_{i}^{\lambda}\;.$
(3)

As a result, the convariant derivative, denoted by $D_{\mu}$, satisfies:

$D_{\mu}e_{\nu}^{i}=\partial_{\mu}e_{\nu}^{i}-\Gamma_{\nu\mu}^{\lambda}e_{\lambda}^{i}=0\;.$
(4)

To describe the difference between Weitzenböck and Levi-Civita connections, a
contorsion tensor $K_{\;\;\mu\nu}^{\rho}$ needs to be introduced:

$K_{\;\;\mu\nu}^{\rho}\equiv\Gamma_{\mu\nu}^{\rho}-\overset{\circ}{\Gamma}{}_{\mu\nu}^{\rho}=\frac{1}{2}(T_{\mu}{}^{\rho}{}_{\nu}+T_{\nu}{}^{\rho}{}_{\mu}-T_{\;\;\mu\nu}^{\rho})\;.$
(5)

Here $T_{\;\;\mu\nu}^{\rho}$ is the torsion tensor

$T_{\;\;\mu\nu}^{\rho}={\Gamma}_{\nu\mu}^{\rho}-{\Gamma}_{\mu\nu}^{\rho}=e_{i}^{\rho}(\partial_{\mu}e_{\nu}^{i}-\partial_{\nu}e_{\mu}^{i})\;,$
(6)

and $\overset{\circ}{\Gamma}{}_{\mu\nu}^{\rho}$ denotes the Levi-Civita
connection

$\overset{\circ}{\Gamma}{}_{\mu\nu}^{\rho}=\frac{1}{2}g^{\rho\sigma}(\partial_{\mu}g_{\sigma\nu}+\partial_{\nu}g_{\sigma\mu}-\partial_{\sigma}g_{\mu\nu}).$
(7)

By defining the super-potential $S_{\sigma}^{\;\;\mu\nu}$

$S_{\sigma}^{\;\;\mu\nu}\equiv
K_{\;\;\;\;\sigma}^{\mu\nu}+\delta_{\sigma}^{\mu}T_{\;\;\;\;\;\alpha}^{\alpha\nu}-\delta_{\sigma}^{\nu}T_{\;\;\;\;\;\alpha}^{\alpha\mu}\;,$
(8)

we obtain the torsion scalar $T$

$T\equiv\frac{1}{2}S_{\sigma}^{\;\;\mu\nu}T_{\;\;\mu\nu}^{\sigma}=\frac{1}{4}T^{\alpha\mu\nu}T_{\alpha\mu\nu}+\frac{1}{2}T^{\alpha\mu\nu}T_{\nu\mu\alpha}-T_{\alpha\mu}^{\;\;\;\;\alpha}T_{\;\;\;\;\;\nu}^{\nu\mu}\;.$
(9)

In TG, the Lagrangian density is given by:

$L_{T}=\frac{eT}{2\kappa^{2}}\;,$ (10)

where, $e=\det(e_{\mu}^{i})=\sqrt{-g}\;,\kappa^{2}{\equiv}8\pi G$.
Generalizing $T$ to be an arbitrary function $f$ of $T$ in the above
expression, we obtain the Lagrangian density of $f(T)$ theory

$L_{T}=\frac{ef(T)}{2\kappa^{2}}\;.$ (11)

Adding a matter Lagrangian density $L_{M}$ to Eq. (11), and varying the action
with respect to the vierbein, one finds the following field equation of $f(T)$
theory:

$\displaystyle[e^{-1}\partial_{\mu}(ee^{\rho}_{i}S^{\;\;\nu\mu}_{\rho})-e^{\lambda}_{i}S^{\rho\mu\nu}T_{\rho\mu\lambda}]f_{T}(T)+e^{\rho}_{i}S^{\;\;\nu\mu}_{\rho}\partial_{\mu}(T)f_{TT}(T)$
(12)
$\displaystyle+\frac{1}{2}e^{\nu}_{i}f(T)=\kappa^{2}e^{\rho}_{i}\overset{em}{T}{}^{\;\;\nu}_{\rho}.$

Here $f_{T}=df(T)/dT$, $f_{TT}=d^{2}f(T)/dT^{2}$, and
$\overset{em}{T}{}^{\nu}_{\rho}$ is the matter energy-momentum tensor. In a
coordinate system, this field equation can be rewritten as

$\displaystyle
A_{\mu\nu}f_{T}(T)+S^{\;\;\;\;\;\sigma}_{\nu\mu}(\nabla_{\sigma}T)f_{TT}(T)+\frac{1}{2}g_{\mu\nu}f(T)=\kappa^{2}\overset{em}{T}{}_{\mu\nu}\;,$
(13)

where

$\displaystyle
A_{\mu\nu}=g_{\sigma\mu}e^{i}_{\nu}[e^{-1}\partial_{\xi}(ee^{\rho}_{i}S^{\;\;\sigma\xi}_{\rho})-e^{\lambda}_{i}S^{\rho\xi\sigma}T_{\rho\xi\lambda}]$
(14)
$\displaystyle\qquad=G_{\mu\nu}-\frac{1}{2}g_{\mu\nu}T=-\nabla^{\sigma}S_{\nu\sigma\mu}-S_{\;\;\;\;\mu}^{\rho\lambda}K_{\lambda\rho\nu}\;,$

$G_{\mu\nu}$ is the Einstein tensor, and $\nabla_{\sigma}$ is the covariant
derivative associated with the Levi-Civita connection. The trace of Eq. (12)
or (13), which can be used to simplify and constrain the field equation, can
be expressed as

$\displaystyle-[2e^{-1}\partial_{\sigma}(eT^{\;\;\rho\sigma}_{\rho})+T]f_{T}(T)+S^{\;\;\rho\sigma}_{\rho}(\partial_{\sigma}T)f_{TT}(T)+2f(T)=\kappa^{2}\overset{em}{T}\;,$
(15)

where,
$\overset{em}{T}=\overset{em}{T}{}^{\mu}_{\;\;\mu}=g^{\mu\nu}\overset{em}{T}{}_{\mu\nu}$
is the trace of the energy-momentum tensor. Clearly, in the case of TG,
$f(T)=T$, and Eq. (15) reduces to

$T-2e^{-1}\partial_{\sigma}(eT_{\rho}^{\;\;\rho\sigma})=\kappa^{2}\overset{em}{T}\;,$
(16)

which shows an equivalence between GR and TG since

$-R=T-2e^{-1}\partial_{\sigma}(eT_{\rho}^{\;\;\rho\sigma})\;.$ (17)

##  III  vierbein for cylindrical symmetry metric

Since the Gödel-type metric is usually expressed in cylindrical coordinates
$[(r,\phi,z)]$, we consider a general cylindrical symmetry metric

$\displaystyle ds^{2}=dt^{2}+2H(r)dtd\phi-dr^{2}-G(r)d\phi^{2}-dz^{2}\;,$ (18)

where $H$ and $G$ are the arbitrary functions of $r$. This metric can be re-
expressed in the following form

$\displaystyle ds^{2}=[dt+H(r)d\phi]^{2}-D^{2}(r)d\phi^{2}-dr^{2}-dz^{2}\;,$
(19)

where

$\displaystyle D(r)=\sqrt{G(r)+H^{2}(r)}\;.$ (20)

Since the local Lorentz invariance is violated in $f(T)$ theory and the
vierbein have six degrees of freedom more than the metric, one should be
careful in choosing a physically reasonable tetrad in terms of Eq.(2). Here,
we choose the tetrad anstaz of the cylindrical symmetry metric to be:

$\displaystyle e^{i}_{\mu}\equiv\left(\begin{array}[]{cccc}1&0&H&0\\\
0&1&0&0\\\ 0&0&D&0\\\
0&0&0&1\end{array}\right)\;,\;\;\;e^{\mu}_{i}\equiv\left(\begin{array}[]{cccc}1&0&-\frac{H}{D}&0\\\
0&1&0&0\\\ 0&0&\frac{1}{D}&0\\\ 0&0&0&1\\\ \end{array}\right)\;.$ (29)

Using Eqs. (3–9), one can find that the Weitzenböck invariant $T$ is

$\displaystyle T=\frac{1}{2}\left(\frac{H^{\prime}}{D}\right)^{2}\;,$ (30)

where a prime presents a derivative with respect to $r$.

Substituting the vierbein given in Eq. (29) into Eq. (13), we obtain the
following non-zero components of the $f(T)$ field equation:
$\nu=0,i=0$

$\displaystyle\bigg{(}T-\frac{D^{\prime\prime}}{D}+\frac{HT^{\prime}}{2H^{\prime}}\bigg{)}f_{T}(T)+\bigg{(}\frac{HT}{H^{\prime}}-\frac{D^{\prime}}{D}\bigg{)}T^{\prime}f_{TT}(T)+\frac{1}{2}f(T)=\kappa^{2}\overset{em}{T}{}^{0}_{\;\;0}$
(31)

$\nu=0,i=2$

$\displaystyle\bigg{(}HT+\frac{T^{\prime}D^{2}}{2H^{\prime}}\bigg{)}f_{T}(T)+\frac{T^{\prime}H^{\prime}}{2}f_{TT}(T)-\frac{H}{2}f(T)=\kappa^{2}\bigg{(}\overset{em}{T}{}^{0}_{\;\;2}-H\overset{em}{T}{}^{0}_{\;\;0}\bigg{)}\;,$
(32)

$\nu=1,i=1$

$\displaystyle-
Tf_{T}(T)+\frac{1}{2}f(T)=\kappa^{2}\overset{em}{T}{}^{1}_{\;\;1}\;,$ (33)

$\nu=2,i=0$

$\displaystyle
T^{\prime}\bigg{[}\frac{1}{2H^{\prime}}f_{T}(T)+\sqrt{\frac{T}{2}}f_{TT}(T)\bigg{]}=\kappa^{2}\overset{em}{T}{}^{2}_{\;\;0}\;,$
(34)

$\nu=2,i=2$

$\displaystyle-
Tf_{T}(T)+\frac{1}{2}f(T)=\kappa^{2}\bigg{(}\overset{em}{T}{}^{2}_{\;\;2}-H\overset{em}{T}{}^{2}_{\;\;0}\bigg{)}\;,$
(35)

$\nu=3,i=3$

$\displaystyle-\frac{D^{\prime\prime}}{D}f_{T}(T)-\frac{T^{\prime}D^{\prime}}{D}T^{\prime}f_{TT}(T)+\frac{1}{2}f(T)=\kappa^{2}\overset{em}{T}{}^{3}_{\;\;3}\;.$
(36)

Apparently, the non-symmetric components of the modified Einstein equation are
consistent with the tetrad anstaz given in Eq. (29). In the above equations,
all other components of $\overset{em}{T}{}^{\mu}_{\;\;\nu}$ must be zero,
which means that, $\overset{em}{T}{}_{\mu\nu}$, has the cylindrical symmetry
as expected. In a Gödel-type spacetime, the energy-momentum tensor in a local
basis, $\overset{em}{T}_{ab}$ given in (50), has a general form:
$\overset{em}{T}{}_{ab}=diag(\rho,p_{1},p_{2},p_{3})$. Using
$\overset{em}{T}{}_{\mu\nu}=e^{a}_{\mu}e^{b}_{\nu}\overset{em}{T}{}_{ab}$, we
have

$\displaystyle\overset{em}{T}{}_{00}=\rho,\;\;\overset{em}{T}{}_{11}=p_{1},\;\;\;\overset{em}{T}{}_{22}=H^{2}\rho+D^{2}p_{2},\;\;\;\overset{em}{T}{}_{33}=p_{3},\;\;\;\overset{em}{T}{}_{02}=\overset{em}{T}{}_{20}=H\rho\;.$
(37)

One can then find easily

$\displaystyle\overset{em}{T}{}^{2}_{\;\;0}=0,\;\;\;\;\overset{em}{T}{}^{0}_{\;\;2}=H\bigg{(}\overset{em}{T}{}^{0}_{\;\;0}-\overset{em}{T}{}^{2}_{\;\;2}\bigg{)}\;\;.$
(38)

Thus, Eq. (32) seems to give an extra constraint on $f(T)$ gravity. This
equation is satisfied automatically in a Gödel-type spacetime, since $T$, as
shown in Eq. (42), is a constant in a Gödel-type universe. Furthermore, it is
easy to see that, in a Gödel-type spacetime, Eq. (32) gives the same
expression as Eq. (35). Four independent field equations are obtained, which
is consistent with the anstaz of tetrad. In addition, one can check that the
field equations (23-28) for the vierbein given in (29) can also be obtained
from an action constructed by replacing the specific form of $T$ (30) with the
general action of $f(T)$ theory. Therefore, the dynamical equations are
consistent, which means that the tetrad anstaz given in Eq. (29) is a good
guess for the Gödel-type spacetime.

##  IV Gödel-type universe in $f(T)$ theory

To show the possibility of existence of the closed timelike curves and the
causality feature in $f(T)$ gravity, we consider the Gödel-type metric, which
has the form of Eq. (18) with $H$ and $G$ being:

$\displaystyle
H(r)=\frac{4\omega}{m^{2}}\sinh^{2}\bigg{(}\frac{mr}{2}\bigg{)}\;,$ (39)
$\displaystyle
G(r)=\frac{4}{m^{2}}\sinh^{4}\bigg{(}\frac{mr}{2}\bigg{)}\bigg{[}\coth^{2}\bigg{(}\frac{mr}{2}\bigg{)}-\frac{4\omega^{2}}{m^{2}}\bigg{]}\;,$
(40)

where $\omega$ and $m$ ($-\infty<m^{2}<+\infty,0<\omega^{2}$) are two constant
parameters used to classify different Gödel-type geometries. Thus, we have

$\displaystyle D(r)=\frac{1}{m}\sinh(mr)\;.$ (41)

Substituting the expressions of $H$ and $D$ into Eq. (30), one can obtain
easily

$\displaystyle T=2\omega^{2}\;,$ (42)

which is a positive constant.

If $G(r)<0$, Eq. (18) shows that one type of closed timelike curve, called
noncausal Gödel circle Godel , exists in the case of $t,z,r=const$. This means
a violation of causality. For a particular case of $0<m^{2}<4\omega^{2}$, the
causality violation region, i.e., $G(r)<0$ region, exists if

$\displaystyle\tanh^{2}\frac{mr}{2}<\frac{m^{2}}{4\omega^{2}}\;.$ (43)

Thus, one can define a critical radius $r_{c}$ Godel ; RebClif05e ; Reb09b1 ;
Reb10ed

$\displaystyle\tanh^{2}\frac{mr_{c}}{2}=\frac{m^{2}}{4\omega^{2}}\;,$ (44)

beyond which, $G(r)<0$ and causality is violated. When $m=0$, the critical
radius is $r_{c}=1/\omega$. When $m^{2}=4\omega^{2}$, $r_{c}=+\infty$, which
means that a breakdown of causality is avoided. Thus, the codomain range of
$r_{c}$ is $r_{c}\in(1/\omega,+\infty)$. Therefore, the condition for non-
violation of causality is $m^{2}\geq 4\omega^{2}$ or $r<r_{c}$. For the case
in which $m^{2}=-\mu^{2}<0$, both
$H(r)=\frac{4\omega}{\mu^{2}}\sin^{2}(\frac{\mu r}{2})$ and
$G(r)=\frac{4}{\mu^{2}}\sin^{4}(\frac{\mu r}{2})[\cot^{2}(\frac{\mu
r}{2})-\frac{4\omega^{2}}{\mu^{2}}]$ are periodic functions. Thus, an infinite
circulation of causal and noncausal ranges appears Reb09b1 ; Reb10ed .

It is easy to see that, if one further defines a set of bases
$\\{\theta^{a}\\}$:

$\displaystyle\theta^{0}=dt+H(r)d\phi,\qquad\theta^{1}=dr,$ (45)
$\displaystyle\theta^{2}=D(r)d\phi,\qquad\theta^{3}=dz,$ (46)

the Goödel-type line element can be simplified to be:

$\displaystyle ds^{2}=\eta_{ab}\theta^{a}\theta^{b}\;,$ (47)

where $\eta_{ab}=diag(1,-1,-1,-1)$ is the Minkowski metric. By choosing
$\\{\theta^{a}\\}$ as basis, the $f(T)$ field equation (13) becomes:

$\displaystyle
A_{ab}f_{T}(T)+\frac{1}{2}\eta_{ab}f(T)=\kappa^{2}\overset{em}{T}_{ab}\;.$
(48)

Here, both $f(T)$ and $f_{T}(T)$ are evaluated at $T=2\omega^{2}$. The second
term of Eq. (13) is discarded in obtaining the above equation since the
torsion scalar $T$ is a constant. We find that the nonzero components of
$A_{ab}$ are

$\displaystyle A_{00}=2\omega^{2}-m^{2},\quad A_{11}=A_{22}=2\omega^{2},\quad
A_{33}=m^{2}\;.$ (49)

Thus, we obtain a very simple form of the field equation in $f(T)$ gravity,
which will help us discuss the causality issue.

##  V Causality Problem in $f(T)$ theory

One can see, from Eq. (48), that, in order to discuss the causality problem,
the matter source is a very important component. As was obtained in RebClif05e
; Reb09b1 ; Reb10ed , different matter sources may lead to different results.
In this paper, we assume that the matter source consists of two different
components: a perfect fluid and a scalar field. Thus, the energy-momentum
tensor $\overset{em}{T}{}_{ab}$ has the form

$\displaystyle\overset{em}{T}_{ab}=\overset{m}{T}_{ab}+\overset{s}{T}_{ab}\;,$
(50)

where, $\overset{m}{T}_{ab}$ and $\overset{s}{T}_{ab}$ correspond to the
energy-momentum tensors of the perfect-fluid and the scalar field,
respectively. In basis $\\{\theta^{a}\\}$, $\overset{m}{T}_{ab}$ and
$\overset{s}{T}_{ab}$ can be expressed as

$\displaystyle\overset{m}{T}_{ab}=(\rho+p)u_{a}u_{b}-p\eta_{ab}\;,$ (51)
$\displaystyle\overset{s}{T}_{ab}=D_{a}\Phi
D_{b}\Phi-\frac{1}{2}\eta_{ab}D_{c}\Phi D_{d}\Phi\eta^{cd}\;,$ (52)

where $u_{a}=(1,0,0,0)$, $\rho$ and $p$ are the energy density and pressure of
the perfect fluid, respectively, and $p=\text{w}\rho$ with w being the
equation of state parameter. $\Phi$ is the scalar field, and $D_{a}$ denotes
the covariant derivative relative to the local basis $\theta^{a}$. The scalar
field equation is $\square\,\Phi=\eta^{ab}\,\nabla_{a}\nabla_{b}\,\Phi\,=0$.
It is easy to prove that $\Phi(z)=\varepsilon z+\text{const}$ with a constant
amplitude $\varepsilon$ satisfies this field equation Reboucas1983 . Using the
solution $\Phi(z)=\varepsilon z+\text{const}$, one can obtain the nonvanishing
components of $\overset{s}{T}_{ab}$

$\overset{s}{T}_{00}=-\overset{s}{T}_{11}=-\overset{s}{T}_{22}=\overset{s}{T}_{33}=\frac{\varepsilon^{2}}{2}\,,$
(53)

Thus, the energy-momentum tensor of matter source becomes

$\displaystyle\overset{em}{T}_{ab}=diag\bigg{(}\rho+\frac{\varepsilon^{2}}{2}\;,\text{w}\rho-\frac{\varepsilon^{2}}{2}\;,\text{w}\rho-\frac{\varepsilon^{2}}{2}\;,\text{w}\rho+\frac{\varepsilon^{2}}{2}\bigg{)}\;.$
(54)

Substituting Eqs. (49) and (54) into the $f(T)$ field equation (Eq. (48)), we
find

$\displaystyle(2\omega^{2}-m^{2})f_{T}(T)+\frac{1}{2}f(T)=\kappa^{2}(\rho+\frac{\varepsilon^{2}}{2})\;;$
(55) $\displaystyle
2\omega^{2}f_{T}(T)-\frac{1}{2}f(T)=\kappa^{2}(\text{w}\rho-\frac{\varepsilon^{2}}{2})\;;$
(56) $\displaystyle
m^{2}f_{T}(T)-\frac{1}{2}f(T)=\kappa^{2}(\text{w}\rho+\frac{\varepsilon^{2}}{2})\;.$
(57)

Since the effective Newton gravity constant in $f(T)$ gravity becomes
$G_{N,eff}=G_{N}/f_{T}(T)$ Zheng2011 , only the case $f_{T}(T)>0$ will be
considered in the following in order to ensure a positive $G_{N,eff}$. From
Eqs. (55) and (56), one can derive a relation between $m$ and $\omega$:

$\displaystyle
m^{2}=2\omega^{2}\bigg{[}1+\frac{\varepsilon^{2}}{\rho(1+\text{w})+\varepsilon^{2}}\bigg{]}\;,$
(58)

which implies that the critical radius of the Gödel’s circle, Eq. (44),
satisfies

$\displaystyle\tanh^{2}\left(\frac{mr_{c}}{2}\right)=1-\frac{\rho(1+\text{w})}{2[\rho(1+\text{w})+\varepsilon^{2}]}\;.$
(59)

Obviously, different matter sources give rise to different critical radii and
therefore different causality structures, e.g. when $\varepsilon\rightarrow
0$, we have a finite $r_{c}$, while for $\rho\rightarrow 0$, $r_{c}=\infty$.
Therefore, a violation of causality may occur for the case of a perfect fluid
as the matter source, whereas causality is preserved in the case of a scalar
field. In order to show the causality feature in more detail and the
conditions for obtaining the Gödel-type solutions, we will divide our
discussion into two special cases: $\varepsilon\rightarrow 0$ and
$\rho\rightarrow 0$. In addition, a concrete $f(T)$ model will be considered.

###  V.1  $\varepsilon^{2}\rightarrow 0$

$\varepsilon^{2}\rightarrow 0$ corresponds to the case that the universe only
contains a perfect fluid. Since $f_{T}(T)>0$, Eqs. (55), (56), and (57) reduce
to:

$\displaystyle m^{2}=2\omega^{2}\;;$ (60) $\displaystyle
Tf_{T}(T)=\kappa^{2}\rho(1+\text{w})\;;$ (61) $\displaystyle
f(T)=2\kappa^{2}\rho\;.$ (62)

From Eqs. (61, 62), it is easy to see that, in the limit of general relativity
without a cosmological constant ($f(T)=T$), $\text{w}=1$ is required to ensure
the existence of the Gödel-type solutions godelnote1978 ; RebClif05e ; Reb09b1
; Reb10ed . This means that a violation of causality in general relativity is
only possible for the so-called stiff fluid ($\text{w}=1$) which is not a
normal fluid in our Universe. In $f(T)$ theory, $Tf_{T}(T)>0$ and $\rho>0$
lead to $\text{w}>-1$. So, the perfect fluid must satisfy the weak energy
condition ($\rho>0$ and $\rho(1+\text{w})>0$). Using the above results, the
equation of state can be expressed as a function of the torsion scalar:

$\displaystyle\text{w}=\frac{2Tf_{T}(T)}{f(T)}-1\;.$ (63)

Different from general relativity that requires $\text{w}=1$ for perfect-fluid
Gödel solutions, the equation of state parameter of the fluid w in $f(T)$
gravity can differ from one and its value is determined by concrete $f(T)$
models. For example, a special $f(T)=\lambda T^{\delta}$ gives
$\text{w}=2\delta-1$, from which one can see that w can be an arbitrary number
for an arbitrary $\delta$. So, even normal matter, such as pressureless matter
or radiation, can lead to a violation of causality in certain $f(T)$ theories.
This indicates that the issue of causality violation seems more severe in
$f(T)$ gravity than in general relativity where only an exotic stiff fluid
allows the existence of Gödel-type solutions. From Eqs. (60), (61) and (62),
and using $T=2\omega^{2}$, we find that the critical radius given in Eq. (59)
becomes

$\displaystyle
r_{c}=2\text{tanh}^{-1}\bigg{(}\frac{1}{\sqrt{2}}\bigg{)}\cdot\sqrt{\frac{f_{T}(T)}{(1+\text{w})\kappa^{2}\rho}}\;,$
(64)

which is dependent both on the specifics of $f(T)$ theory and the properties
of the perfect fluid.

Now, let us consider a concrete power law $f(T)$ model Linder

$\displaystyle f(T)=T-\alpha T_{*}\left(\frac{T}{T_{*}}\right)^{n}\;,$ (65)

where $\alpha$ and $n$ are model parameters, and $T_{*}$ is a special value of
the torsion scalar, which is introduced to make $\alpha$ dimensionless.
$|n|\ll 1$ is required in order to obtain an usual early cosmic evolution
Wu2010b . The current cosmic observations give that
$\alpha=-0.79^{+0.35}_{-0.79}$ and $n=0.04^{+0.22}_{-0.33}$ at the $68.3\%$
confidence level Wu2010a . Thus, a negative $\alpha$ is favored by
observations. In term of Eq. (63), the equation of state of the perfect fluid
becomes

$\displaystyle\text{w}=1-\frac{2\alpha(n-1)T^{1-n}_{*}}{T^{1-n}-\alpha
T^{1-n}}\;.$ (66)

The equation above can be re-expressed as

$\displaystyle\frac{\alpha(2n-1-\text{w})}{1-\text{w}}=\left(\frac{T}{T_{*}}\right)^{1-n}>0\;,$
(67)

where a positive $T/T_{*}$ is considered. Recalling $\alpha<0$ and
$\text{w}>-1$, from Eq. (67) one can obtain the possible ranges of w for the
Gödel-type universes

$\displaystyle 1>\text{w}>-1+2n\quad(1>n>0)\;,\qquad
1>\text{w}>-1\quad(n<0)\;.$ (68)

For this power law model, the critical radius has the form

$\displaystyle
r_{c}=2\left[\frac{\alpha(2n-1-\text{w})}{1-\text{w}}\right]^{\frac{1}{2(n-1)}}\text{tanh}^{-1}(1/\sqrt{2})\;,$
(69)

which is determined completely by the model parameters and the equation of
state of the perfect fluid.

###  V.2  $\rho\rightarrow 0$

This is the case of a scalar field as the matter source. Eqs. (55), (56), and
(57) now reduce to

$\displaystyle m^{2}=4\omega^{2}\;,$ (70) $\displaystyle
Tf_{T}(T)=\kappa^{2}\varepsilon^{2}\;,$ (71) $\displaystyle
f(T)=3\kappa^{2}\varepsilon^{2}\;.$ (72)

Note that (71) and (72) combined together admit a relation between $T$ and
$f(T)$:

$\displaystyle 3Tf_{T}(T)-f(T)=0\;,$ (73)

which constrains the class of solutions with no violation of causality. For
the power law model, the causal Gödel-type solution gives that the torsion
scalar should satisfy

$\displaystyle
T=2\omega^{2}=\left[-\frac{(1-3n)\alpha}{2}\right]^{\frac{1}{1-n}}T_{*}\;.$
(74)

Thus, $n<1/3$ is required if the numerator of $\frac{1}{1-n}$ is not even
since the observations show $\alpha<0$.

##  VI Conclusions

$f(T)$ theory, a new modified gravity, provides an alternative way to explain
the present accelerated cosmic acceleration with no need of dark energy. Some
problems, including large scale structure, local Lorentz invariance, and so
on, of this modified gravity have been discussed. In this paper, we study the
issue of causality in $f(T)$ theory by examining the possibility of existence
of the closed timelike curves in the Gödel metric. Assuming that the matter
source is a scalar field or a perfect fluid, we examine the existence of the
Gödel-type solutions. For the scalar field case, we find that $f(T)$ gravity
allows a particular Gödel-type solution with $r_{c}\rightarrow\infty$, where
$r_{c}$ is the critical radius beyond which the causality is broken down.
Thus, the violation of causality can be forbidden. In the case of a perfect
fluid as the matter source, we find that the fluid must have an equation of
state parameter greater than minus one and this parameter should satisfy Eq.
(63) for the Gödel-type solutions to exist. For certain $f(T)$ models, the
perfect fluid that allows the Gödel-type solutions can even be normal matter,
such as pressureless matter or radiation. Since the critical radius $r_{c}$ of
perfect fluid Gödel-type solutions which depends on both matter and gravity is
finite, the issue of causality violation seems more severe in $f(T)$ gravity
than in general relativity where only an exotic stiff fluid allows the
existence of Gödel-type solutions.

###### Acknowledgements.

PXW would like to thank Prof. Qingguo Huang for helpful discussions. This work
was supported by the National Natural Science Foundation of China under Grants
Nos. 10935013, 11175093 and 11075083, Zhejiang Provincial Natural Science
Foundation of China under Grants Nos. Z6100077 and R6110518, the FANEDD under
Grant No. 200922, the National Basic Research Program of China under Grant No.
2010CB832803, the NCET under Grant No. 09-0144, and K.C. Wong Magna Fund in
Ningbo University.

## References

  * (1) A. Einstein, Sitzungsber. Preuss. Akad. Wiss. Phys. Math. Kl., 217 (1928).
  * (2) A. Einstein, Math. Ann. 102, 685 (1930).
  * (3) C. Pellegrini and J. Pl.ebanski, K. Dan. Vidensk. Selsk.Math. Fys. Skr. 2, 4 (1962).
  * (4) C. Møller, K. Dan. Vidensk. Selsk. Math. Fys. Skr. 89, 13 (1978).
  * (5) K. Hayashi and T. Nakano. Prog. Theor. Phys. 38, 491(1967).
  * (6) K. Hayashi and T. Shirafuji, Phys. Rev. D 19, 3524 (1979).
  * (7) K. Hayashi and T. Shirafuji, Phys. Rev. D 24, 3312 (1981).
  * (8) G. R. Bengochea and R. Ferraro, Phys. Rev. D 79, 124019 (2009).
  * (9) R. Ferraro and F. Fiorini, Phys. Rev. D 75, 084031 (2007).
  * (10) R. Ferraro and F. Fiorini, Phys. Rev. D 78, 124019 (2008).
  * (11) R. Ferraro and F. Fiorini, arXiv:1106.6349.
  * (12) P. Wu and H. Yu, Phys. Lett. B 703, 223 (2011).
  * (13) P. Wu and H. Yu, Eur. Phys. J. C 71, 1552 (2011).
  * (14) K. Bamba, C. Geng, C. Lee and L. Luo, JCAP 1101, 021 (2011).
  * (15) P. Wu and H. Yu, Phys. Lett. B 693, 415 (2010).
  * (16) G. R. Bengochea, Phys. Lett. B 695, 405 (2011).
  * (17) P. Wu and H. Yu, Phys. Lett. B 692, 176 (2010).
  * (18) Y. Zhang, H. Li, Y. Gong and Z. Zhu, JCAP 07, 015(2011).
  * (19) R. Ferraro, F. Fiorini, Phys. Lett. B 702, 75 (2011).
  * (20) E. V. Linder, Phys. Rev. D 81, 127301 (2010).
  * (21) K. K. Yerzhanov, S. R. Myrzakul, I. I. Kulnazarov, R. Myrzakulov, arXiv:1006.3879.
  * (22) R. Yang, Eur. Phys. J. C 71, 1797 (2011).
  * (23) P. Y. Tsyba, I. I. Kulnazarov, K. K. Yerzhanov, R. Myrzakulov, Int. J. Theor. Phys. 50, 1876 (2011).
  * (24) S. Chen, J. B. Dent, S. Dutta, E. N. Saridakis, Phys.Rev.D 83, 023508 (2011).
  * (25) K. Bamba, C. Geng, C. Lee, arXiv:1008.4036.
  * (26) R. Myrzakulov, arXiv:1008.4486.
  * (27) K. Karami, A. Abdolmaleki, arXiv:1009.2459.
  * (28) K. Karami, A. Abdolmaleki, arXiv:1009.3587.
  * (29) R. Yang, Europhys. Lett. 93, 60001 (2011).
  * (30) J. B. Dent, S. Dutta, E. N. Saridakis, JCAP 1101, 009 (2011).
  * (31) T. Wang, Phys. Rev. D 84, 024042 (2011).
  * (32) Y. Cai, S. Chen, J. B. Dent, S. Dutta, E. N. Saridakis, Class. Quantum Grav. 28, 215011 (2011).
  * (33) S. Chattopadhyay, U. Debnath, Int. J. Mod. Phys. D 20, 1135 (2011).
  * (34) H. Wei, X. Ma, H. Qi, Phys. Lett. B 703, 74 (2011).
  * (35) X. Meng, Y. Wang, Eur. Phys. J. C 71 1755 (2011).
  * (36) C. G. Boehmer, A. Mussa, N. Tamanini, Class. Quant. Grav. 28 245020 (2011).
  * (37) H. Wei, H. Qi, X. Ma, arXiv:1108.0859.
  * (38) S. Capozziello, V. F. Cardone, H. Farajollahi, A. Ravanpak, Phys. Rev. D 84, 043527 (2011).
  * (39) C. Geng, C. Lee, E. N. Saridakis, and Y. Wu, Phys. Lett. B 704, 384 (2011).
  * (40) K. Bamba, C. Geng, JCAP 1111, 008 (2011).
  * (41) R. Ferraro, F. Fiorini, Phys. Rev. D 84, 083518 (2011).
  * (42) H. Wei, Phys. Lett. B 712 430 (2012).
  * (43) Y. Wu, C. Geng, arXiv:1110.3099.
  * (44) P. A. Gonzalez, E. N. Saridakis, and Y. Vasquez, JHEP 1207 053 (2012).
  * (45) C. G. Boehmer, T. Harko, and F. S. N. Lobo, Phys. Rev. D 85, 044033 (2012).
  * (46) K. Karami, A. Abdolmaleki, arXiv:1111.7269.
  * (47) H. Wei, X. Guo, L. Wang, Phys. Lett. B 707, 298 (2012).
  * (48) K. Atazadeh, F. Darabi, Eur. Phys. J. C 72 2016 (2012).
  * (49) H. Farajollahi, A. Ravanpak, P. Wu, Astrophys. Space Sci. 338 23 (2012).
  * (50) K. Karami, A. Abdolmaleki, JCAP 1204 007 (2012).
  * (51) J. Yang, Y. Li, Y. Zhong, Y. Li, Phys. Rev. D 85 084033 (2012).
  * (52) W. Zhang, et al., arXiv:1202.0892.
  * (53) Z. Haghani, T. Harko, H. Sepangi, S. Shahidi, arXiv:1202.1879.
  * (54) K. Karami, A. Abdolmaleki, arXiv:1202.2278.
  * (55) C. Xu, E. N. Saridakis, G. Leon, JCAP 1207 005 (2012).
  * (56) M. R. Setare, M. J. S. Houndjo arXiv:1203.1315.
  * (57) M. Hamani Daouda, Manuel E. Rodrigues and M. J. S. Houndjo Eur. Phys. J. C 71: 1817 (2011).
  * (58) M. Hamani Daouda, Manuel E. Rodrigues and M. J. S. Houndjo Eur. Phys. J. C 72 1890 (2012).
  * (59) M. Hamani Daouda, Manuel E. Rodrigues and M. J. S. Houndjo Eur. Phys. J. C 72 1893 (2012).
  * (60) M. Hamani Daouda, Manuel E. Rodrigues and M. J. S. Houndjo arXiv:1202.1147.
  * (61) B. Li, T. P. Sotiriou and J. D. Barrow, Phys. Rev. D 83, 064035 (2011).
  * (62) T. P. Sotiriou, B. Li, J. D. Barrow, Phys. Rev. D 83, 104030 (2011).
  * (63) M. Li, R. Miao, Y. Miao, JHEP 07, 108 (2011).
  * (64) R. Miao, M. Li, Y. Miao, JCAP 11, 033 (2011).
  * (65) B. Li, T. P. Sotiriou, J. D. Barrow, Phys. Rev. D 83, 104017 (2011).
  * (66) R. Zheng and Q. Huang, JCAP 03, 002 (2011).
  * (67) A. G. Riess, et al., Astron. J. 116, 1009 (1998).
  * (68) S. Perlmutter, et al., Astrophys. J. 517, 565 (1999).
  * (69) D. N. Spergel, et al., ApJS, 148, 175 (2003).
  * (70) D. N. Spergel, et al., ApJS, 170, 377 (2007).
  * (71) M. Tegmark, et al., Phys. Rev. D 69, 103501 (2004).
  * (72) D. J. Eisenstein, et al., Astrophys. J. 633, 560 (2005).
  * (73) S. Nojiri, S. D. Odintsov, arXiv:0807.0685.
  * (74) T. P. Sotiriou, V. Faraoni, Rev. Mod. Phys. 82, 451 (2010).
  * (75) A. De Felice and S. Tsujikawa, Living Rev. Rel. 13, 3 (2010).
  * (76) S. Nojiri, S. D. Odintsov, Phys. Rept. 505, 59 (2011).
  * (77) T. Clifton, P. G. Ferreira, A. Padilla, C. Skordis, Phys. Rept. 513, 1 (2012).
  * (78) K. Gödel, Rev. Mod. Phys. 21, 447 (1949).
  * (79) M. M Som and A. K. Raychaudhuri, Proc. R. Soc. London A 304, 81 (1968).
  * (80) M. J. Reboucas, Phys. Lett. A 70, 161 (1979).
  * (81) A. K. Raychaudhuri and S. N. Guha Thakurta, Phys. Rev. D 22, 802 (1980).
  * (82) A. F. F. Teixeira, M. J. Rebouças and J. E. Aman, Phys. Rev. D 32 3309 (1985).
  * (83) M. J. Rebouças and A. F. F. Teixeira, J. Math. Phys. 33, 2855 (1992).
  * (84) H. L. Carrion, M. J. Rebouças, and A. F. F. Teixeira, J. Math. Phys. 40, 4011 (1999).
  * (85) W. A. Hiscock, Phys. Rev. D 17, 1497 (1978).
  * (86) S. K. Chakraborty and N. Bandyopadhyay, J. Math. Phys. 24, 129 (1983).
  * (87) V. F. Panov, Sov. Phys. J. 34, 130 (1991).
  * (88) F. Bampi and C. Zordan, Gen. Rel. Grav. 9, 393 (1978).
  * (89) V. M. Villalba, Mod. Phys. Lett. A 8, 3011 (1993).
  * (90) L. O. Pimentel, A. Camacho and A. Macias, Mod. Phys. Lett. A 9, 3703 (1994).
  * (91) V. G. Krechet and I. V. Sandina, Sov. Phys. J. 24, 950 (1982).
  * (92) D. A. Leahy, Int. J. Theor. Phys. 21, 703 (1982).
  * (93) A. Herrera and G. N. Shikin, Grav. Cosmol. 2, 155 (1996).
  * (94) M. J. Rebouças and J. Tiomno, Phys. Rev. D 28, 1251(1983).
  * (95) M. J. Rebouças, J. E. Aman and A. F. F. Teixeira, J. Math. Phys. 27, 1370 (1986).
  * (96) Y. N. Obukhov, T. Vargas Phys. Lett. A 327, 365 (2004).
  * (97) T. Clifton and J. D. Barrow, Phys. Rev. D 72, 123003 (2005).
  * (98) M. J. Rebouças, J. Santos, Phys. Rev. D 80, 063009 (2009).
  * (99) J. Santos, M. J. Rebouças, T. B. R. F. Oliveira, Phys. Rev. D 81, 123017 (2010).
  * (100) P. Kanti and C. E. Vayonakis, Phys. Rev. D 60, 103519 (1999) .
  * (101) J. D. Barrow and M. P. Da̧browski Phys. Rev. D 58 103502 (1998)


# Stopping a reaction-diffusion front

Jean-Guy Caputo  Laboratoire de Mathématiques, INSA de Rouen - B.P. 8, Avenue
de l’Université, 76801 Saint-Etienne du Rouvray, France. E-mail: caputo@insa-
rouen.fr  Gustavo Cruz-Pacheco  Depto. Matemáticas y Mecánica,
I.I.M.A.S.-U.N.A.M., Apdo. Postal 20–726, 01000 México D.F., México  Benoît
Sarels  Sorbonne Université, CNRS, Université de Paris, Laboratoire Jacques-
Louis Lions (LJLL), F-75005 Paris, France

###### Abstract

We revisit the problem of pinning a reaction-diffusion front by a defect, in
particular by a reaction-free region. Using collective variables for the front
and numerical simulations, we compare the behaviors of a bistable and
monostable front. A bistable front can be pinned as confirmed by a pinning
criterion, the analysis of the time independant problem and simulations.
Conversely, a monostable front can never be pinned, it gives rise to a
secondary pulse past the defect and we calculate the time this pulse takes to
appear. These radically different behaviors of bistable and monostable fronts
raise issues for modelers in particular areas of biology, as for example, the
study of tumor growth in the presence of different tissues.

##  I Introduction

Reaction-diffusion equations are general models for active chemical reactions
like combustion zeldovich , see the book by Alwyn Scott scott  for examples,
or the review xin00 . In biology, while these equations are ubiquitous in
population dynamics or population genetics, they are also used in a variety of
problems such as nerve impulse propagation in an axon or growth of a cancer
tumor cancer , see the general book by Murray murray . These equations have a
number of stationary solutions. In the physics context, the models have three
stationary solutions, two stable and one unstable, they are termed bistable.
In biology, the bistable model can be used to describe the growth of a
population. Still, the main model remains the one that describes the wave of
advance of an advantageous gene like in the seminal paper fisher . It only has
two stationary states, one stable and another unstable; it is called
monostable. The normal forms of these two different reaction terms are cubic
for the bistable and quadratic for the monostable. Mostly, the study of these
two nonlinearities gives the qualitative picture for almost all other
nonlinearities.

Important dynamical solutions are the fronts (in 1D) that connect two such
stationary solutions. In 2D and 3D, assuming radial symmetry, fronts describe
the interface of ”blobs” where the inside has one state and the outside
another. Fronts have a speed that is proportional to the square root of the
diffusion times the reaction rate. Many studies dealt with the stability of
such fronts, see scott  for references. Assuming a front is stable, important
parameters are its position and width i.e. the spatial extension of the
transition region separating the two stationary states. These fronts have
exact forms for the bistable scott  and monostable ablowitz  nonlinearities.
Using these solutions as Rayleigh-Ritz type ansatzes, we and other authors
have derived ordinary differential equations giving the evolution of the
position and width -termed collective coordinates- whose solutions match
remarkably well the full dynamics, see cs1 ; susanto  for the bistable case
and bcp14  for the monostable case.

For many applications, it is useful to understand how fronts respond to
perturbations of the environment. Such perturbations can be temporal such as
the action of therapy on a cancer cell cancer2 . They can also be spatial like
a defect in a burning candle. Such a geographic defect can act on scales
larger or smaller than the front width. In the latter case, the front
adiabatically adapts to the defect and changes its width and speed
accordingly, see for example our study cs1  on bistable fronts. When the
defect is narrow compared to the front, a bistable front can be pinned by the
defect cs1 . This effect can also be seen in waveguides whose transverse width
abruptly changes, see the works bbc16  and thesebouhours .

In this article, we revisit the issue of front pinning comparing bistable and
monostable fronts. As expected, wide defects cannot pin any of the two types
of fronts. Narrow defects will pin a bistable front, see the numerics and
analysis of cs1 . To understand further the phenomenon, we consider as defect
a no reaction region of a given extension. Using collective variables for the
front, we find a pinning criterion which indicates that static bistable fronts
exist in such a region if it is large enough. This approach is confirmed by
analyzing the time independant problem and the solution of the 1D partial
differential equation. Monostable fronts behave very differently, they can
never be pinned. Numerical simulations of monostable fronts appear to stop the
front but the wave continues to propagate and gives rise to a secondary pulse,
past the defect. We calculated the time of delay, i.e. the difference between
the arrival at the defect and the appearance of the secondary pulse, when its
maximum reaches 0.5. This delay time scales approximately linearly with the
extension of the zero reaction zone, in agreement with a simple calculation
based on the diffusion kernel and the instability rate of the zero stationary
solution.
The article is organized as follows. The bistable and monostable models and
their exact solutions are recalled in section 2. Section 3 presents and
discusses the collective variable differential equations for the front
position and width. Sections 4 and 5 detail the front motion in a no reaction
zone for the bistable and monostable models respectively. Section 6 concludes
the article.

##  II The model

In the following we are concerned with the equation

$\left\\{\begin{aligned} u_{t}&=u_{xx}+s(x)~{}R(u)\\\ u(x,0)&=u_{0}(x)\\\
\end{aligned}\right.$ (1)

for $(x,t)\in\mathbb{R}\times\mathbb{R}_{+}^{*}$ and $u_{0}$ a given function
with sufficient regularity. We will study the two canonical types of non-
linearities: the monostable one $R(u)=u(1-u)$ and the bistable one
$R(u)=u(1-u)(u-a)$ for some $0<a<1$. To account for the variable growth rate
of the quantity $u$ (chemical density, population density…), we make use of a
reactivity $s$ that is space-dependent and remains positive over the
considered range: $s(x)>0$ for every real $x$. Let us first recall the theory
when $s$ is constant. We recall that in the former case, the state $u^{*}=1$
is a stable equilibrium of the associated ordinary differential equation while
$u^{*}=0$ is an unstable equilibrium. In the latter, the states $u^{*}=0,1$
are stable equilibria while $u^{*}=a$ is unstable.

To carry out our analysis in the next sections, we will make extensive use of
the fact that exact solutions of equation (1) are known when $s$ is constant.
In the bistable case, all traveling wave solutions are fronts connecting the
two stable equilibria. We choose to consider in the rest of the paper only
fronts going from $1$ at $-\infty$ to $0$ at $+\infty$. Then all traveling
wave solutions - up to a translation - are of the form

$u(x,t)=\tilde{U}_{bi}(x-ct)=\frac{1}{1+\exp{\left(\sqrt{\frac{s}{2}}(x-ct)\right)}},$
(2)

where the speed $c$ is known scott  and related to the parameter $a$ _via_ the
formula

$c=\sqrt{\frac{s}{2}}(1-2a).$ (3)

To consider only positive speeds, we will restrain - without loss of
generality - our analysis to the range $0<a<\frac{1}{2}$. In this case, the
propagation of the front translates as an invasion of the state $0$ by the
state $1$.

At each time, the front is centered around $x=ct$, this means
$u(ct,t)=\tilde{U}_{bi}(0)=\frac{1}{2}$. We can then define the width $w$ of a
given front by the relation

$\tilde{U}_{bi}(-w)=\frac{1}{1+e^{-1}}\approx 0.73$

or equivalently $\tilde{U}_{bi}(w)=\frac{1}{1+e}\approx 0.27$. With that
definition we have here

$w=\sqrt{\frac{2}{s}}$ (4)

In the monostable case, the situation is somewhat different. There are
traveling wave solutions for a continuum of speeds $c\geq 2\sqrt{s}$, and here
again traveling wave solutions are fronts going from $1$ at $-\infty$ to $0$
at $+\infty$. It is known that asymptotically, a great number of solutions
converge to the front of minimal speed $c=2\sqrt{s}$. That is notably the case
when $u_{0}$ is compactly supported roquejoffre1994 . Only one family of exact
solutions is known ablowitz , they take the form

$u(x,t)=\tilde{U}_{mono}(x-ct)=\frac{1}{\left(1+K~{}\exp\left(\sqrt{\frac{s}{6}}(x-ct)\right)\right)^{2}}$
(5)

where $K$ is a constant and the speed is given by

$c=5\sqrt{\frac{s}{6}}.$ (6)

We will use this solution with $K=\sqrt{2}-1$ to ensure that
$\tilde{U}_{mono}(0)=\frac{1}{2}$. Here we define the width $w$ of a given
front by the relation

$\tilde{U}_{mono}(-w)=\frac{1}{(1+(\sqrt{2}-1)e^{-1})^{2}}\approx 0.75$

or equivalently $\tilde{U}_{mono}(w)=\frac{1}{(1+(\sqrt{2}-1)e)^{2}}\approx
0.22$. With that definition we have

$w=\sqrt{\frac{6}{s}}$ (7)

Some differences between the bistable case and the monostable one are better
understood when we see waves in the bisable case as being pushed by the bulk
of the population distribution and waves in the monostable case wave as being
pulled by the leading edge of the distribution garnier2012 .

##  III Collective variables

In the model we are interested in, the reactivity $s$ is allowed to be space-
dependent. When the variations of $s$ are small, it is reasonable to expect
that the solutions of (1) will remain close to the solutions in the
homogeneous case. This approach is sometimes called the use of collective
variables. First, we can expect that the front remains close to its original
profile but will move with a modulation of its speed:

$u(x,t)\approx\tilde{U}(x-c(t)t)$

Actually, the numerical simulations indicate that the solution behavior is
better captured through a modulation of its center $x_{0}$ and width $w$:

$u(x,t)\approx\tilde{U}\left(\frac{x-x_{0}(t)}{w(t)}\right)$

When investigating the evolution of a compactly supported initial condition,
for instance growing from near $0$ values, it can also be useful to allow for
a modulation of the amplitude in the form of the following ansatz:

$u(x,t)\approx A(t)\tilde{U}\left(\frac{x-x_{0}(t)}{w(t)}\right)$

This kind of approach is used for instance in bcp14 , where evolution
equations for the collective variables are obtained through balance laws, just
like in cs1 . As our goal is to investigate the dynamics of established
fronts, and not the evolution from an initial condition to a generalized
traveling front, we will make no use here of this third collective variable.

That being said, in the context of this paper, we will assume that at order
zero we have

$u(x,t)=U(x,x_{0}(t),w(t)):=(\tilde{U}\circ z)(x,x_{0}(t),w(t))$ (8)

where $\tilde{U}$ is a given profile and we define
$z(x,x_{0}(t),w(t))=\frac{x-x_{0}(t)}{w(t)}$. We can compute the time
derivative and the second order space derivative

$\left\\{\begin{aligned} u_{t}&=\dot{x_{0}}\frac{\partial U}{\partial
x_{0}}+\dot{w}\frac{\partial U}{\partial w}\\\
u_{xx}&=\frac{\partial^{2}U}{\partial x^{2}}\end{aligned}\right.$

The next step is to obtain time evolution equations for the collective
variables $x_{0}$ and $w$. To that end, we follow the procedure exposed in cs1
and susanto . After the analysis is carried over (details of the computations
are given in appendix A), we obtain the system of ordinary differential
equations:

$\displaystyle\dot{x_{0}}$
$\displaystyle=\frac{\alpha_{1}}{w}+\alpha_{2}wK_{0}(x_{0},w,s,R)+\alpha_{3}wK_{1}(x_{0},w,s,R)$
(9a) $\displaystyle\dot{w}$
$\displaystyle=\frac{\alpha_{4}}{w}+\alpha_{5}wK_{0}(x_{0},w,s,R)+\alpha_{6}wK_{1}(x_{0},w,s,R)$
(9b)

where the $\alpha_{i}$ are numbers and the integrals $K_{n}$ are given by

$K_{n}(x_{0},w,s,R)=\int_{-\infty}^{\infty}z^{n}s(wz+x_{0})~{}R(\tilde{U}(z))\tilde{U}^{\prime}(z)~{}dz.$
(10)

The integrals $K_{n}$ are the main driver of the front time evolution, as we
will see in the next sections.

###  III.1 Time-evolution equations for the collective variables

In the case of a bistable reaction term $R=R_{b}$, we are always able to
compute the integrals because there is unicity of the profile

$\tilde{U}(z)=\frac{1}{1+\mathrm{exp}(z)}$ (11)

and the equations (9) reduce to:

$\displaystyle\dot{x_{0}}$ $\displaystyle=-6wK_{0}(x_{0},w,s,R_{b}),$ (12a)
$\displaystyle\dot{w}$
$\displaystyle=\frac{3}{2(\pi^{2}-6)w}-\frac{18}{\pi^{2}-6}wK_{1}(x_{0},w,s,R_{b}),$
(12b)

These equations were obtained by Dawes and Susanto susanto .

In the case of a monostable reaction term $R=R_{m}$, we are able to explicitly
compute the integrals only when the profile $\tilde{U}$ is known. This is
especially the case if we consider the simplest following ansatz

$\tilde{U}(z)=\frac{1}{(1+(\sqrt{2}-1)\mathrm{exp}(z))^{2}}$ (13)

which is an exact solution of (1) when $s$ is constant. Then (9) give a system
silimar to (12) though a bit more intricate:

$\displaystyle\dot{x_{0}}$
$\displaystyle=\frac{\alpha_{1}}{w}+\alpha_{2}wK_{0}(x_{0},w,s,R_{m})+\alpha_{3}wK_{1}(x_{0},w,s,R_{m})$
(14a) $\displaystyle\dot{w}$
$\displaystyle=\frac{\alpha_{4}}{w}+\alpha_{5}wK_{0}(x_{0},w,s,R_{m})+\alpha_{6}wK_{1}(x_{0},w,s,R_{m}),$
(14b)

and the constants are (see appendix B for the definitions of the $I_{k}$):

$\displaystyle\alpha_{1}$
$\displaystyle=\frac{90I_{1}}{85-12\pi^{2}}\approx-0.026$
$\displaystyle\alpha_{2}$
$\displaystyle=\frac{900I_{2}}{85-12\pi^{2}}\approx-5.0$
$\displaystyle\alpha_{3}$ $\displaystyle=\frac{900I_{1}}{12\pi^{2}-85}\approx
0.26$ $\displaystyle\alpha_{4}$ $\displaystyle=\frac{18}{12\pi^{2}-85}\approx
0.54$ $\displaystyle\alpha_{5}$
$\displaystyle=\frac{900I_{1}}{12\pi^{2}-85}\approx 0.26$
$\displaystyle\alpha_{6}$ $\displaystyle=\frac{180}{85-12\pi^{2}}\approx-5.4$

###  III.2 Defect wide compared to the front width

For a defect that varies on a scale longer than $w$, $s(x_{0}+wz)\approx
s(x_{0})$, so that $s$ goes out of the integral, allowing us to take the
computation one step further. We obtain for the bistable case

$\left\\{\begin{aligned} \dot{x_{0}}&=\frac{1-2a}{2}ws(x_{0})\\\
\dot{w}&=\frac{3}{2(\pi^{2}-6)w}-\frac{3}{4(\pi^{2}-6)}ws(x_{0})\end{aligned}\right.$
(15)

and for the monostable case

$\left\\{\begin{aligned}
\left(\frac{12\pi^{2}-85}{90}\right)\dot{x_{0}}&=\frac{5-6\ln(1+\sqrt{2})}{30w}\\\
&+\frac{30\pi^{2}-220+9\ln(1+\sqrt{2})}{270}ws(x_{0})\\\
\left(\frac{12\pi^{2}-85}{90}\right)\dot{w}&=\frac{1}{5w}-\frac{1}{30}ws(x_{0})\\\
\end{aligned}\right.$ (16)

These equations justify the notion of a local speed and width of the front.
When $s(x)=s$ is constant, we recover the speeds and widths of the kinks
respectively (3) (6) and (4) (7) in the homogeneous situation.

For both the bistable and monostable cases, the right hand side of the
$\dot{x_{0}}$ equation is always positive so that no pinning of the front
occurs for wide defects. The numerical solutions match well the predictions
given by (15) and (16), see cs1 ; susanto  for the bistable case and bcp14
for the monostable case.

###  III.3 Defect narrow compared to the front width

We now turn our attention to defects that are narrow in comparison to the
width of the incident front. Namely, we consider here the case of a Heaviside
function. The defect function $s(x)$ can be written

$\boxed{s(x)=s_{0}+s_{1}H(x).}$ (17)

The integrals $K_{i}(x_{0},w,s,R)$ can be broken down into two parts:

$K_{i}(x_{0},w,s,R)=s_{0}\int_{-\infty}^{+\infty}+s_{1}\int_{\frac{-x_{0}}{w}}^{+\infty}$

We introduce the expressions

$\displaystyle K_{i}^{l}(u,R)$
$\displaystyle=\int_{-\infty}^{u}z^{i}R(\tilde{U}(z))\tilde{U}^{\prime}(z)~{}dz$
(18a) $\displaystyle K_{i}^{r}(u,R)$
$\displaystyle=\int_{u}^{+\infty}z^{i}R(\tilde{U}(z))\tilde{U}^{\prime}(z)~{}dz$
(18b)

with the exponents $l$ for left and $r$ for right. Then $K_{i}(x_{0},w,s,R)$
for the defect (17) can be written as

$K_{i}(x_{0},w,s,R)=s_{0}K_{i}(x_{0},w,1,R)+s_{1}K_{i}^{r}\left(\frac{-x_{0}}{w},R\right)$
(19)

Then the system (9) reads

$\displaystyle\dot{x_{0}}=\frac{\alpha_{1}}{w}$
$\displaystyle+ws_{0}\left(\alpha_{2}K_{0}(x_{0},w,1,R)+\alpha_{3}K_{1}(x_{0},w,1,R)\right)$
$\displaystyle+ws_{1}\left(\alpha_{2}K_{0}^{r}\left(\frac{-x_{0}}{w},R\right)+\alpha_{3}K_{1}^{r}\left(\frac{-x_{0}}{w},R\right)\right)$
(20a) $\displaystyle\dot{w}=\frac{\alpha_{4}}{w}$
$\displaystyle+ws_{0}\left(\alpha_{5}K_{0}(x_{0},w,1,R)+\alpha_{6}K_{1}(x_{0},w,1,R)\right)$
$\displaystyle+ws_{1}\left(\alpha_{5}K_{0}^{r}\left(\frac{-x_{0}}{w},R\right)+\alpha_{6}K_{1}^{r}\left(\frac{-x_{0}}{w},R\right)\right)$
(20b)

Up to here, everything is reaction-agnostic. We can, as we did earlier in
subsection III.1, compute the different terms in this system. It has been done
before for the bistable case susanto , and the constants for the monostable
one are given here in appendix B. We write them here for the sake of
completeness, first in the bistable case:

$\displaystyle\dot{x_{0}}$
$\displaystyle=\frac{1-2a}{2}ws_{0}-6ws_{1}K_{0}^{r}\left(\frac{-x_{0}}{w},R\right)$
(21a) $\displaystyle\dot{w}$
$\displaystyle=\frac{3}{2(\pi^{2}-6)w}+\frac{3}{4(6-\pi^{2})}ws_{0}$ (21b)
$\displaystyle+\frac{18}{6-\pi^{2}}ws_{1}K_{1}^{r}\left(\frac{-x_{0}}{w},R\right)$
(21c)

and second in the monostable case:

$\displaystyle\dot{x_{0}}$ $\displaystyle=-\frac{0.026}{w}+0.83ws_{0}$ (22a)
$\displaystyle-5ws_{1}K_{0}^{r}\left(\frac{-x_{0}}{w},R\right)+0.26ws_{1}K_{1}^{r}\left(\frac{-x_{0}}{w},R\right)$
(22b) $\displaystyle\dot{w}$ $\displaystyle=\frac{0.54}{w}-0.09ws_{0}$ (22c)
$\displaystyle+0.26ws_{1}K_{0}^{r}\left(\frac{-x_{0}}{w},R\right)-5.4ws_{1}K_{1}^{r}\left(\frac{-x_{0}}{w},R\right)$
(22d)

To illustrate how well the collective variables match with the PDE solution we
consider a Heaviside defect (9) such that

$s(x)=0.2,~{}~{}~{}x<0~{}~{}~{}~{}~{}~{}s(x)=0.1,~{}~{}~{}x>0.$

Fig. 1 shows $x_{0}(t),~{}w(t),~{}w(x_{0})$ from left to right for both the
monostable PDE and equations (20). The values $x_{0},w$ are obtained by
fitting the PDE solution by the ansatz (13) just as was done for the bistable
solution in cs1 . One can see the good agreement between the two solutions.

Figure 1: Motion of a monostable front in a tanh like defect, $s=0.2$ for
negative $x$ and $s=0.1$ for positive $x$. From left to right:
$x_{0}(t),~{}w(t),~{}w(x_{0})$.

###  III.4 Reaction-free zone

We consider here a reaction-free zone. This occurs for example in forest fires
where a trench with no trees is made to prevent fire propagation. In
combustion in a duct, this corresponds to a region where there is no fuel. The
defect function $s(x)$ can be written

$\boxed{s(x)=0,~{}~{}~{}0\leq x\leq d,~{}~{}~{}s(x)=s_{0}~{}~{}~{}{\rm
elsewhere}.}$ (23)

With that defect, the integrals $K_{i}(x_{0},w,s,R)$ break down into two
parts:

$K_{i}(x_{0},w,s,R)=\int_{-\infty}^{\frac{-x_{0}}{w}}+\int_{\frac{d-x_{0}}{w}}^{+\infty}$

with the naming conventions of (18). We get

$K_{i}(x_{0},w,s,R)=s_{0}K_{i}^{l}\left(\frac{-x_{0}}{w},R\right)+s_{0}K_{i}^{r}\left(\frac{d-x_{0}}{w},R\right)$
(24)

so that the system (9) reads

$\displaystyle\dot{x_{0}}=\frac{\alpha_{1}}{w}$
$\displaystyle+ws_{0}\left(\alpha_{2}K_{0}^{l}\left(\frac{-x_{0}}{w},R\right)+\alpha_{3}K_{1}^{l}\left(\frac{-x_{0}}{w},R\right)\right)$
$\displaystyle+ws_{0}\left(\alpha_{2}K_{0}^{r}\left(\frac{d-x_{0}}{w},R\right)+\alpha_{3}K_{1}^{r}\left(\frac{d-x_{0}}{w},R\right)\right)$
(25a) $\displaystyle\dot{w}=\frac{\alpha_{4}}{w}$
$\displaystyle+ws_{0}\left(\alpha_{5}K_{0}^{l}\left(\frac{-x_{0}}{w},R\right)+\alpha_{6}K_{1}^{l}\left(\frac{-x_{0}}{w},R\right)\right)$
$\displaystyle+ws_{0}\left(\alpha_{5}K_{0}^{r}\left(\frac{d-x_{0}}{w},R\right)+\alpha_{6}K_{1}^{r}\left(\frac{d-x_{0}}{w},R\right)\right)$
(25b)

Again, up to here, the computation doesn’t depend on the shape of the reaction
term $R$.

##  IV Reaction-free zone for bistable

In the bistable case, the major result is that the front can be stopped for a
wide enough reaction-free zone. We first illustrate this on a 2D example and
then go on to analyze it in 1D.

Consider the 2D reaction-diffusion equation

$u_{t}=u_{xx}+u_{yy}+s(x)~{}R(u)$ (26)

in a rectangular domain $|x|\leq 150,~{}~{}~{}|y|\leq 150$ with homogeneous
Neuman boundary conditions. We assume the bistable nonlinearity
$R(u)=u(1-u)(a-u)$ with $a=0.3$. The reaction free region is

$0\leq x\leq 40.$

The initial condition is a square pulse centered at $x=-20,~{}~{}y=0$ of
extension 12.

Figure 2: Motion of a two-dimensional bistable front reaching the reaction-
free region $0\leq x\leq 40$. The four panels show the contour lines of the
solution for increasing times.

The time evolution of $u(x,y,t)$ was computed using the Comsol software comsol
, snapshots are shown in Fig. 2 as contour lines ranging from $0$ (dark blue
online) to $1$ (red online). After having assumed a radial shape, the pulse
hits the defect. It then extends along it in the $y$ direction. There is not
much difference in the region $x\geq 0$ between the profiles at times $t=130$
and $250$. This indicates that the wave cannot enter the region $x\geq 0$.

To understand this effect quantitatively, we return to the 1D case. The
analysis is developed in the following subsections.

###  IV.1 Criterion for front stopping

In this section, we analyze the stopping of the front. For simplicity of the
analysis, we assume $s_{0}=1$ because the amplitude $s_{0}$ of the defect
$s(x)$ can be scaled out of the problem by the following change of variables

$x={x^{\prime}\over\sqrt{s_{0}}},~{}~{}t={t^{\prime}\over
s_{0}},~{}~{}s=s_{0}s^{\prime}(x^{\prime}).$ (27)

When the front stops, $x_{0}$ and $w$ are stationary. The collective variable
odes (25) yield

$\displaystyle-\frac{\alpha_{1}}{w}$
$\displaystyle=ws_{0}\left(\alpha_{2}K_{0}^{l}\left(\frac{-x_{0}}{w},R\right)+\alpha_{3}K_{1}^{l}\left(\frac{-x_{0}}{w},R\right)\right)$
$\displaystyle+ws_{0}\left(\alpha_{2}K_{0}^{r}\left(\frac{d-x_{0}}{w},R\right)+\alpha_{3}K_{1}^{r}\left(\frac{d-x_{0}}{w},R\right)\right)$
(28a) $\displaystyle-\frac{\alpha_{4}}{w}$
$\displaystyle=ws_{0}\left(\alpha_{5}K_{0}^{l}\left(\frac{-x_{0}}{w},R\right)+\alpha_{6}K_{1}^{l}\left(\frac{-x_{0}}{w},R\right)\right)$
$\displaystyle+ws_{0}\left(\alpha_{5}K_{0}^{r}\left(\frac{d-x_{0}}{w},R\right)+\alpha_{6}K_{1}^{r}\left(\frac{d-x_{0}}{w},R\right)\right)$
(28b)

then

$\displaystyle-\frac{\alpha_{1}\alpha_{4}}{w^{2}s_{0}}$
$\displaystyle=\alpha_{2}\alpha_{4}K_{0}^{l}\left(\frac{-x_{0}}{w},R\right)+\alpha_{3}\alpha_{4}K_{1}^{l}\left(\frac{-x_{0}}{w},R\right)$
$\displaystyle+\alpha_{2}\alpha_{4}K_{0}^{r}\left(\frac{d-x_{0}}{w},R\right)+\alpha_{3}\alpha_{4}K_{1}^{r}\left(\frac{d-x_{0}}{w},R\right)$
(29a) $\displaystyle-\frac{\alpha_{1}\alpha_{4}}{w^{2}s_{0}}$
$\displaystyle=\alpha_{1}\alpha_{5}K_{0}^{l}\left(\frac{-x_{0}}{w},R\right)+\alpha_{1}\alpha_{6}K_{1}^{l}\left(\frac{-x_{0}}{w},R\right)$
$\displaystyle+\alpha_{1}\alpha_{5}K_{0}^{r}\left(\frac{d-x_{0}}{w},R\right)+\alpha_{1}\alpha_{6}K_{1}^{r}\left(\frac{d-x_{0}}{w},R\right)$
(29b)

The system implies an equation that can be written as

$F^{l}\left(\frac{x_{0}}{w},R\right)=F^{r}\left(\frac{d-x_{0}}{w},R\right),$
(30)

Equation (30) is verified for a critical distance

$d=d_{c}.$

The feasibility of pinning can then easily be checked by studying the
functions $F^{l}$ and $F^{r}$; these are plotted in Fig. 3 for the bistable
nonlinearity.

Figure 3: Left: the functions $F^{l}(u,R_{b})$ and $F^{r}(u,R_{b})$. Right:
zoom on the $[-1,4]$ region.

For the critical distance $d_{c}$, we need the maximum of $F_{r}$ to cancel
out the contribution from $F_{l}$. If we define $x_{r}=\rm{Argmax}(F_{r})$,
$x_{l}=F_{l}^{-1}(\max(F_{r}))$, then the following equations hold:

$\displaystyle\frac{d_{c}-x_{0}}{w}=x_{r}$ (31a)
$\displaystyle\frac{x_{0}}{w}=x_{l}$ (31b)

They give in turn $w$ using (28b) (a bit of caution is necessary as
$\alpha_{1}=0$ in the bistable case) and finally $d_{c}$. As an illustration,
we plot in Fig. 4 the width of the front when pinning occurs (left panel), and
the critical distance (right panel), as functions of the parameter $a$ of the
bistable reaction term. When the parameter $a=0.3$ for instance, we have a
critical distance of $6.32$ and a width at pinning of $2.19$. They compare
very well with what we obtain through a simulation of the original PDE.

Figure 4: Left: width of the front when pinning occurs as a function of the
parameter $a$. Right: critical distance of the non-reaction zone as a function
of the parameter $a$.

###  IV.2 Stationary solution when the front stops

Let us consider the stationary case. The problem is

$\displaystyle u_{xx}+u(1-u)(u-a)$ $\displaystyle=0,~{}~{}x\leq 0,$ (32)
$\displaystyle u_{xx}$ $\displaystyle=0,~{}~{}0\leq x\leq d,$ (33)
$\displaystyle u_{xx}+u(1-u)(u-a)$ $\displaystyle=0,~{}~{}x\geq d,$ (34)

together with the boundary conditions $u=1,~{}u_{x}=0,~{}~{}x\to-\infty$ and
$u=0,~{}u_{x}=0,~{}~{}x\to\infty$. The interface conditions at $x=0,d$ are
continuity of $u$ and of $u_{x}$. Multiplying (32,34) by $u_{x}$ and
integrating, we get

${u_{x}^{2}\over 2}+\left(-{u^{4}\over 4}+{u^{3}\over 3}(a+1)-{u^{2}\over
2}a\right)=C_{1},~{}~{}~{}x\leq 0.$ ${u_{x}^{2}\over 2}+\left(-{u^{4}\over
4}+{u^{3}\over 3}(a+1)-{u^{2}\over 2}a\right)=0,~{}~{}~{}x\geq d.$

The constant $C_{1}$ is obtained from the first expression evaluated for
$x\to-\infty$ where $u=1$ and $u_{x}=0$. We get

$C_{1}={1-2a\over 12}.$ (35)

Inside the strip $0\leq x\leq d$, the solution is

$u=\alpha x+\beta.$

Using the interface conditions at $x=0,d$, we get the two relations for the
unknowns $\alpha,\beta$.

$\displaystyle{\alpha^{2}\over 2}$ $\displaystyle-{1\over 4}\beta^{4}+{1\over
3}(a+1)\beta^{3}-{a\over 2}\beta^{2}=C_{1},$ (36)
$\displaystyle{\alpha^{2}\over 2}$ $\displaystyle-{1\over 4}(\alpha
d+\beta)^{4}+{1\over 3}(a+1)(\alpha d+\beta)^{3}-{a\over 2}(\alpha
d+\beta)^{2}=0.$ (37)

To solve this algebraic system we used a graphical method, introducing the
values $u_{l}=\beta$ and $u_{r}=\alpha d+\beta$. In Fig. 5 we plot the $0$
contour lines of the functions $g(u_{l},u_{r})$ and $h(u_{l},u_{r})$ defined
by

$\displaystyle g(u_{l},u_{r})={(u_{l}-u_{r})^{2}\over 2d^{2}}+f(u_{l})-C_{1},$
(38) $\displaystyle h(u_{l},u_{r})={(u_{l}-u_{r})^{2}\over 2d^{2}}+f(u_{r}),$
(39) $\displaystyle f(x)=-{1\over 4}x^{4}+{a+1\over 3}x^{3}-{a\over 2}x^{2},$
(40)

in the square $(u_{l},u_{r})\in[0;1]\times[0;1]$.

Figure 5: Bistable case, zero reaction region $0\leq x\leq d$. Plots of the
functions $g(u_{l},u_{r})$ (curve) and $h(u_{l},u_{r})$ (line) in the
$(u_{l},u_{r})$ plane for two different values of $d$, $d=6<d_{c}$ (left) and
$d=7>d_{c}$ (right); the parameter is $a=0.3$.

In the left panel $d<d_{c}$ so that there is no solution. In the right panel
$d>d_{c}$ and there are 2 solutions $(u_{l},u_{r})$. One of them is stable and
the other unstable, as expected from standard bifurcation theory. For $a=0.3$
we find $d_{c}\approx 6.5$ close to the value obtained from the pinning
criterion of section IV.1.

To compare with these estimates, we analyze the PDE solutions for two
different values of $d$. Fig. 6 shows snapshots of the solution $u(x,t)$ for
different times for $d=6$ (left panel) and $d=7$ (right panel). As expected
the front is stopped by the zero reaction region for $d>6.5$. Therefore the
PDE solution agrees with the analytical estimates.

Figure 6: Snapshots of $u$ for $t=50,100,150,200,250$ and $300$ when the front
crosses the zero reaction region, $d=6$ (left panel) and when the front is
stopped by the zero reaction region, $d=7$ (right panel).

##  V Reaction-free zone for monostable

###  V.1 Non existence of stationary solution

Following the same strategy as for the bistable case, we look for a stationary
solution of the Fisher problem with a strip $[0;d]$. The formalism is the same
as above. We then get

${u_{x}^{2}\over 2}+\left({u^{2}\over 2}-{u^{3}\over 3}\right)={s\over
6},x\leq 0.$ ${u_{x}^{2}\over 2}+\left({u^{2}\over 2}-{u^{3}\over
3}\right)=0,x\geq d.$

There is no solution to the problem because the second expression is always
greater than 0, for $u\in[0;1]$. This is consistent with the instability of
the $u=0$ stationary state for the monostable case.

###  V.2 Appearance of secondary front

As for the bistable model, a 2D simulation is useful to illustrate what is
happening. We consider the 2D reaction-diffusion equation (26) with the
monostable nonlinearity. Two snapshots are shown in Fig. 7, $t=90$ (left
panel) and $t=100$ (right panel). For $t=90$, the front seems to be trapped at
the left interface $x=0$ of the reaction-free zone.

Figure 7: Monostable case, zero reaction region $0\leq x\leq 40$ and
$s_{0}=0.3$. Contour lines of the solution at times $t=90$ (left) and $t=100$
(right).

However, a new structure appears close to the right interface $x=40$. The
right panel of Fig. 7 shows that this new structure is a secondary pulse.

The 1D problem can be analyzed in more details. Fig. shows snapshots $u(x,t)$
of the 1D solution for seven successive times from $t=27.5$ up to 192.5 for a
reaction-free region $|x|\leq 35$. The left panel shows the solution in linear
scale. One can see the advancing profile for small $u$. It reaches the right
interface and becomes visible at $t=165$.

Figure 8: Motion of a monostable front in a zero defect region of extension
$d=70$: snapshots of $u(x,t)$ for $t=27.5,55,82.5,110,137.5,165$ and 192.5 for
a no reaction zone of extension $d=70$ in linear (left) and log (right) scales
for $u$.

The increase of $u$ with time is clearer in the right panel (log plot). There,
one sees the diffusion of $u$ in the reaction-free zone followed by its
amplification for $x\geq 35$.

The time of appearance of the secondary front can be estimated from the
solution. Assume a zero reaction-region

$0\leq x\leq d.$

In Fig. 9 we plot the time interval $t_{2}-t_{1}$ as a function of the half-
width $d$ of the zero reaction region. The instant $t_{1}$ corresponds to the
front reaching the left edge of the strip, when its position is $x=0$. The
time $t_{2}$ is when the secondary pulse reaches the value $0.5$.

Figure 9: Time of appearance of secondary front $t_{2}-t_{1}$ for the
monostable case for different widths of the zero reaction region. The blue
line corresponds to the numerical results and the black straight line is for
the approximation (44).

It is surprising that $t_{2}-t_{1}$ depends linearly on $d$. To explain this
effect, we suggest the following simple model. In the linear (no reaction)
zone, the propagation is governed by the heat kernel giving the solution at
$(x,t)$ for a source at $(y,t=0)$

$K(x,y,t)={1\over\sqrt{4\pi t}}e^{-{(x-y)^{2}\over 4t}}.$ (41)

Then for $x=d$, the solution of the heat equation $u_{h}$ is
$u_{h}(d,t)={1\over\sqrt{4\pi t}}e^{-d^{2}\over 4t}$. The $u=0$ solution of
the Fisher equation is unstable and any plane wave perturbation $e^{ikx}$ will
grow as

$\delta u=e^{ikx+(s_{0}-k^{2})t}.$ (42)

Combining (41) and (42), the solution of the Fisher equation at $x=d$ will be

$u(d,t)={1\over\sqrt{4\pi t}}e^{-{d^{2}\over 4t}+s_{0}t},$ (43)

where $\xi$ is a coefficient smaller than 1 accounting for the average of the
growth rate.

We can calculate when $u(d,t)=0.5$. Using (43) we obtain

${1\over\sqrt{4\pi t}}e^{-{d^{2}\over 4t}+s_{0}t}=0.5.$

This implies

$d=\sqrt{4s_{0}t^{2}-2t\log(\pi t)}$ (44)

In Fig. 9 we plotted the values of $d,t$ such that $u(d,t)=0.5$. These follow
a straight line (red) that is not far from to the numerical results (blue).
Despite appearances, the formula (44) gives a linear dependance because of the
interval in t $[20,180]$ on which it is plotted.

##  VI Conclusion

We analyzed the stopping of a 1D reaction-diffusion front by a reaction-free
region for bistable and monostable nonlinearities.

Bistable fronts can be stopped and -using a collective variable description-
we obtained a stopping criterion linking the width of the region and the
parameter of the nonlinearity. This criterion is in good agreement with the
analysis of the time independant problem and the PDE solutions.

The monostable nonlinearity is more complex. If the front is accelerated by
the defect, the collective variables agree well with the PDE numerical
solution, otherwise a secondary pulse appears and the collective variable
description is wrong.

For a reaction-free region, we can predict the time of appearance of the
secondary pulse using a simple model based on the diffusion kernel and the
growth rate of the zero unstable state.

This reaction-free region can be generalized into a damped-reaction region to
model the action of (chemo/radio) therapy on a cancer tumor cancer2 . The
function $s(x)$ becomes

$s(x)=-s_{1},~{}~{}~{}0\leq x\leq d,~{}~{}~{}s(x)=s_{0}~{}~{}~{}{\rm
elsewhere}.$ (45)

Using the same approach as for (43), we find the form of the solution in the
damped-reaction region as

$u(d,t)={1\over 4\pi t}e^{{d^{2}\over 4t}+(s_{0}-s_{1})t},$ (46)

so that we can predict the time of crossing

$d=\sqrt{4(s_{0}-s_{1})t^{2}-2t\log(\pi t)}.$ (47)

The radically different behaviors between a monostable and a bistable front
raises important questions for modelers. Does it make sense that a monostable
front front can cross a no reaction zone of arbitrary width? Maybe tumor
researchers should use a $u^{2}(1-u)$ nonlinearity instead of the standard
logistic. Then $u=0$ would be stable and there would be a critical distance
that the front could not cross.

##  Appendix A Derivation of the time evolution equations for the collective
variables

We multiply the equation

$\dot{x_{0}}\frac{\partial U}{\partial x_{0}}+\dot{w}\frac{\partial
U}{\partial w}=\frac{\partial^{2}U}{\partial
x^{2}}+s(x)~{}R(U(x,x_{0}(t),w(t)))$

by the test functions $\frac{\partial U}{\partial x_{0}}$ and $\frac{\partial
U}{\partial w}$ respectively and we integrate over $x$ going from $-\infty$ to
$+\infty$. We obtain

$\left\\{\begin{aligned}
&{}\dot{x_{0}}\int_{-\infty}^{\infty}\left(\frac{\partial U}{\partial
x_{0}}\right)^{2}+\dot{w}\int_{-\infty}^{\infty}\frac{\partial U}{\partial
x_{0}}\frac{\partial U}{\partial w}\\\
&=\int_{-\infty}^{\infty}\frac{\partial^{2}U}{\partial x^{2}}\frac{\partial
U}{\partial
x_{0}}+\int_{-\infty}^{\infty}s(x)~{}R(U(x,x_{0}(t),w(t)))\frac{\partial
U}{\partial x_{0}}\\\ &{}\dot{x_{0}}\int_{-\infty}^{\infty}\frac{\partial
U}{\partial x_{0}}\frac{\partial U}{\partial
w}+\dot{w}\int_{-\infty}^{\infty}\left(\frac{\partial U}{\partial
w}\right)^{2}\\\ &=\int_{-\infty}^{\infty}\frac{\partial^{2}U}{\partial
x^{2}}\frac{\partial U}{\partial
w}+\int_{-\infty}^{\infty}s(x)~{}R(U(x,x_{0}(t),w(t)))\frac{\partial
U}{\partial w}\end{aligned}\right.$

Thanks to the chain rule, we have

$\left\\{\begin{aligned} \frac{\partial U}{\partial
x}&=\frac{1}{w}\tilde{U}^{\prime}(z)\quad\mbox{and}\quad\frac{\partial^{2}U}{\partial
x^{2}}=\frac{1}{w^{2}}\tilde{U}^{\prime\prime}(z)\\\ \frac{\partial
U}{\partial x_{0}}&=-\frac{1}{w}\tilde{U}^{\prime}(z)\\\ \frac{\partial
U}{\partial w}&=-\frac{z}{w}\tilde{U}^{\prime}(z)\\\ \end{aligned}\right.$

We remark that it is convenient to make the change of variable
$z=\frac{x-x_{0}(t)}{w(t)}$ in each integral to get (after multiplication by
$w$):

$\left\\{\begin{aligned}
&{}\dot{x_{0}}\int_{-\infty}^{\infty}\tilde{U}^{\prime}(z)^{2}~{}+\dot{w}\int_{-\infty}^{\infty}z\tilde{U}^{\prime}(z)^{2}\\\
&=-\frac{1}{w}\int_{-\infty}^{\infty}\tilde{U}^{\prime}(z)\tilde{U}^{\prime\prime}(z)-w\int_{-\infty}^{\infty}s(wz+x_{0})~{}R(\tilde{U}(z))\tilde{U}^{\prime}(z)\\\
&{}\dot{x_{0}}\int_{-\infty}^{\infty}z\tilde{U}^{\prime}(z)^{2}+\dot{w}\int_{-\infty}^{\infty}z^{2}\tilde{U}^{\prime}(z)^{2}\\\
&=-\frac{1}{w}\int_{-\infty}^{\infty}z\tilde{U}^{\prime}(z)\tilde{U}^{\prime\prime}(z)-w\int_{-\infty}^{\infty}zs(wz+x_{0})~{}R(\tilde{U}(z))\tilde{U}^{\prime}(z)\\\
\end{aligned}\right.$

The system is written in abstract form:

$\left\\{\begin{aligned}
I_{0}\dot{x_{0}}+I_{1}\dot{w}&=-\frac{J_{0}}{w}-wK_{0}(x_{0},w,s,R)\\\
I_{1}\dot{x_{0}}+I_{2}\dot{w}&=-\frac{J_{1}}{w}-wK_{1}(x_{0},w,s,R)\\\
\end{aligned}\right.$

where

$\left\\{\begin{aligned}
I_{n}&=\int_{-\infty}^{\infty}z^{n}\tilde{U}^{\prime}(z)^{2}~{}dz\\\
J_{n}&=\int_{-\infty}^{\infty}z^{n}\tilde{U}^{\prime}(z)\tilde{U}^{\prime\prime}(z)~{}dz\\\
K_{n}(x_{0},w,s,R)&=\int_{-\infty}^{\infty}z^{n}s(wz+x_{0})~{}R(\tilde{U}(z))\tilde{U}^{\prime}(z)~{}dz\end{aligned}\right.$

The integrals $I_{n}$ and $J_{n}$ are numbers, they only depend on the profile
$\tilde{U}$. In the cases that we consider, the determinant of the system is
non zero and we get the following equations governing $\dot{x_{0}}$ and
$\dot{w}$:

$\left\\{\begin{aligned}
\dot{x_{0}}&=\frac{\alpha_{1}}{w}+\alpha_{2}wK_{0}(x_{0},w,s,R)+\alpha_{3}wK_{1}(x_{0},w,s,R)\\\
\dot{w}&=\frac{\alpha_{4}}{w}+\alpha_{5}wK_{0}(x_{0},w,s,R)+\alpha_{6}wK_{1}(x_{0},w,s,R)\end{aligned}\right.$

where

$\left\\{\begin{aligned}
\alpha_{1}&=\frac{I_{1}J_{1}-I_{2}J_{0}}{I_{0}I_{2}-I_{1}^{2}}\\\
\alpha_{2}&=\frac{I_{2}}{I_{1}^{2}-I_{0}I_{2}}\\\
\alpha_{3}&=\frac{I_{1}}{I_{0}I_{2}-I_{1}^{2}}\\\
\alpha_{4}&=\frac{I_{1}J_{0}-I_{0}J_{1}}{I_{0}I_{2}-I_{1}^{2}}\\\
\alpha_{5}&=\frac{I_{1}}{I_{0}I_{2}-I_{1}^{2}}\\\
\alpha_{6}&=\frac{I_{0}}{I_{1}^{2}-I_{0}I_{2}}\\\ \end{aligned}\right.$

##  Appendix B Exact values of integrals present in the paper

For the bistable case,

$\displaystyle I_{0}$
$\displaystyle=\int_{-\infty}^{\infty}\tilde{U}^{\prime}(z)^{2}~{}dz=\frac{1}{6}$
$\displaystyle I_{1}$
$\displaystyle=\int_{-\infty}^{\infty}z\tilde{U}^{\prime}(z)^{2}~{}dz=0$
$\displaystyle I_{2}$
$\displaystyle=\int_{-\infty}^{\infty}z^{2}\tilde{U}^{\prime}(z)^{2}~{}dz=\frac{\pi^{2}-6}{18}\approx
0.21$ $\displaystyle J_{0}$
$\displaystyle=\int_{-\infty}^{\infty}\tilde{U}^{\prime}(z)\tilde{U}^{\prime\prime}(z)~{}dz=0$
$\displaystyle J_{1}$
$\displaystyle=\int_{-\infty}^{\infty}z\tilde{U}^{\prime}(z)\tilde{U}^{\prime\prime}(z)~{}dz=-\frac{1}{12}$
$\displaystyle K_{0}(x_{0},w,1,R_{b})$
$\displaystyle=\int_{-\infty}^{\infty}R_{b}(\tilde{U}(z))\tilde{U}^{\prime}(z)~{}dz=\frac{2a-1}{12}$
$\displaystyle K_{1}(x_{0},w,1,R_{b})$
$\displaystyle=\int_{-\infty}^{\infty}zR_{b}(\tilde{U}(z))\tilde{U}^{\prime}(z)~{}dz=\frac{1}{24}$

and the constants in (9) read

$\displaystyle\alpha_{1}$ $\displaystyle=0$ $\displaystyle\alpha_{2}$
$\displaystyle=-6$ $\displaystyle\alpha_{3}$ $\displaystyle=0$
$\displaystyle\alpha_{4}$ $\displaystyle=\frac{1}{12I_{2}}\approx 0.39$
$\displaystyle\alpha_{5}$ $\displaystyle=0$ $\displaystyle\alpha_{6}$
$\displaystyle=-\frac{1}{I_{2}}\approx-4.7$

For the monostable case,

$\displaystyle I_{0}$
$\displaystyle=\int_{-\infty}^{\infty}\tilde{U}^{\prime}(z)^{2}~{}dz=\frac{1}{5}$
$\displaystyle I_{1}$
$\displaystyle=\int_{-\infty}^{\infty}z\tilde{U}^{\prime}(z)^{2}~{}dz=\frac{\ln(1+\sqrt{2})}{5}-\frac{1}{6}$
$\displaystyle\approx 0.0096$

$\displaystyle I_{2}$
$\displaystyle=\int_{-\infty}^{\infty}z^{2}\tilde{U}^{\prime}(z)^{2}~{}dz$
$\displaystyle=-\frac{1}{3}+\frac{\pi^{2}}{15}-\frac{\ln(1+\sqrt{2})}{3}+\frac{(\ln(1+\sqrt{2}))^{2}}{5}$
$\displaystyle\approx 0.19$ $\displaystyle J_{0}$
$\displaystyle=\int_{-\infty}^{\infty}\tilde{U}^{\prime}(z)\tilde{U}^{\prime\prime}(z)~{}dz=0$
$\displaystyle J_{1}$
$\displaystyle=\int_{-\infty}^{\infty}z\tilde{U}^{\prime}(z)\tilde{U}^{\prime\prime}(z)~{}dz=-\frac{1}{10}$
$\displaystyle K_{0}(x_{0},w,1,R_{m})$
$\displaystyle=\int_{-\infty}^{\infty}R_{m}(\tilde{U}(z))\tilde{U}^{\prime}(z)~{}dz=-\frac{1}{6}$
$\displaystyle K_{1}(x_{0},w,1,R_{m})$
$\displaystyle=\int_{-\infty}^{\infty}zR_{m}(\tilde{U}(z))\tilde{U}^{\prime}(z)~{}dz$
$\displaystyle=\frac{7}{45}-\frac{\ln(1+\sqrt{2})}{6}\approx 0.0087$

and the constants in (9) read

$\displaystyle\alpha_{1}$
$\displaystyle=\frac{90I_{1}}{85-12\pi^{2}}\approx-0.026$
$\displaystyle\alpha_{2}$
$\displaystyle=\frac{900I_{2}}{85-12\pi^{2}}\approx-5.0$
$\displaystyle\alpha_{3}$ $\displaystyle=\frac{900I_{1}}{12\pi^{2}-85}\approx
0.26$ $\displaystyle\alpha_{4}$ $\displaystyle=\frac{18}{12\pi^{2}-85}\approx
0.54$ $\displaystyle\alpha_{5}$
$\displaystyle=\frac{900I_{1}}{12\pi^{2}-85}\approx 0.26$
$\displaystyle\alpha_{6}$ $\displaystyle=\frac{180}{85-12\pi^{2}}\approx-5.4$

###### Acknowledgements.

Part of this work was performed using computing resources of CRIANN (Normandy,
France).

## References

  * (1) Ya.B. Zeldovich and D.A. Frank-Kamenetsky, K teorii ravnomernogo rasprostraneniya plameni, Dokladi Akademii Nauk SSSR, Vol. 19(9):693-697 (1938).
  * (2) A. Scott, Nonlinear science, emergence and dynamics of coherent structures, Oxford University Press (2nd edition, 2003).
  * (3) J. Xin, Front propagation in heterogeneous media. SIAM Review, 42, 161–230, (2000).
  * (4) K. R. Swanson, R. Rostomily, E. C. Alvord Jr, Predicting survival of patients with glioblastoma by combining a mathematical model and pre-operative MR imaging characteristics: a proof of principle. British J. Cancer , 98, 113–9, (2008).
  * (5) J. D. Murray, Mathematical biology I: an introduction. 3rd ed. New York: Springer; 2002.
  * (6) R. Fisher, The wave of advance of advantageous genes, Annals of Eugenics, Vol. 7:355-369 (1937).
  * (7) M. J. Ablowitz and A. Zeppetella, Explicit solutions of Fisher’s equation for a special wave speed, Bulletin of Mathematical Biology, Vol. 41:835-840 (1979).
  * (8) J.-G. Caputo and B. Sarels, Reaction-diffusion front crossing a local defect, Phys. Rev. E 84, 041108 (2011).
  * (9) J. H. P. Dawes and H. Susanto, Variational approximation and the use of collective coordinates, Phys. Rev. E 87, 063202 (2013).
  * (10) J. Belmonte-Beitia, G. F. Calvo, V. M. Perez-Garcia, Effective particle methods for Fisher–Kolmogorov equations: Theory and applications to brain tumor dynamics Commun Nonlinear Sci Numer Simulat 19 (2014) 3267–3283.
  * (11) H. Mi, C. Petitjean, B. Dubray, P. Vera, and S. Ruan, Prediction of Lung Tumor Evolution During Radiotherapy in Individual Patients With PET, IEEE transactions on medical imaging, vol. 33, 4, (2014).
  * (12) H. Berestycki, J. Bouhours and G. Chapuisat, ”Front blocking and propagation in cylinders with varying cross section”, Calculus of Variations and Partial Differential Equations, Springer, (2016).
  * (13) Juliette Bouhours, ”Reaction diffusion equation in heterogeneous media : persistance, propagation and effect of the geometry”, General Mathematics, Université Pierre et Marie Curie - Paris VI, (2014).
https://tel.archives-ouvertes.fr/tel-01070608

  * (14) J.M. Roquejoffre, Convergence to Travelling Waves for Solutions of a Class of Semilinear Parabolic Equations, Journal of Differential Equations, Vol. 108(2):262-295, (1994).
  * (15) J. Garnier and T. Giletti and F. Hamel and L. Roques, Inside dynamics of pulled and pushed fronts, Journal de Mathématiques Pures et Appliquées, Vol. 98(4):428-449 (2012).
  * (16) https://www.comsol.com/