# Emergence of disconnected clusters in heterogeneous complex systems

István A. Kovács<EMAIL_ADDRESS>Northwestern University,
Department of Physics and Astronomy, Evanston, 60208, USA  Wigner Research
Centre for Physics, Institute for Solid State Physics and Optics, Budapest,
1121, Hungary  Róbert Juhász  Wigner Research Centre for Physics, Institute
for Solid State Physics and Optics, Budapest, 1121, Hungary

###### Abstract

Percolation theory dictates an intuitive picture depicting correlated regions
in complex systems as densely connected clusters. While this picture might be
adequate at small scales and apart from criticality, we show that highly
correlated sites in complex systems can be inherently disconnected. This
finding indicates a counter-intuitive organization of dynamical correlations,
where functional similarity decouples from physical connectivity. We
illustrate the phenomena on the example of the Disordered Contact Process
(DCP) of infection spreading in heterogeneous systems. We apply numerical
simulations and an asymptotically exact renormalization group technique (SDRG)
in $1,2$ and $3$ dimensional systems as well as in two-dimensional lattices
with long-ranged interactions. We conclude that the critical dynamics is well
captured by mostly one, highly correlated, but spatially disconnected cluster.
Our findings indicate that at criticality the relevant, simultaneously
infected sites typically do not directly interact with each other. Due to the
similarity of the SDRG equations, our results hold also for the critical
behavior of the disordered quantum Ising model, leading to quantum correlated,
yet spatially disconnected, magnetic domains.

## Introduction

Correlated clusters emerge in a broad range of systems, ranging from magnetic
models to out-of-equilibrium systems FK ; CI ; MD . Apart from the critical
point, the correlation length is finite, limiting the spatial separation
between highly correlated sites, leading to spatially localized, finite
clusters. Although at the critical point the correlation length diverges, our
traditional intuition is driven by percolation processes, indicating that
highly correlated critical clusters remain connected, while broadly varying in
their size. As an alternative scenario, we show that at the critical point,
highly correlated clusters can break up into spatially disconnected regions in
both classical and quantum systems. As a proof-of-concept example, we focus on
the Disordered Contact Process (DCP), a simple (out-of-equilibrium) infection
spreading model noest ; noest2 ; moreira . The critical behavior of the model
is well understood in the presence of disorder, at least at the level of
statistical properties averaged over a large numbers of samples moreira ; vd ;
vfm ; vojta_3d ; vojta_rev . The detailed simulation of individual samples is
much more challenging, due to large dynamical fluctuations and extremely slow
dynamics around criticality hiv ; hiv2 .

As an alternative, efficient approach, the strong disorder renormalization
group (SDRG) method provides asymptotically exact results for the critical
contact process in the presence of disorder, at least below $d=4$ spatial
dimensions hiv ; hiv2 ; im ; im2 . Besides being computationally efficient,
the SDRG method also offers some counter-intuitive insights into the
underlying correlation structure. Namely, the SDRG predicts that highly
correlated sites at long time scales (simultaneously infected individuals)
typically do not know each other directly, only via indirect connections
through the rest of the contact network. According to the SDRG, only a few of
these highly correlated, but essentially disconnected, clusters govern the
large-scale behavior of the system. Observing these clusters in simulations is
notoriously difficult due to extremely slow dynamics, inducing anomalously
large fluctuations and poor statistical properties. As a key step, here we
show how to find the highly correlated sites efficiently via a quasi-
stationary simulation. The simulated density profile is then found to be in a
good agreement with the SDRG predictions, confirming that the asymptotic
dynamics is governed by spatially disconnected clusters in stark contrast to
traditional intuition.

### The disordered contact process

The contact process, in the most general case, is defined on a network given
by an adjacency matrix $A_{ij}$. As a special case, the network is often
chosen to be a d-dimensional hypercubic lattice with nearest-neighbour edges.
The state of the system is given by a set of binary variables, $n_{i}=0,1$,
characterizing the sites of the network. As the contact process is frequently
interpreted as a simple epidemic spreading model, sites with $n_{i}=1$ are
referred as ’infected’, while sites with $n_{i}=0$ are ’healthy’ or
’susceptible’. The contact process is a continuous-time Markov process on this
state space, specified by the rates of possible (independent) transitions,
which are the following (for an illustration see Fig. 1).

Figure 1:  Illustration of allowed transitions in the contact process. Allowed
transitions of an active (infected) site and two of its neighbours, connected
by two links ($12$ and $23$). Active (inactive) sites are depicted by red
(white) dots. The corresponding transition rates are written on the arrows.

First, infected sites become spontaneously healthy with a rate $\mu_{i}$,
which may be site-dependent. Second, infected sites attempt to infect their
adjacent sites with rates $\lambda_{ij}$, and the trial is successful if the
target site $j$ was healthy. Again, the infection rates $\lambda_{ij}$ can
vary from link to link. We will assume that the infection rate from site $i$
to site $j$ are the same as that from site $j$ to site $i$, i.e.
$\lambda_{ij}=\lambda_{ji}$. This variant of the model is also known as the
susceptible-infected-susceptible (SIS) model barrat . This technical
restriction, which is necessary for the applicability of the SDRG method, is
irrelevant from the point of view of universal critical properties juhasz2013
.

The contact process exhibits a continuous, non-equilibrium phase transition
between a non-fluctuating (absorbing) phase in which all sites are healthy and
a fluctuating phase, where the density of infected sites is non-zero liggett .
For regular lattices and uniform transition rates the transition falls into
the robust universality class of directed percolation md ; odor ; odor2 ; hhl
.

If the transition rates are independent, random variables, the universality
class is well characterized in one, two- and three-dimensional regular
lattices due to large-scale Monte Carlo simulations moreira ; vd ; vfm ;
vojta_3d ; vojta_rev . The observed critical behavior is in line with the
results of the strong-disorder renormalization group method hiv ; hiv2 ; im ;
im2 . According to this, the dynamical scaling of average quantities is ultra-
slow, characterized by power-laws of the logarithm of time rather than the
time itself, whereas the static critical exponents are also altered compared
to the directed percolation class hiv ; hiv2 ; 2d . Most interestingly, the
critical exponents are independent of the form of the distribution of
transition rates. Besides low dimensional regular lattices, a similar behavior
has been found in spatially embedded networks with long-range connections
munoz ; munoz2 ; jk . These networks consist of a regular, $d$-dimensional,
hypercubic lattice and a set of long links, which exist with a probability
$p_{ij}\sim\beta d(i,j)^{-s}$, where $d(i,j)$ is the Euclidean distance
between $i$ and $j$. Note, that even for uniform transition rates these
networks contain a so-called topological disorder due to the random
connectivity of sites by long links. In the case $s=2d$, the topological
dimension is finite and varies with $\beta$ bb ; coppersmith ; juhasz2012 . As
it was demonstrated in Refs. munoz ; munoz2 ; jk  for $d=1$ and $s=2$, the
contact process shows an ultra-slow scaling and the critical exponents vary
with $\beta$, while an additional disorder in the transition rates is
irrelevant. According to a general scaling argument presented in Ref. jki , a
qualitatively similar behavior is expected in higher dimensions for $s=2d$.

### Strong-disorder renormalization group

The strong-disorder renormalization group is the key method to understand the
long-time behavior of the disordered contact process (DCP). For a general
review and a detailed introduction we refer the reader to Ref. im ; im2 . The
first application of the method to the DCP was in Ref. hiv ; hiv2 . Next, we
recapitulate the essential features of the method.

The SDRG is a sequence of iterative steps operating on blocks of sites
containing the largest rate (either an infection rate or a healing rate) of
the model, see Fig. 2.

Figure 2:  Illustration of the SDRG method. Two types of reduction steps of
the SDRG procedure: formation of a cluster (a) and elimination of a site (b).
The thickness of lines indicates the magnitude of the infection rate on the
corresponding link. The link (a) and the site (b) to be decimated is shaded in
green. For details, see the text.

If the largest rate is an infection rate, $\lambda_{ij}$, the block of sites
$i$ and $j$ is merged to a cluster, characterized by an effective healing
rate:

$\ln\tilde{\mu}_{ij}=\ln\mu_{i}+\ln\mu_{j}-\ln\lambda_{ij}+\ln 2.$ (1)

This simplification is a good approximation if $\mu_{i}$ and $\mu_{j}$ are
much smaller than $\lambda_{ij}$. If the largest parameter is a healing rate,
$\mu_{i}$, site $i$ is deleted, and effective infection rates are generated
between all neighbors of site $i$:

$\ln\tilde{\lambda}_{jk}=\ln\lambda_{ji}+\ln\lambda_{ik}-\ln\mu_{i}.$ (2)

This approximation is justified under the condition
$\lambda_{ji},\lambda_{ik}\ll\mu_{i}$. Apart from a one-dimensional lattice,
it may happen that a newly formed cluster has, through its constituents, two
infection rates to a third site. Similarly, if a new infection rate is
generated, there may be a pre-existing infection rate between those two sites.
These situations are usually treated by discarding the smaller infection rate,
which is known as the maximum rule.

The above steps are applied sequentially, lowering thereby gradually the
number of degrees of freedom, as well as the rate scale
$\Omega=\max\\{\mu_{i},\lambda_{ij}\\}$. Regarding the set of clusters present
in the system, the SDRG procedure can also be viewed as a special coagulation-
annihilation process with extremal dynamics. The critical behavior of the DCP
is governed by the so called infinite-randomness fixed-point (IRFP) of the
SDRG transformation hiv ; hiv2 ; fisher ; fisher2 , at least for sufficiently
strong initial disorder hiv ; hiv2 ; vd ; nft ; nft2 ; hoyos . As the IRFP is
approached along the critical renormalization trajectory, the distribution of
logarithmic rates broadens without limits, and both types of reduction steps
become asymptotically exact. Furthermore, this feature also justifies the
applicability of the maximum rule, and even the neglection of the $\ln{2}$
term in Eq. (1). In the numerical SDRG calculations we therefore resorted to
these approximations, which, as a further advantage, enable a computationally
very efficient implementation of the SDRG procedure 3d ; 3d2 .

Following the SDRG transformation of the DCP at a given realization of random
rates down to a rate scale $\Omega\ll 1$ provides an effective DCP with a
smaller set of degrees of freedom, which approximates well the dynamics of the
original system for times $t\gg\Omega^{-1}$. Consider, for instance, the state
of the system, evolved from a fully occupied initial state, at some late time
$t$. Applying the SDRG procedure down to rate scale $\Omega=t^{-1}$, provides
the set of sites with an $O(1)$ occupation as the constituents of the clusters
of the renormalized system, while the sites eliminated during the procedure
will have an almost zero occupation probability.

In the case of a finite sample, the true steady state is always the absorbing,
empty lattice state. Accordingly, all clusters are removed during the SDRG
procedure, the last one at some rate $\Omega_{1}$, the second last at
$\Omega_{2}>\Omega_{1}$, in general, the $n$th cluster becomes irrelevant at
$\Omega_{n}>\Omega_{n-1}$. The inverse rates,
$\Omega_{1}^{-1},\Omega_{2}^{-1},\dots$, are interpreted as the mean lifetimes
of the corresponding clusters removed during the procedure. Clearly, the last
one, $\Omega_{1}^{-1}$, gives the mean lifetime of the sample needed to be
absorbed in the empty state. In finite but large systems in the vicinity of
the critical point, the last few lifetimes such as $\Omega_{1}^{-1}$ and
$\Omega_{2}^{-1}$ are typically well separated, meaning that they differ by
orders of magnitude, and this time-scale separation becomes more and more
pronounced with increasing system size. Consequently, in typical samples there
is a considerable time span, $\Omega_{1}^{-2}\ll t\ll\Omega_{1}^{-1}$, within
which practically no sites but the constituents of the lastly removed cluster
are occupied.

The structure of the lastly removed cluster can also be used to determine a
sample-dependent pseudo-critical point: making a double sized system by
glueing together two copies of the original one, the onset of the active phase
is indicated when the last cluster is different from that of the original
system 2d .

## Quasi-stationary simulation

The simulations were implemented for regular lattices in the following way. An
occupied site ($i$) is randomly chosen and, with a probability $\mu/(\mu+1)$
it is made unoccupied, or, with a probability $1/(\mu+1)$, one of its $n$
neighbors ($j$) is picked with a uniform probability and infected with the
probability $\lambda_{ij}$, provided it was healthy. The time increment
related to such a move is $\Delta t=1/N_{s}$, where $N_{s}$ is the total
number of infected sites.

For non-regular networks, where the sites may have different coordination
numbers, a different implementation is needed in order to correctly simulate
the process. Here, besides the list of infected sites, also a list of active
links is stored, which contains all directed links with an infected source
site. The number of elements of this list is denoted by $N_{e}$. Then, with a
probability $\mu N_{s}/(\mu N_{s}+N_{e})$, a healing event occurs on one of
the infected sites, chosen equiprobably, whereas, with the complementary
probability, $N_{e}/(\mu N_{s}+N_{e})$, an infection event is attempted. In
this case, an active link ($ij$) is chosen with a uniform probability from the
corresponding list, and the target site ($j$) is infected with a probability
$\lambda_{ij}$, provided it was healthy. The time increment of such an
elementary move is $\Delta t=1/(\mu N_{s}+N_{e})$.

We apply a quasi-stationary simulation with a reflecting boundary condition.
This means that, at the point where there is only a single active site in the
system, the healing event is rejected. This way we only need to validate that
the system reached the quasi-stationary state, easily done by checking whether
the mean density and its variance remain unchanged under increasing the
relaxation time. The total simulation time was typically $10^{28}$, the first
half of which is left for relaxation, while the measurements are performed in
the second half.

Figure 3:  Results in 1D systems. Asymptotic probability of being active in
the quasi-stationary simulations as the function of the spatial coordinate for
$N=1000$ sites with periodic boundary conditions and strength of disorder
$\alpha=1$. As illustrated by the purple and red clusters the SDRG is able to
efficiently capture the activity profile even with a few clusters. The
correlated clusters are disconnected fractals with a fractal dimension
$d_{f}=\frac{1+\sqrt{5}}{4}\approx 0.819$ fisher ; fisher2 , corresponding to
a highly uneven activity profile. Gaps of strongly reduced density inside the
clusters can be seen already at small scales, indicating asymptotically
disconnected clusters of activity.

First, in each sample, we determine an individual pseudo-critical point,
which, for the system sizes we use, can significantly deviate from the
ensemble average. Following the method proposed in Ref. ferreira , we identify
the pseudo-critical point with the maximum of susceptibility binder  defined
as

$\chi=N\frac{\langle\rho^{2}\rangle-\langle\rho\rangle^{2}}{\langle\rho\rangle},$
(3)

where $N$ is the total number of sites, $\rho$ denotes the global density of
active sites, and $\langle\cdot\rangle$ stands for the expected value in the
quasi-stationary state.

Having estimated the pseudo-critical point, we performed here quasi-stationary
simulations for a longer time, $10^{30}$, and measured the mean local
densities (Figs. 3-6). In addition, we started a number of independent
simulations and averaged over them, in order to avoid the (rather improbable)
possibility that the long-lasting activity is stuck at a cluster other than
the one with the longest lifetime.

For a satisfying agreement with the SDRG method for moderate system sizes, the
strength of the disorder needs to be sufficiently large, i.e. the distribution
of the logarithmic infection rates needs to be sufficiently broad. In
practice, the infection rates were chosen from a power-law distribution,
$P_{<}(\lambda)=\lambda^{1/\alpha}$, where the exponent $\alpha>0$ controls
the strength of disorder.

## Comparison of the cluster structure with the simulations

Figure 4:  Results in 2D systems. A good agreement is found between the
asymptotic density profile in the simulations and the largest clusters
obtained by the SDRG method. In 2D the fractal dimension of the clusters is
$d_{f}=1.02(2)$ 2d , illustrated here for $N=200\times 200$ sites with
periodic boundary conditions and strength of disorder $\alpha=4$. The
grayscale indicates the site probability of being active in the simulations
and the cluster size in the SDRG calculations, respectively.  Figure 5:
Results in 3D systems. We again find a good agreement between the asymptotic
density profile in the simulations and the clusters obtained by the SDRG
method. In 3D the fractal dimension is $d_{f}=1.16(2)$ 3d ; 3d2 , illustrated
here for $N=50\times 50\times 50$ sites with periodic boundary conditions and
strength of disorder $\alpha=4$. For better visibility, only sites with at
least $10\%$ of the maximum probability are indicated in the simulations,
while we show only clusters of size above the square root of the size of the
main cluster.  Figure 6:  Results in 2D with long-range interactions. Here we
show our results for algebraically decaying interaction probability with an
exponent $s=4$, for $N=100\times 100$ sites with periodic boundary conditions.
The number of short-cuts is $N/16$ and the strength of disorder is $\alpha=4$.
Red links indicate long-range connections within the main cluster. The
critical clusters are even more strongly disconnected objects, with a formally
zero fractal dimension jki .

The average order parameter of the DCP is related to the fraction of original
sites comprised by the clusters of the renormalized system, which decreases
gradually as the SDRG proceeds. Its critical scaling properties can then be
inferred from its dependence on the inverse time scale $\Omega$ and inverse
length scale (the number density of clusters) along the critical trajectory of
the SDRG im ; im2 . Instead of this, here we compare the spatial structure of
SDRG clusters with the map of local occupancies obtained by simulations in
individual samples.

Ideally, the comparison of the two methods should be done under the same
circumstances: i) for the same random sample (set of random rates), ii) at the
same time (given by $\Omega^{-1}$ in the SDRG method). The first requirement
is, however, not the appropriate choice for a fair comparison of the methods.
The reason is that the SDRG transformation, due to its approximative nature at
early stages, does not preserve the position of the critical point. For
instance, a truly critical initial system will depart from the critical
trajectory and, in order to arrive at the critical IRFP, a compensatory
initial shift in the control parameter is needed to be imposed. Therefore,
rather than the initial point of the renormalization trajectory, its end point
has to match the sample used in simulations. For this reason, the first
requirement, i.e. the complete identity of rates must be relaxed, but a notion
of “equivalence of the random environment” must still be preserved. Our
approach to overcome this controversy was the following. We chose the
infection rates randomly, while kept the healing rates constant. The latter
can then serve as a control parameter and can be used to compensate the shift
induced in the SDRG. The amount of the shift is quantitatively a priori
unknown, therefore it is natural to choose the critical point as a “common
point” to which the system can be tuned in both methods by using an indicator
of criticality internal to that method. By this construction, the transition
rates used in the two methods are not completely identical, but the difference
is only in the uniform healing rate, while the random infection rates, which
form the “random environment”, and which govern the shaping of SDRG clusters,
are identical.

Concerning the second requirement, we implemented the SDRG procedure up to the
stage at which only one cluster remains in the system. Provided the lifetimes
of clusters are well separated, this corresponds to the time span between the
lifetime of the second last and the last cluster. The state of the system in
this time span can be conveniently simulated, even without knowing the above
lifetimes, by the above described quasi-stationary simulation, which prevents
the system from getting absorbed in the empty state.

In Fig. 3, we illustrate for three random samples that the main SDRG cluster
(purple) not only contains all sites of nearly maximal probabilities, but also
captures the majority of the probability weight. For such, moderately sized,
systems smaller clusters still play a role, as shown by probability
distribution captured by the secondary cluster. Interestingly, both the main
and secondary clusters are fractured by deep gaps in the probability
distribution, as signatures of asymptotically disconnected clusters. Overall,
we see a good agreement already for moderate system sizes as illustrated in
Figs. 3-5 for $d=1,2$ and $3$ dimensional lattices as well as for a two-
dimensional lattice with long-range interactions (Fig. 6).

## Discussion

In this paper, we studied the disordered contact process (DCP) in a random
environment. By applying the combination of quasi-static simulations and an
efficient renormalization group method, we showed that the critical behavior
of the DCP is dominated by one strongly correlated cluster. Yet, as opposed to
equilibrium systems, the governing clusters are predominantly disconnected
objects, indicating that individuals who are infected at the same time
typically do not know each other directly. According to the SDRG method, such
disconnected clusters emerge from a strong, positive feedback mechanism, in
which remote sites can effectively infect each other over and over again for a
prolonged amount of time through indirect paths in the contact network phd .

Besides infection spreading, variants of the DCP have recently gained
increasing interest also in functional modeling of the brain moretti . Our
results suggest that, as opposed to traditionally expected functional ’blobs’,
strongly correlated brain regions are not necessarily directly connected. In
other words, functional connectivity might qualitatively deviate from physical
connectivity. This expectation exists in addition to the challenge that in the
brain correlations do not necessarily decay with increasing physical distance
brain_review . We leave the extension of our results to dynamical models on
experimentally obtained brain connectome datasets for future studies.

## Acknowledgements

This paper was supported by the National Research, Development and Innovation
Office - NKFIH under grants No. K128989 and No. K131458. IAK was supported by
the ’Nation’s Young Talents Scholarship’, Ministry of National Resources,
Hungary under Grant No. NTP-NFTÖ-17-C-0247.

## Author contributions statement

I.A.K initiated the project and performed the SDRG calculations. R.J.
performed the MC simulations. All authors contributed to writing the
manuscript.

## Competing interests

The authors declare no competing interests.

## Data availability

All data generated or analysed during this study are included in this
published article.

## References

  * (1) Fortuin, C.M. & Kasteleyn, P.W. On the random-cluster model : I. Introduction and relation to other models. Physica, 57, (4): 536-564 (1972).
  * (2) Carey, R. & Isaac, E.D. Magnetic Domains And Techniques For Their Observation The English University Press Ltd, London, (1966).
  * (3) Mansfield, M.L. & Douglas, J.F. Shape characteristics of equilibrium and non-equilibrium fractal clusters. AIP The Journal of Chemical Physics, 139, 044901 (2013).
  * (4) Noest, A.J. New universality for spatially disordered cellular automata and directed percolation. Phys. Rev. Lett. 57, 90 (1986).
  * (5) Noest, A.J. Power-law relaxation of spatially disordered stochastic cellular automata and directed percolation. Phys. Rev. B 38, 2715 (1988).
  * (6) Moreira, A.G. & Dickman, R. Critical dynamics of the contact process with quenched disorder. Phys. Rev. E 54, R3090 (1996).
  * (7) Vojta, T., Farquhar, A. & Mast, J. Infinite-randomness critical point in the two-dimensional disordered contact process. Phys. Rev. E 79, 011111 (2009).
  * (8) Vojta, T. & Dickison, M. Critical behavior and Griffiths effects in the disordered contact process. Phys. Rev. E 72, 036126 (2005).
  * (9) Vojta, T. Monte Carlo simulations of the clean and disordered contact process in three dimensions. Phys. Rev. E 86, 051137 (2012).
  * (10) Vojta, T. Rare region effects at classical, quantum and nonequilibrium phase transitions. J. Phys. A 39, R143 (2006).
  * (11) Hooyberghs, J., Iglói, F. & Vanderzande, C. Strong disorder fixed point in absorbing-state phase transitions. Phys. Rev. Lett. 90, 100601 (2003).
  * (12) Hooyberghs, J., Iglói, F. & Vanderzande. C. Absorbing state phase transitions with quenched disorder. Phys. Rev. E 69, 066140 (2004).
  * (13) Iglói, F. & Monthus, C. Strong disorder RG approach of random systems. Phys. Rep. 412, 277 (2005).
  * (14) Iglói, F. & Monthus, C. Strong disorder RG approach–a short review of recent developments. Eur. Phys. J. B 91, 290 (2018).
  * (15) Barrat, A., Barthélemy, M. & Vespignani, A. Dynamical Processes On Complex Networks Cambridge Univ. Press (Cambridge) (2008).
  * (16) Juhász, R. Disordered contact process with asymmetric spreading. Phys. Rev. E 87, 022133 (2013).
  * (17) Liggett, T.M. Stochastic Interacting Systems: Contact, Voter, And Exclusion Processes Springer (Berlin) (2005).
  * (18) Marro, J. & Dickman, R. Non-equilibrium Phase Transitions In Lattice Models Cambridge Univ. Press (Cambridge) (1999).
  * (19) Ódor, G. Universality In Nonequilibrium Lattice Systems World Scientific (Singapore) (2008).
  * (20) Ódor, G. Universality classes in nonequilibrium lattice systems. Rev. Mod. Phys. 76, 663 (2004).
  * (21) Henkel, M., Hinrichsen, H. & Lübeck, S. Non-equilibrium Phase Transitions Springer (Berlin) (2008).
  * (22) Kovács, I.A. & Iglói, F. Renormalization group study of the two-dimensional random transverse-field Ising model. Phys. Rev. B 82, 054437 (2010).
  * (23) Juhász, R. & Kovács, I.A. Infinite randomness critical behavior of the contact process on networks with long-range connections. J. Stat. Mech. P06003 (2013).
  * (24) Muñoz, M.A., Juhász, R., Castellano, C. & Ódor, G. Griffiths phases on complex networks. Phys. Rev. Lett. 105, 128701 (2010).
  * (25) Juhász, R., Ódor, G., Castellano, C. & Muñoz, M.A. Rare-region effects in the contact process on networks. Phys. Rev. E 85, 066125 (2012).
  * (26) Benjamini, I. & Berger, N. The diameter of long-range percolation clusters on finite cycles. Rand. Struct. Alg. 19, 102 (2001).
  * (27) Coppersmith, D., Gamarnik, D. & Sviridenko, M. The diameter of a long-range percolation graph. Rand. Struct. Alg. 21, 1 (2002).
  * (28) Juhász, R. Competition between quenched disorder and long-range connections: a numerical study of diffusion. Phys. Rev. E 85, 011118 (2012).
  * (29) Juhász, R., Kovács, I.A. & Iglói, F. Long-range epidemic spreading in a random environment. Phys. Rev. E 91, 032815 (2015).
  * (30) Fisher, D.S. Random transverse field Ising spin chains. Phys. Rev. Lett. 69, 534 (1992)
  * (31) Fisher, D.S. Critical behavior of random transverse-field Ising spin chains. Phys. Rev. B 51, 6411 (1995).
  * (32) Neugebauer, C.J., Fallert, S.V. & Taraskin, S.N. Contact process in heterogeneous and weakly disordered systems, Phys. Rev. E 74, 040101(R) (2006).
  * (33) Fallert, S.V. & Taraskin, S.N. Scaling behavior of the disordered contact process. Phys. Rev. E 79, 042105 (2009).
  * (34) Hoyos, J.A. Weakly disordered absorbing-state phase transitions. Phys. Rev. E 78, 032101 (2008).
  * (35) Kovács I.A. & Iglói, F. Infinite-disorder scaling of random quantum magnets in three and higher dimensions. Phys. Rev. B 83, 174207 (2011).
  * (36) Kovács I.A. & Iglói, F. Renormalization group study of random quantum magnets. J. Phys.: Condens. Matter 23, 404204 (2011).
  * (37) Ferreira, S.C., Castellano, C. & Pastor-Satorras, R. Epidemic thresholds of the susceptible-infected-susceptible model on networks: A comparison of numerical and theoretical results. Phys. Rev. E 86, 041125 (2012).
  * (38) Binder, K. & Heermann, D.W. Monte Carlo Simulation In Statistical Physics (Springer, Berlin) (2010).
  * (39) Kovács, I.A. Infinitely Disordered Critical Behavior In Higher Dimensional Quantum Systems (Ph.D. dissertation) Eötvös Loránd University (2013).
  * (40) Moretti, P. & Muñoz, M.A. Griffiths phases and the stretching of criticality in brain networks. Nat. Commun. 4, 2521 (2013).
  * (41) Bowman, F.D. Brain imaging analysis. Annu. Rev. Stat. Appl. 1, 61-85 (2014).


# Fractional time differential equations as a singular limit of the
Kobayashi–Warren–Carter system

Yoshikazu Giga  Graduate School of Mathematical Sciences, The University of
Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8914, Japan.
<EMAIL_ADDRESS>,  Ayato Kubo  Department of Mathematics, Faculty of
Science, Hokkaido University, Kita 10, Nishi 8, Kita-Ku, Sapporo, Hokkaido,
060-0810, Japan<EMAIL_ADDRESS>,  Hirotoshi Kuroda
Department of Mathematics, Faculty of Science, Hokkaido University, Kita 10,
Nishi 8, Kita-Ku, Sapporo, Hokkaido, 060-0810, Japan.
<EMAIL_ADDRESS>,  Jun Okamoto  Institute for the Advanced Study
of Human Biology, Kyoto University, Yoshida-Konoe-Cho, Sakyo-ku, Kyoto
606-8501, Japan<EMAIL_ADDRESS>,  Koya Sakakibara  Faculty of
Mathematics and Physics, Institute of Science and Engineering, Kanazawa
University, Kakuma-machi, Kanazawa-shi, Ishikawa 920-1192, Japan; RIKEN
iTHEMS, 2-1 Hirosawa, Wako-shi, Saitama 351-0198, Japan.
<EMAIL_ADDRESS>and  Masaaki Uesaka  Graduate School of
Mathematical Sciences, The University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo
153-8914, Japan; Arithmer Inc., ONEST Hongo Square 3F, 1-24-1 Hongo, Bunkyo-
ku, Tokyo 113-0033, Japan.

###### Abstract.

This paper is concerned with a singular limit of the Kobayashi–Warren–Carter
system, a phase field system modelling the evolutions of structures of grains.
Under a suitable scaling, the limit system is formally derived when the
interface thickness parameter tends to zero. Different from many other
problems, it turns out that the limit system is a system involving fractional
time derivatives, although the original system is a simple gradient flow. A
rigorous derivation is given when the problem is reduced to a gradient flow of
a single-well Modica–Mortola functional in a one-dimensional setting.

###### Key words and phrases:

Kobayashi–Warren–Carter system, gradient flow, singular limit, fractional time
derivative

###### 2020 Mathematics Subject Classification:

35R11, 35Q74, 35K20, 74N99

##  1\. Introduction

We consider the Kobayashi–Warren–Carter system, introduced in [11, 12, 13], to
model evolutions of structures in a multi-grain problem. It is a kind of
phase-field system in a domain $\Omega$ in $\mathbf{R}^{n}$, formally a
gradient flow of the energy

(1.1) $\displaystyle
E_{\mathrm{KWC}}^{\varepsilon}(u,v):=\int_{\Omega}\alpha_{0}(v)|\nabla
u|+E_{\mathrm{sMM}}^{\varepsilon}(v),$ (1.2) $\displaystyle
E_{\mathrm{sMM}}^{\varepsilon}(v):=\int_{\Omega}\frac{\varepsilon}{2}|\nabla
v|^{2}\,dx+\int_{\Omega}\frac{1}{2\varepsilon}F(v)\,dx.$

Here, $\alpha_{0}(v)\geq 0$ is a given function, typically
$\alpha_{0}(v)=sv^{2}$, with a constant $s>0$, and $F(v)$ is a single-well
potential, typically $F(v)=a^{2}(v-1)^{2}$, with a constant $a>0$. The
functional $E_{\mathrm{sMM}}^{\varepsilon}$ is often called a single-well
Modica–Mortola functional. The Kobayashi–Warren–Carter system is regarded as a
gradient flow with respect to $L^{2}$-inner product

$\left((u_{1},v_{1}),(u_{2},v_{2})\right)=\int_{\Omega}\alpha_{w}u_{1}u_{2}\,dx+\tau\int_{\Omega}v_{1}v_{2}\,dx,\\\
(u_{i},v_{i})\in L^{2}(\Omega)\times L^{2}(\Omega),\quad i=1,2,$

where $\alpha_{w}\geq 0$ and $\tau>0$ are weights. The function
$\alpha_{w}=\alpha_{w}(v)$ is given, but it may depend on $v\in
L^{2}(\Omega)$, so the above inner product is a Riemann metric on the tangent
bundle $TL^{2}(\Omega)$. A typical form of $\alpha_{w}(v)$ equals
$\alpha_{w}(v)=\tau_{0}v^{2}$, where $\tau_{0}$ is a positive constant. We
consider the gradient flow of $E_{\mathrm{KWC}}^{\varepsilon}$ under this
metric, and its explicit form is 2 τv_t &= εΔv - 12ε F’(v) - α’_0(v)—∇u—,
α_w(v) u_t = div ( α_0(v) ∇u—∇u— ). An explicit form in [11] corresponds to
the case when $s=\varepsilon^{2}$, $\tau_{1}=\tau\varepsilon$,
$F(v)=(v-1)^{2}$, $\alpha_{0}(v)=\varepsilon v^{2}$,
$\alpha_{w}(v)=\tau_{0}v^{2}/\varepsilon$ with $\tau_{0}>0$. In other words, 2
τ_1 v_t &= sΔv + (1-v) -2sv—∇u—,
τ_0 v^2 u_t = sdiv ( v^2 ∇u—∇u— ). The function $v$ represents an order
parameter. Where $v=1$ corresponds to a grain region, and where $v$ is away
from $1$ corresponds to grain boundaries. The function $u$ represents a
structure-like averaged angle in each grain.

We are interested in a singular limit problem for (1), (1) as
$\varepsilon\downarrow 0$. It turns out that the correct scaling of time
should be $\tau=\tau_{1}/\varepsilon$, while $\tau_{1}$ is independent of
$\varepsilon$. Since the system (1)–(1) is regarded as a gradient flow of
$E_{\mathrm{KWC}}^{\varepsilon}(u,v)$ of (1.1), we are tempting to expect that
the limit flow is the gradient flow of its limit energy $E_{\mathrm{KWC}}^{0}$
which was obtained in our papers [9, 8]. Surprisingly, this conjecture is
wrong. The limit flow contains a fractional time derivative. In this paper, we
consider the problem in a one-dimensional setting. Moreover, we consider a
special but typical case when the problem is essentially reduced to a single
equation for $v$ in (1) because handing (1) is technically involved since it
is a total variation flow type equation. This reduced problem becomes a linear
problem and is easy to discuss.

We consider (1)–(1), where $\Omega$ is an interval $\Omega=(-L,L)$ and impose
the Dirichlet boundary condition for $u$ and the Neumann boundary condition
for $v$. More precisely,

(1.3) $u(-L,t)=0,\quad u(L,t)=b>0$

while

(1.4) $v_{x}(\pm L,t)=0\quad\text{for}\quad t>0.$

We set

(1.5) $F(v)=a^{2}(v-1)^{2}\quad\text{with}\quad a\geq
0,\quad\alpha_{0}(v)=v^{2}.$

We expect that the function

$\displaystyle u^{b}(x)=\left\\{\begin{array}[]{ll}b,&x>0,\\\
0,&x<0\end{array}\right.$

with $b>0$ solves (1). Since equation (1) is of total variation flow type, the
definition of a solution is not obvious. Fortunately, under a suitable
assumption of $v$, the function $u^{b}$ solves (1) under (1.3), as shown in
the following lemma by setting $\beta=\alpha_{0}(v)$.

######  Lemma 1.1 (A stationary solution).

Assume that $\beta\in C[-L,L]$ satisfies

$\beta(0)\leq\beta(x)\quad\text{for all}\quad x\in(-L,L).$

Then $u^{b}$ solves

$\left(\beta\frac{u_{x}}{|u_{x}|}\right)_{x}=0\quad\text{in}\quad(-L,L)$

under (1.3).

We stress that the notion of a solution of the equation for $u$ in Lemma 1.1
and also (1) under the Dirichlet condition (1.3) is not obvious and will be
discussed in Section 4.

Problem (1)–(1) is reduced to

(1.7) $\frac{\tau_{1}}{\varepsilon}v_{t}=\varepsilon
v_{xx}-\frac{a^{2}(v-1)}{\varepsilon}-2b\partial_{x}(1_{x>0})v\quad\text{in}\quad\Omega\times(0,\infty)$

under the boundary condition

(1.8) $v_{x}(\pm L,t)=0,\quad t>0$

and the initial condition

(1.9) $v(x,0)=v_{0}^{\varepsilon}(x),\quad x\in\Omega,$

where $1_{x>0}$ is a characteristic function of $(0,+\infty)$, i.e., the
Heaviside function, so that its distributional derivative equals the Dirac
$\delta$ function. By a scaling transformation $y=x/\varepsilon$, (1.7)
becomes

(1.10) $\tau_{1}V_{t}=V_{yy}-a^{2}(V-1)-2b\partial_{y}(1_{y>0})V$

in $(-L/\varepsilon,L/\varepsilon)\times(0,\infty)$ for
$V=V^{\varepsilon}(y,t)=v^{\varepsilon}(\varepsilon y,t)$, where
$v^{\varepsilon}$ is a solution of (1.7). Thus, we expect this limit
$V=\lim_{\varepsilon\to 0}V^{\varepsilon}$ solves (1.10) on $\mathbf{R}$ and
is bounded. Since the solution $v^{\varepsilon}$ of (1.7) is expected to
converge to $1$ except at $x=0$, we are interested in the behaviour of
$\xi^{\varepsilon}(t)=v^{\varepsilon}(0,t)$. More precisely, we would like to
find the equation which $\xi=\lim_{\varepsilon\to 0}\xi^{\varepsilon}$ solves.
Let $E_{\mathrm{sMM}}^{0,b}(\zeta)$ be a limit energy obtained by [9] for a
set-valued function $\Xi$, defined as $\Xi(x)=\\{1\\}$ for $x\neq 0$ and
$\Xi(0)=[\zeta,1]$, $\zeta\in(0,1)$. In other words,

$\displaystyle E_{\mathrm{sMM}}^{0,b}(\zeta)$
$\displaystyle=b\zeta^{2}+2G(\zeta)=b\zeta^{2}+a(\zeta-1)^{2},\quad\text{where}\quad
G(\zeta)=\int_{1}^{\zeta}\sqrt{F(\rho)}\,d\rho.$

A key observation is to derive an equation for $\xi$. For
$\alpha\in\mathbf{R}$, $\beta>0$, we set

$f_{\beta}^{\alpha}(t):=\frac{e^{-\alpha t}t^{\beta-1}}{\Gamma(\beta)}.$

We consider well-prepared initial data in the sense that it is continuous and
solves (1.10) outside $y=0$.

######  Lemma 1.2 (Limit equation).

Assume that $\tau_{1}=1$. Let $V$ be the bounded solution of (1.10) in
$\mathbf{R}\times(0,\infty)$ with well-prepared initial data
$V(y,0)=1-ce^{-a|y|}$ with some $c\in\mathbf{R}$. Then $\xi(t)=V(0,t)$ solves

(1.11)
$\int_{0}^{t}m_{a}(t-s)\xi_{s}(s)\,ds=-\operatorname{grad}E_{\mathrm{sMM}}^{0,b}(\xi)$

with

$m_{a}(t)=2\left\\{f_{1/2}^{a^{2}}(t)+a^{2}\int_{0}^{t}f_{1/2}^{a^{2}}(s)\,ds-a\right\\}.$

Moreover, $m_{a}^{\prime}(t)<0$ and $m_{a}(t)>0$, and
$\lim_{t\to\infty}m_{a}(t)=0$.

The assumption $\tau_{1}=1$ is just for the convenience of presentation. The
formula for general $\tau_{1}$ is obtained by rescaling the time variable $t$
by $t^{\prime}=t\tau_{1}$.

Note that in case $a=0$, the left-hand side of (1.11) becomes the Caputo
derivative $\partial_{t}^{1/2}$. For even initial data, the equation (1.10) in
$\mathbf{R}\times(0,\infty)$ is reduced to the Robin boundary problem in
$(-\infty,0)\times(0,\infty)$ with

$\partial_{y}V(0,t)+bV(0,t)=0.$

The Caputo derivative appears in the equation of the boundary value.

######  Corollary 1.3.

Let $w=w(x,t)$ be the bounded solution of the heat equation

$w_{t}-w_{xx}=0\quad\text{in}\quad(-\infty,0)\times(0,\infty)$

with the Robin boundary condition

$w_{x}(0,t)+bw(0,t)=0\quad\text{for}\quad t>0.$

Assume that $w(x,0)=-c$ for some $c\in\mathbf{R}$; the boundary value
$\xi(t)=w(0,t)$ solves

$\int_{0}^{t}f_{1/2}^{0}(t-s)\xi_{s}(s)\,ds=-b\xi(t),\quad t>0.$

In other words, $\partial_{t}^{1/2}\xi=-b\xi$, where $\partial_{t}^{1/2}$ is
the Caputo half derivative.

It is well known that the fractional Laplace operator $(-\Delta)^{1/2}$ arises
as the Dirichlet–Neumann map of the Laplace equation. Here, the Caputo
derivative $\partial_{t}^{1/2}$ is obtained as the Dirichlet–Neumann map of
the heat equation. Formally, it is easy to guess since the Robin boundary
condition yields $\partial_{t}^{1/2}\xi+b\xi=0$ by replacing $\partial_{x}$
with $\partial_{t}^{1/2}$, which is natural since
$\partial_{t}=\partial_{x}^{2}$ for $w$. In a seminal paper, Caffarelli and
Silvestre [3] show that $(-\Delta)^{\gamma}$ ($0<\gamma<1$) is obtained as the
Dirichlet–Neumann map for the degenerate Laplace equation. We remark that
$\partial_{t}^{\gamma}$ is obtained as the Dirichlet–Neumann map of a
degenerate heat equation $w_{t}-(-x)^{\alpha}w_{xx}=0$ with
$\gamma=1/(2-\alpha)$ as in [3]; see Remark 2.4 at the end of Section 2. We do
not pursue this problem in this paper.

Since the equation is linear and of constant coefficients, we can use the
Laplace transform to obtain the desired equation. As we expect, proving the
convergence of $V^{\varepsilon}$ and $v^{\varepsilon}$ is not tricky.

######  Lemma 1.4 (Convergence).

Let $V$ be the bounded solution of (1.10) in $\mathbf{R}\times(0,\infty)$ with
initial data $V_{0}\in C(\mathbf{R})\cap L^{\infty}(\mathbf{R})$. Let
$v^{\varepsilon}$ be the solution of (1.7) under (1.8) and (1.9). Assume that
$V^{\varepsilon}(x,0)\in C\left[-L/\varepsilon,L/\varepsilon\right]$ converges
to $V_{0}$ in the sense that

$\lim_{\varepsilon\downarrow 0}\sup_{|x|\leq
L/\varepsilon}\left|V^{\varepsilon}(x,0)-V_{0}(x)\right|=0.$

Then $V^{\varepsilon}$ converges to $V$ locally uniformly in
$\mathbf{R}\times[0,\infty)$. In particular, $v^{\varepsilon}$ converges to
$1$ locally uniformly in
$\left(\mathbf{R}\backslash\\{0\\}\right)\times[0,\infty)$, and
$\xi^{\varepsilon}$ converges to $\xi(t)=V(0,t)$ locally uniformly in
$[0,\infty)$.

Applying Lemma 1.2 and Lemma 1.4, we can obtain a characterisation of the
limit equation.

######  Theorem 1.5.

Assume that $\tau_{1}=1$. Let $v^{\varepsilon}$ be the solution of (1.7) under
(1.8) and (1.9). Assume that $v_{0}^{\varepsilon}$ is well-prepared in the
sense

$\sup_{x\in\Omega}\left|v_{0}^{\varepsilon}(x)-\left(1-ce^{-a|x|/\varepsilon}\right)\right|\rightarrow
0$

as $\varepsilon\to 0$ with some $c$ independent of $\varepsilon$. Then
$\xi^{\varepsilon}(t)$ converges to $\xi$ locally uniformly in $[0,\infty)$,
and $\xi$ solves (1.11). Moreover, the graph of $v_{\varepsilon}$ converges to
a set-valued function $\Xi$ of the form

$\Xi(x,t)=\left\\{\begin{array}[]{l}\\{1\\}\quad\text{for}\quad x=0,\\\
\text{closed interval between}\ \xi(t)\ \text{and}\ 1\ \text{for}\ x\neq
0.\end{array}\right.$

The convergence is in the sense of the Hausdorff distance of graphs over
$[-L,L]\times[0,T]$ for any $T>0$.

We also handle initial data not necessarily well-prepared. In this case, the
equation (1.11) is altered because there are a few lower-order terms. If
$V(y,0)=1-ce^{-\mu|x|}$ with $\mu\neq a$, we still get an explicit form
corresponding to (1.11). In these cases, the solution $\xi$ is explicitly
represented by using the error function. We shall discuss these extensions for
non-well-prepared data and the proof of Lemma 1.2 in Section 2. We also
calculate numerically how the solution of the $\varepsilon$-problem converges
by comparing it with the explicit solution of (1.11) and more general
equations (2.16).

We now come back to the singular limit problem of the Kobayashi–Warren–Carter
system (1)–(1) with $\tau=\tau_{1}/\varepsilon$, $F(v)=a^{2}(v-1)^{2}$,
$\alpha_{0}(v)=v^{2}$ as $\varepsilon\to 0$. The equation (1) is of the total
variation flow type, and its well-posedness for given $v$ is known when $v$ is
independent of time and $\inf\alpha_{w}(v)>0$; in this case, (1) is the
gradient flow of the weighted total variation $\int\alpha_{0}(v)|\nabla u|$
under $\alpha_{w}(v)$-weighted $L^{2}$ inner product if we impose the natural
boundary condition like Neumann boundary condition. As in [4], we consider (1)
for piecewise constant functions. For an interval $I=(p,q)$ and its given
division $p=p_{0}<p_{1}<\cdots<p_{m}=q$, we consider a piecewise constant
function of form

$u(x,t)=h_{j}(t),\quad p_{j-1}<x\leq p_{j}\quad(j=1,\ldots,m).$

We interpret a solution by mimicking the notion of a solution when $v$ is
independent of time and $\inf\alpha_{w}(v)>0$ (under the periodic boundary
condition for simplicity). It is of the form

(1.12)
$\left\\{\begin{array}[]{l}\alpha_{w}(v)\partial_{t}u=\operatorname{div}(v^{2}z),\quad|z|\leq
1,\\\ z(p_{j})=\chi_{j},\quad j=0,1,\ldots,m-1,\end{array}\right.$

where $\chi_{i}=\operatorname{sgn}(h_{j+1}-h_{j})$; we identify $p_{m}=p_{0}$
by periodicity. See Figure 1.

Figure 1. graph of $u$
$\chi_{1}=1$, $\chi_{2}=-1$

The function $z$ is called a Cahn-Hoffman vector field.

For a given $v$, it is unclear whether the solution stays in a class of
piecewise constant functions [4], known as a facet-splitting problem. If $v$
is constant, it is well known that the solution stays in a class of piecewise
constant functions.

We postulate that our system (1.12) has a (spatially) piecewise constant
solution. Integrating the first equation of (1.12) in $(p_{j-1},p_{j})$ yields

$\int_{p_{j-1}}^{p_{j}}\alpha_{w}\left(v(x,t)\right)\,dx\frac{d}{dt}h_{j}(t)=v^{2}(p_{j},t)\chi_{j}-v^{2}(p_{j-1},t)\chi_{j-1}.$

If $v_{\varepsilon}\to\Xi$ in the graph topology as $\varepsilon\to 0$, we
arrive at

(1.13)
$\alpha_{w}(1)\frac{d}{dt}h_{j}=\frac{\xi_{j}^{2}\chi_{j}-\xi_{j-1}^{2}\chi_{j-1}}{p_{j}-p_{j-1}}.$

Since jump $b_{j}=|h_{j+1}-h_{j}|$ determines the $\xi_{j}$-equation as
indicated in Theorem 1.5, the equation for $\xi_{j}$ (assuming $\tau_{1}=1$)
is expected to be

(1.14)
$M_{a}\partial_{t}\xi_{j}=-\operatorname{grad}E_{\mathrm{sMM}}^{0,b_{j}}(\xi_{j})\left(=-2\left(\left(|h_{j+1}-h_{j}|+a\right)\xi_{j}-a\right)\right),$

where $M_{a}f=\int_{0}^{t}m_{a}(t-s)f(s)\,ds$ for well-prepared initial data.
Thus, the singular limit equation of the Kobayashi–Warren–Carter equation
(1)–(1) as $\varepsilon\to 0$ (under the periodic boundary condition) is
expected to be the system (1.13)–(1.14) for $h_{j}$ and $\xi_{j}$. If we
consider other boundary conditions, we always impose the homogeneous Neumann
boundary condition for $v$, like (1.4). Considering the Dirichlet boundary
condition for $u$, $h_{1}$ and $h_{m}$ are fixed time-independent constants.
We can impose the Neumann condition for $u$; in this case, imposing
$z(p_{0})=z(p_{m})=0$ in (1.12) is natural. It is rather standard [21] to
construct a unique local-in-time solution for a system of a fractional
differential equation and an ordinary differential equation.

We note that the solvability of the initial value problem for the original
Kobayashi–Warren–Carter system (1)–(1) is still an open problem, even in a
one-dimensional setting. If we write it in the form of (1)–(1), the difficulty
stems from the fact that the weights $\alpha_{w}$ and $\alpha_{0}$ can vanish
somewhere. In the literature, $\alpha_{w}$ is assumed to be away from zero. If
$\alpha_{0}$ is allowed to vanish, $\Delta u$ term is added in the right-hand
side of (1). For example, instead of considering (1), we consider

(1.15)
$\tau_{0}(v^{2}+\delta)u_{t}=s\operatorname{div}\left((v^{2}+\delta^{\prime})\nabla
u/|\nabla u|+\mu\nabla u\right)$

with $\delta>0$, $\delta^{\prime}\geq 0$ and $\mu\geq 0$ such that
$\delta^{\prime}+\mu>0$. The existence of a solution to (1) and (1.15) with
its large time behaviour is established in [10, 17, 18, 24, 25, 27] under
several homogeneous boundary conditions. Unfortunately, the uniqueness of
their solution is only known in a one-dimensional setting under the relaxation
term $\mu>0$ [10, Theorem 2.2]. The extension of these results to
inhomogeneous boundary condition is not difficult. In [19], under non-
homogeneous Dirichlet boundary conditions, structured patterns of stationary
(i.e., time-independent) solutions were studied. In a one-dimensional setting,
they thoroughly characterised all stationary solutions. In this paper, we do
not discuss convergence problems as $\varepsilon\to 0$.

This paper is organised as follows. In Section 2, we derive the limit equation
both for well-prepared and non-well-prepared initial data. Most of the
calculations there are very explicit. In Section 3, we prove Lemma 1.4.
Section 4 gives a rigorous definition of the Dirichlet problem for (1),
assuming that $v$ is given and $\alpha_{w}(v)\equiv 1$. In Section 5, we give
several numerical tests.

##  2\. Derivation of equations with fractional time derivative

We begin with recalling several elemental properties of the Laplace transform

$\mathcal{L}[g](\lambda):=\int_{0}^{\infty}e^{-\lambda
s}g(s)\,ds,\quad\lambda>0$

for a locally integrable function $g$ in $[0,\infty)$. By definition,

(2.1) $\mathcal{L}\left[e^{-\mu
s}g(s)\right](\lambda)=\mathcal{L}[g](\lambda+\mu),$

and by definition of the Gamma function, we see

$\mathcal{L}[s^{\beta-1}](\lambda)=\Gamma(\beta)\lambda^{-\beta}.$

We now arrive at a well-known formula

(2.2) $\mathcal{L}[f_{\beta}^{\alpha}](\lambda)=\mathcal{L}[e^{-\alpha
s}f_{\beta}^{0}]=(\lambda+\alpha)^{-\beta}$

for $\alpha,\beta>0$.

######  Lemma 2.1.

For $a>0$, let $m_{a}(t)$ be

$m_{a}(t)=2\left\\{f_{1/2}^{a^{2}}(t)+a^{2}\int_{0}^{t}f_{1/2}^{a^{2}}(s)\,ds-a\right\\}.$

Then

$\mathcal{L}[m_{a}](\lambda)=2\left(\frac{\sqrt{\lambda+a^{2}}-a}{\lambda}\right).$

Moreover, $m_{a}^{\prime}(t)<0$ and $m_{a}(t)>0$ for $t>0$ with
$\lim_{t\to\infty}m_{a}(t)=0$.

###### Proof.

We recall that the Laplace transform of the convolution
$(g_{1}*g_{2})(t)=\int_{0}^{t}g_{1}(t-s)g_{2}(s)\,ds$ is given by

(2.3) $\mathcal{L}[g_{1}*g_{2}]=\mathcal{L}[g_{1}]\mathcal{L}[g_{2}].$

If we take $g_{1}\equiv 1$, we see

(2.4)
$\mathcal{L}\left[\int_{0}^{t}g_{2}(s)\,ds\right](\lambda)=\mathcal{L}[g_{2}](\lambda)\lambda^{-1}$

since $\mathcal{L}[g_{1}](\lambda)=\lambda^{-1}$ by (2.2). Setting
$\alpha=a^{2}$, we apply (2.2) and (2.4) to get

$\displaystyle\mathcal{L}\left[\frac{m_{a}}{2}\right]$
$\displaystyle=(\lambda+\alpha)^{-1/2}+\alpha(\lambda+\alpha)^{-1/2}\lambda^{-1}-a\lambda^{-1}$
$\displaystyle=\left((\lambda+\alpha)^{1/2}-a\right)\lambda^{-1},$

which yields the desired formula for $\mathcal{L}[m_{a}]$.

We differentiate $m_{a}/2$ to get

$\displaystyle\frac{1}{2}m_{a}^{\prime}(t)$
$\displaystyle=(-\alpha)f_{1/2}^{\alpha}(t)-\frac{1}{2t}f_{1/2}^{\alpha}(t)+\alpha
f_{1/2}^{\alpha}(t)$
$\displaystyle=-\frac{1}{2t}f_{1/2}^{\alpha}(t)<0\quad\text{for}\quad t>0.$

We observe that

$\lim_{t\to\infty}\int_{0}^{t}f_{1/2}^{\alpha}(s)\,ds=\frac{1}{\Gamma\left(1/2\right)}\int_{0}^{\infty}e^{-\alpha
s}s^{\frac{1}{2}-1}\,ds=\alpha^{-1/2}.$

Thus, $\lim_{t\to\infty}m_{a}(t)=0$ since
$\lim_{t\to\infty}f_{1/2}^{\alpha}(t)=0$, which this implies that $m_{a}(t)>0$
for $t>0$ since $m_{a}^{\prime}(t)<0$. ∎

###### Proof of Lemma 1.2.

Since the property of $m_{a}$ has been proved in Lemma 2.1, it suffices to
derive (1.11). Studying the equation for $w=V-1$ instead of (1.10) is more
convenient. The equation (1.10) for $w$ (with $\tau_{1}=1$) becomes

(2.5) $w_{t}-w_{xx}+a^{2}w+2b(w+1)\delta=0,\quad x\in\mathbf{R},\quad t>0.$

The initial data $w(x,0)$ equals

$w(x,0)=w_{0}(x)$

with $w_{0}(x)=-ce^{-a|x|}$. Let $\widehat{w}(x,\lambda)$ be the Laplace
transform of $w$ in the $t$ variable, i.e.,

$\widehat{w}(x,\lambda)=\mathcal{L}[w](x,\lambda)=\int_{0}^{\infty}e^{-\lambda
t}w(x,t)\,dt.$

We note that

$\widehat{w_{t}}=\lambda\widehat{w}-w_{0},\quad\left.w\right|_{t=0}=w_{0}.$

Taking the Laplace transform of (2.5), we arrive at

(2.6)
$\lambda\widehat{w}-w_{0}-\widehat{w}_{xx}+a^{2}\widehat{w}+2b\left(\widehat{w}+\lambda^{-1}\right)\delta=0,$

a second-order linear ordinary differential equation with a jump of the
derivatives $\widehat{w}_{x}$. Since the coefficients are constants, we can
solve (2.6) explicitly.

A bounded solution satisfying

$\lambda\widehat{w}-w_{0}-\widehat{w}_{xx}+a^{2}\widehat{w}=0\quad\text{for}\quad
x>0\quad\text{and}\quad x<0$

is of form

$\widehat{w}=Ae^{-\sqrt{\lambda+a^{2}}|x|}-\frac{c}{\lambda}e^{-a|x|}$

with $A\in\mathbf{R}$. One should determine $A$ so that (2.6) holds. We
observe

$\widehat{w}_{x}=\left(-A\sqrt{\lambda+a^{2}}e^{-\sqrt{\lambda+a^{2}}|x|}+\frac{ca}{\lambda}e^{-a|x|}\right)(\operatorname{sgn}x)$

so that

$\widehat{w}_{xx}=\left(-2A\sqrt{\lambda+a^{2}}+\frac{2ca}{\lambda}\right)\delta+A(\lambda+a^{2})e^{-\sqrt{\lambda+a^{2}}|x|}-\frac{ca^{2}}{\lambda}e^{-a|x|}$

since $(\operatorname{sgn}x)^{\prime}=2\delta$. The equation (2.6) imposes
that

$\text{the}\ \delta\ \text{part of}\
\widehat{w}_{xx}=2b(\widehat{w}+\lambda^{-1})\delta,$

which implies

$-2A\sqrt{\lambda+a^{2}}+\frac{2ca}{\lambda}=2b\left(\widehat{w}(0,\lambda)+\lambda^{-1}\right).$

Since $\widehat{w}(0,\lambda)=A-c\lambda^{-1}$, we end up with

$\left(\sqrt{\lambda+a^{2}}+b\right)A=-\frac{b(1-c)}{\lambda}+\frac{ca}{\lambda}.$

We set $\eta(t)=w(0,t)$ and observe that

$\widehat{\eta}=A-\frac{c}{\lambda}=\frac{1}{\sqrt{\lambda+a^{2}}+b}\frac{1}{\lambda}(ca+bc-b)-\frac{c}{\lambda}.$

We next extract the time derivative of $\eta$. We note that

$\widehat{\eta_{t}}=\lambda\widehat{\eta}-\eta(0)$

so that $\widehat{\eta_{t}}=\lambda\widehat{\eta}+c$ since $\eta(0)=-c$. Then

$\widehat{\eta_{t}}=\frac{1}{\sqrt{\lambda+a^{2}}+b}(ca+bc-b).$

Multiplying $\sqrt{\lambda+a^{2}}+b$ and subtracting $a\widehat{\eta_{t}}$
from both sides, we get

$\displaystyle\left(\sqrt{\lambda+a^{2}}-a\right)\widehat{\eta_{t}}$
$\displaystyle=-(b+a)\widehat{\eta_{t}}+(ca+bc-b)$
$\displaystyle=-(b+a)\lambda\widehat{\eta}-(b+a)c+(ca+bc-b)$
$\displaystyle=\lambda\left(-(b+a)\widehat{\eta}-\frac{b}{\lambda}\right).$

The energy for $\xi$ equals

$E_{\mathrm{sMM}}^{0,b}(\xi)=b\xi^{2}+a(\xi-1)^{2}.$

Since $\eta=\xi-1$, we see
$\operatorname{grad}E_{\mathrm{sMM}}^{0,b}(\xi)=2(b+a)\eta+2b$. Thus

(2.7)
$2\left(\frac{\sqrt{\lambda+a^{2}}-a}{\lambda}\right)\widehat{\eta_{t}}=-2(b+a)\widehat{\eta}-2b\widehat{1}.$

By Lemma 2.1 and (2.3), we now conclude that $\xi=\eta+1$ solves (1.11) by
taking the inverse Laplace transform of (2.7). ∎

In Lemma 1.2, we only considered well-prepared initial data, meaning that
$V(y,t)$ is a stationary solution of (1.10) in $\mathbf{R}\backslash\\{0\\}$.
We take this opportunity to write a general equation corresponding to (1.11)
starting from general initial data. We set

$(M_{a}f)(t)=m_{a}*f(t)=\int_{0}^{t}m_{a}(t-s)f(s)\,ds,$

and the equation (1.11) is

(2.8) $M_{a}\xi_{s}=-\operatorname{grad}E_{\mathrm{sMM}}^{0,b}(\xi).$

For general initial data, our equation corresponding to (2.8) becomes more
complicated than (2.8).

We set

$G_{\lambda}^{a}(y)=\frac{e^{-\sqrt{\lambda+a^{2}}|y|}}{2\sqrt{\lambda+a^{2}}},$

which is the Green function of $-\partial_{y}^{2}+\sigma$ with
$\sigma=\lambda+a^{2}$, i.e.,

$-\partial_{y}^{2}G_{\lambda}^{a}+(\lambda+a^{2})G_{\lambda}^{a}=\delta.$

For $w_{0}\in L^{\infty}(\mathbf{R})$, we set

$g^{a}(\lambda,w_{0})=(G_{\lambda}^{a}*_{x}w_{0})(0),$

where $*_{x}$ denotes the convolution in space, i.e.,

$(G_{\lambda}^{a}*_{x}w_{0})(y)=\int_{-\infty}^{\infty}G_{\lambda}^{a}(y-z)w_{0}(z)\,dz.$

######  Lemma 2.2.

Let $\tau_{1}=1$. Let $V$ be the bounded solution of (1.10) in
$\mathbf{R}\times(0,\infty)$ with initial data $V(y,0)=1+w_{0}$, where $w_{0}$
is bounded and Lipschitz continuous in $\mathbf{R}$. Then $\xi(t)=V(0,t)$
solves

(2.9)
$M_{a}\xi_{t}+m_{a}\left(\xi(0)-1\right)-\mathcal{L}^{-1}\left[2\sqrt{\lambda+a^{2}}g^{a}(\lambda,w_{0})\right]=-\operatorname{grad}E_{\mathrm{sMM}}^{0,b}(\xi).$

###### Proof.

If $w_{0}$ is Lipschitz, then we easily see that $|w_{t}|t^{1/2}$ is bounded
(near $t=0$) so that the integrand of

$\mathcal{L}[w_{t}](\lambda)=\int_{0}^{\infty}e^{-\lambda s}w_{s}(x,s)\,ds$

is integrable for all $x\in\mathbf{R}$ near $s=0$; see Lemma 3.1.

We argue in the same way to prove Lemma 1.2. We begin with (2.6), where
$w=V-1$. Its solution outside $x=0$ is given as

$\widehat{w}=Ae^{-\sqrt{\lambda+a^{2}}|x|}+g^{a}(\lambda,w_{0})$

with $A\in\mathbf{R}$. As before, we determine $A$ and obtain that
$\eta=\xi-1$ satisfies

(2.10)
$\widehat{\eta}(\lambda)=\frac{-b}{\sqrt{\lambda+a^{2}}+b}\left(g^{a}(\lambda,w_{0})+\frac{1}{\lambda}\right)+g^{a}(\lambda,w_{0}).$

Since $\widehat{\eta_{t}}=\lambda\widehat{\eta}-\eta(0)$, we proceed with

$\displaystyle\left(\sqrt{\lambda+a^{2}}-a\right)\widehat{\eta_{t}}$
$\displaystyle=-(b+a)\widehat{\eta_{t}}-\lambda
b\left(g^{a}(\lambda,w_{0})+\frac{1}{\lambda}\right)+\left(\sqrt{\lambda+a^{2}}+b\right)\left(\lambda
g^{a}(\lambda,w_{0})-\eta(0)\right)$
$\displaystyle=-(b+a)\lambda\widehat{\eta}+(b+a)\eta(0)-\lambda
b\left(g^{a}(\lambda,w_{0})+\frac{1}{\lambda}\right)$ $\displaystyle\hskip
150.0pt+\left(\sqrt{\lambda+a^{2}}+b\right)\left(\lambda
g^{a}(\lambda,w_{0})-\eta(0)\right)$
$\displaystyle=-(b+a)\lambda\widehat{\eta}-\lambda\frac{b}{\lambda}+\left(a-\sqrt{\lambda+a^{2}}\right)\eta(0)+\sqrt{\lambda+a^{2}}\lambda
g^{a}(\lambda,w_{0}).$

Thus

$\frac{2\left(\sqrt{\lambda+a^{2}}-a\right)}{\lambda}\widehat{\eta_{t}}=-2(b+a)\widehat{\eta}\\\
-2\frac{b}{\lambda}-\frac{2\left(\sqrt{\lambda+a^{2}}-a\right)}{\lambda}\eta(0)+2\sqrt{\lambda+a^{2}}g^{a}(\lambda,w_{0}).$

Taking the inverse Laplace transform, we obtain

$M_{a}\eta_{t}=-\operatorname{grad}E_{\mathrm{sMM}}^{0,b}(\eta+1)-m_{a}(t)\eta(0)+\mathcal{L}^{-1}\left(2\sqrt{\lambda+a^{2}}g^{a}(\lambda,w_{0})\right),$

the same as (2.9). ∎

The term $\mathcal{L}^{-1}\left[2\sqrt{\lambda+a^{2}}g^{a}\right]$ has a more
explicit form. Let $E(x,t)$ be the Gauss kernel, i.e.,

$E(x,t)=\frac{1}{\sqrt{4\pi t}}e^{-|x|^{2}/4t}.$

We know its Laplace transform (as a function of $t$) is

$\mathcal{L}[E](\lambda)=\frac{1}{2\lambda^{1/2}}e^{-|x|\lambda^{1/2}}.$

Suppose we set $E^{a}=e^{-a^{2}t}E$, then

$\mathcal{L}[E^{a}](\lambda)=\frac{1}{2\sqrt{\lambda+a^{2}}}e^{-\sqrt{\lambda+a^{2}}|x|}=G_{\lambda}^{a}(x).$

Since $\mathcal{L}\left[f_{1/2}^{a^{2}}\right]=(\lambda+a^{2})^{-1/2}$ and
$\mathcal{L}[g_{1}*g_{2}]=\mathcal{L}[g_{1}]\mathcal{L}[g_{2}]$, we end up
with

$\mathcal{L}\left[f_{1/2}^{a^{2}}*E^{a}\right](\lambda)=\frac{1}{2(\lambda+a^{2})}e^{-\sqrt{\lambda+a^{2}}|x|}.$

Thus

$\mathcal{L}\left[\partial_{t}\left(f_{1/2}^{a^{2}}*E^{a}\right)+a^{2}\left(f_{1/2}^{a^{2}}*E^{a}\right)\right]=\frac{1}{2}e^{-\sqrt{\lambda+a^{2}}|x|}.$

Since
$\sqrt{\lambda+a^{2}}g^{a}=\frac{1}{2}\left(e^{-\sqrt{\lambda+a^{2}}|x|}*_{x}w_{0}\right)(0,t)$,
we conclude that

$\mathcal{L}^{-1}\left[2\sqrt{\lambda+a^{2}}g^{a}\right]=2\partial_{t}\left(f_{1/2}^{a^{2}}*(E^{a}*_{x}w_{0})(0)\right)+2f_{1/2}^{a^{2}}*(E^{a}*_{x}w_{0})(0).$

If $w_{0}=-ce^{-\mu|x|}$, we observe that

(2.11)
$\mathcal{L}^{-1}\left[2\sqrt{\lambda+a^{2}}g^{a}\right]=-2c\left(f_{1/2}^{a^{2}}-\mu
e^{(\mu^{2}-a^{2})t}\operatorname{erfc}\left(\mu\sqrt{t}\right)\right),$

where $\operatorname{erfc}$ denotes the error function, i.e.,

$\operatorname{erfc}(s)=\frac{2}{\sqrt{\pi}}\int_{s}^{\infty}e^{-\tau^{2}}\,d\tau.$

Indeed, we proceed with

(2.12) $\displaystyle\begin{aligned}
g^{a}(\lambda)&=G_{\lambda}^{a}*_{x}\left(-ce^{-\mu|x|}\right)(0)=-c\frac{1}{2\sqrt{\lambda+a^{2}}}\int_{-\infty}^{\infty}e^{-\sqrt{\lambda+a^{2}}|0-y|}e^{-\mu|y|}\,dy\\\
&=-\frac{c}{\sqrt{\lambda+a^{2}}}\int_{0}^{\infty}e^{-\left(\mu+\sqrt{\lambda+a^{2}}\right)y}\,dy=-\frac{c}{\sqrt{\lambda+a^{2}}\left(\mu+\sqrt{\lambda+a^{2}}\right)}.\end{aligned}$

By a direct calculation, we observe that

$\displaystyle\mathcal{L}\left[\operatorname{erfc}\left(\mu\sqrt{t}\right)\right]=\int_{0}^{\infty}e^{-\lambda
t}\left(\frac{2}{\sqrt{\pi}}\int_{\mu\sqrt{t}}^{\infty}e^{-s^{2}}\,ds\right)dt$
$\displaystyle\quad=\frac{2}{\sqrt{\pi}}\int_{0}^{\infty}\left[\int_{0}^{s^{2}/\mu^{2}}e^{-s^{2}}e^{-\lambda
t}\,dt\right]ds=\frac{2}{\sqrt{\pi}}\int_{0}^{\infty}e^{-s^{2}}\left[\int_{0}^{s^{2}/\mu^{2}}e^{-\lambda
t}\,dt\right]ds$
$\displaystyle\quad=\frac{2}{\sqrt{\pi}}\int_{0}^{\infty}e^{-s^{2}}\left(\frac{1-e^{\lambda
s^{2}/\mu^{2}}}{\lambda}\right)ds=\frac{2}{\sqrt{\pi}}\frac{1}{\lambda}\left(\frac{\sqrt{\pi}}{2}-\frac{\sqrt{\pi}}{2}\sqrt{\frac{\mu^{2}}{\mu^{2}+\lambda}}\right)$
$\displaystyle\quad=\frac{1}{\lambda}\left(1-\frac{\mu}{\sqrt{\lambda+\mu^{2}}}\right).$

Thus,

(2.13) $\displaystyle\begin{aligned} \mathcal{L}\left[f_{1/2}^{0}-\mu
e^{\mu^{2}t}\operatorname{erfc}\left(\mu\sqrt{t}\right)\right]&=\frac{1}{\sqrt{\lambda}}-\frac{\mu}{\lambda-\mu^{2}}\left(1-\frac{\mu}{\sqrt{\lambda}}\right)=\frac{\lambda-\mu^{2}-\sqrt{\lambda}\mu+\mu^{2}}{\sqrt{\lambda}(\lambda-\mu^{2})}\\\
&=\frac{\sqrt{\lambda}-\mu}{\lambda-\mu^{2}}=\frac{1}{\sqrt{\lambda}+\mu},\end{aligned}$

where we invoked (2.2) and (2.1). We set $q^{\mu}(t)=f_{1/2}^{0}(t)-\mu
e^{\mu^{2}t}\operatorname{erfc}\left(\mu\sqrt{t}\right)$ so that (2.13)
becomes

(2.14) $\mathcal{L}[q^{\mu}]=\frac{1}{\sqrt{\lambda}+\mu}.$

Using (2.1) again, we now obtain (2.11).

The formula (2.13) gives another representation of $m_{a}$. Indeed, since

$\frac{1}{2}\mathcal{L}[m_{a}](\lambda)=\frac{\sqrt{\lambda+a^{2}}-a}{\lambda}=\frac{1}{\sqrt{\lambda+a^{2}}+a},$

we see

(2.15)
$\frac{1}{2}m_{a}(t)=e^{-a^{2}t}q^{a}(t)=f_{1/2}^{a^{2}}(t)-a\operatorname{erfc}\left(a\sqrt{t}\right),$

in particular, implies that $m_{a}(t)\leq 2f_{1/2}^{a^{2}}(t)$ and

$\operatorname{erfc}\left(a\sqrt{t}\right)=1-a\int_{0}^{t}f_{1/2}^{a^{2}}(s)\,ds.$

Moreover, by (2.15), we see

$q^{a}(t)=e^{a^{2}t}m_{a}(t)/2.$

Therefore, the positivity of $m_{a}$ (Lemma 2.1) implies $q^{a}(t)>0$ for
$t>0$. The formula (2.9) in Lemma 2.2 becomes

(2.16)
$M_{a}\xi_{t}+2cae^{-a^{2}t}\operatorname{erfc}\left(a\sqrt{t}\right)-2c\mu
e^{(\mu^{2}-a^{2})t}\left(\operatorname{erfc}\left(\mu\sqrt{t}\right)\right)=-\operatorname{grad}E_{\mathrm{sMM}}^{0,b}(\xi)$

if $w_{0}=-ce^{-\mu|x|}$. If $\mu=a$, this is reduced to (1.11) (or (2.8)).

One can give an explicit form of a solution of (2.16) starting from
$\xi(0)=1-c$. We plug (2.12) into (2.10) to get

(2.17) $\displaystyle\begin{aligned}
\widehat{\eta}(\lambda)&=\frac{\sqrt{\lambda+a^{2}}}{\sqrt{\lambda+a^{2}}+b}g^{a}(\lambda)+\frac{-b}{\lambda\left(\sqrt{\lambda+a^{2}}+b\right)}\\\
&=\frac{1}{\sqrt{\lambda+a^{2}}+b}\left(\frac{-c}{\mu+\sqrt{\lambda+a^{2}}}-\frac{b}{\lambda}\right).\end{aligned}$

Since

$\mathcal{L}\left[e^{-a^{2}t}q^{\mu}\right]=\frac{1}{\sqrt{\lambda+a^{2}}+\mu}$

by (2.13) and (2.1), we have

$\displaystyle\eta$
$\displaystyle=-ce^{-a^{2}t}q^{b}*e^{-a^{2}t}q^{\mu}-b\int_{0}^{t}e^{-a^{2}s}q^{b}(s)\,ds$
$\displaystyle=-ce^{-a^{2}t}(q^{b}*q^{\mu})-b\int_{0}^{t}e^{-a^{2}s}q^{b}(s)\,ds.$

However, the calculation of $q^{b}*q^{\mu}$ is quite involved, and it is
easier to calculate $\widehat{\eta}$ in (2.17) more. We proceed

$\displaystyle\widehat{\eta}(\lambda)$
$\displaystyle=\frac{1}{\sqrt{\lambda+a^{2}}+b}\left(\frac{c\left(\sqrt{\lambda+a^{2}}-\mu\right)}{\lambda+a^{2}-\mu^{2}}-\frac{b}{\lambda}\right)$
$\displaystyle=\frac{1}{\sqrt{\lambda+a^{2}}+b}\left(\frac{-c\left(\left(\sqrt{\lambda+a^{2}}+b\right)-b-\mu\right)}{\lambda+(a^{2}-\mu^{2})}-\frac{b}{\lambda}\right)$
$\displaystyle=\frac{1}{\sqrt{\lambda+a^{2}}+b}\frac{(b+\mu)c}{\lambda+(a^{2}-\mu^{2})}-\frac{c}{\lambda+(a^{2}-\mu^{2})}-\frac{1}{\sqrt{\lambda+a^{2}}+b}\frac{b}{\lambda}$
$\displaystyle\equiv I+II+III.$

It is easy to see that

$\mathcal{L}^{-1}(II)=-ce^{-(a^{2}-\mu^{2})t}.$

As we already observed,

$\mathcal{L}^{-1}(III)=-b\int_{0}^{t}e^{-a^{2}s}q^{b}(s)\,ds.$

For the first term,

$\mathcal{L}^{-1}(I)=-c(b+\mu)e^{-(a^{2}-\mu^{2})t}*e^{-a^{2}t}q^{b}.$

By definition,

$\displaystyle e^{-(a^{2}-\mu^{2})t}*e^{-c^{2}t}q^{b}$
$\displaystyle=\int_{0}^{t}e^{-(a^{2}-\mu^{2})(t-s)}e^{-a^{2}s}q^{b}(s)\,ds$
$\displaystyle=e^{-(a^{2}-\mu^{2})t}\int_{0}^{t}e^{-\mu^{2}s}q^{b}(s)\,ds.$

We thus conclude that

(2.18)
$\eta(t)=c\left[(\mu+b)\int_{0}^{t}e^{-\mu^{2}s}q^{b}(s)\,ds-1\right]e^{-(a^{2}-\mu^{2})t}-b\int_{0}^{t}e^{-a^{2}s}q^{b}(s)\,ds.$

Thus, $\xi=\eta+1$ is the solution of (2.16) with $\xi(0)=1-c$.

From this solution formula, we can establish the solution’s large-time
behaviour. We set

$\eta=\bar{\eta}(t)+\eta_{e}(t),\quad\bar{\eta}(t)=-b\int_{0}^{t}e^{-a^{2}s}q^{b}(s)\,ds.$

######  Lemma 2.3.

  1. (i)

The function $\bar{\eta}$ is positive, monotonically decreasing, and
converging to $-b/(b+c)$ as $t\to\infty$.

  2. (ii)

The estimate

$\left|\frac{\eta_{e}}{c(\mu+b)}\right|\leq
e^{-a^{2}t}\int_{t}^{\infty}q^{b}(s)\,ds\leq e^{-a^{2}t}\frac{1}{b}$

holds for $t>0$. In particular,

$\lim_{t\to\infty}\eta_{e}e^{a^{2}t}=0.$

###### Proof.

  1. (i)

Since $q^{b}\geq 0$, the monotonicity is clear. We observe that

$\int_{0}^{\infty}e^{-a^{2}s}q^{b}(s)\,ds=\mathcal{L}\left[q^{b}\right](0+a^{2})=\frac{1}{\sqrt{a^{2}}+b}$

by (2.13). Thus $\lim_{t\to\infty}\bar{\eta}(t)=-b/(b+a)$.

  2. (ii)

Since

$\int_{0}^{\infty}e^{-\mu^{2}s}q^{b}(s)\,ds=\frac{1}{\mu+b},$

we observe that

$\displaystyle\eta_{e}$
$\displaystyle=-c(\mu+b)\int_{t}^{\infty}e^{-\mu^{2}s}q^{b}(s)\,ds\,e^{-(a^{2}-\mu^{2})t}$
$\displaystyle=-c(\mu+b)\int_{t}^{\infty}e^{\mu^{2}(t-s)}q^{b}(s)\,ds\,e^{-a^{2}t}.$

Since $e^{\mu^{2}(t-s)}\leq 1$ for $s\geq t$, this implies
$\left|\eta_{e}/c(\mu+b)\right|\leq
e^{-a^{2}t}\int_{t}^{\infty}q^{b}(s)\,ds\leq
e^{-a^{2}t}\int_{0}^{\infty}q^{b}(s)\,ds=e^{-a^{2}t}/b$. The proof is now
complete.

∎

######  Remark 2.4.

We consider the initial-boundary value problem for

$\left\\{\begin{aligned}
w_{t}-(-x)^{\alpha}w_{xx}&=0&&\text{in}\quad(-\infty,0)\times(0,\infty),\\\
w(0,t)&=\eta(t)&\quad&\text{for}\quad t>0,\quad\eta(0)=0,\\\
\lim_{x\to-\infty}w(x,t)&=0&&\text{for}\quad t>0\\\ w(x,0)&=0&&\text{for}\quad
x<0.\end{aligned}\right.$

Then, as in [3], we obtain that

$w_{x}(0,t)=c_{\gamma}\partial_{t}^{\gamma}\eta,$

where $\gamma=1/(2-\alpha)$ with some constant $c_{\gamma}>0$ provided that
$\alpha<1$. Indeed as in [3], let $\psi$ be a solution of

$\left\\{\begin{aligned}
\psi-(-x)^{\alpha}\psi_{xx}&=0\quad\text{in}\quad(-\infty,0),\\\
\psi(0)&=1,\\\ \lim_{x\to-\infty}\psi(x)&=0.\end{aligned}\right.$

Since the Laplace transform $\widehat{w}$ of $w$ satisfies

$\lambda\widehat{w}-(-x)^{\alpha}\widehat{w}_{xx}=0,\quad\widehat{w}(0,\lambda)=\widehat{\eta}(\lambda)$

we see, by scaling, that

$\widehat{w}(x,\lambda)=\widehat{\eta}(\lambda)\psi\left(\lambda^{1/(2-\alpha)}x\right).$

Thus

$\displaystyle\partial_{x}\widehat{w}(0,\lambda)$
$\displaystyle=\lambda^{1/(2-\alpha)}\psi^{\prime}(0)\widehat{\eta}(\lambda)$
$\displaystyle=\psi^{\prime}(0)\lambda^{\gamma-1}\widehat{\eta_{t}}(\lambda)$

since $f(0)=0$. Thus

$w_{x}(0,t)=\psi^{\prime}(0)f_{1-\gamma}^{0}*\eta_{t}.$

If $\gamma<1$, then $f_{1-\gamma}^{0}$ is integrable. As noted in [3],
$\psi^{\prime}(0)$ exists (even for the degenerate case, i.e., $\alpha>0$) and
$\psi^{\prime}(0)>0$. Thus,

$w_{x}(0,t)=c_{\gamma}\partial_{t}^{\gamma}\eta\quad\text{with}\quad
c_{\gamma}=\psi^{\prime}(0)$

at least for $\gamma=1/(2-\alpha)<1$, i.e., $\alpha<1$.

######  Remark 2.5.

The reader might be interested in how fractional partial differential
equations like fractional diffusion equations are derived. We consider

$\left\\{\begin{aligned}
w_{t}-(-x)^{\alpha}w_{xx}&=0&\quad\text{in}\quad&(-\infty,0)\times\mathbf{R}^{n-1}\times(0,\infty),\\\
w_{x}-\Delta_{y}w&=f&\quad\text{on}\quad&(-\infty,0)\times\mathbf{R}^{n-1}\times(0,\infty),\\\
w(x,y,0)&=0&\quad\text{on}\quad&(-\infty,0)\times\mathbf{R}^{n-1},\end{aligned}\right.$

where $f=f(y,t)$ is a given function. Then by Remark 2.4, the equation for
$\eta(y,t)=w(0,y,t)$ is formally obtained as

(2.19)
$c_{\gamma}\partial_{t}^{\gamma}\eta-\Delta_{y}\eta=f\quad\text{in}\quad\mathbf{R}^{n-1}\times(0,\infty).$

This type of equation is a kind of fractional diffusion that has been well-
studied; see [15, 29]. Here, we briefly recall only the well-posedness of its
initial-boundary value problem for (2.19) in a domain. In the framework of
distributions, the well-posedness of its initial-boundary value problems has
been established in [23, 28] by using the Galerkin method. The unique
existence of viscosity solutions for (2.19), including general nonlinear
problems, has been established in [7, 20] and also in [26] for the whole space
$\mathbf{R}^{n-1}$. The scope of equations these theories apply is different.
However, it has been proved in [6] that two notions of solutions (viscosity
solution and distributional solution) agree for (2.19) when we consider the
Dirichlet problem in a smooth bounded domain.

##  3\. Convergence

We shall prove Lemma 1.4.

###### Proof of Lemma 1.4.

For a function $f$ on $(-L,L)$, we decompose it into its odd and even parts,
i.e.,

$f_{\mathrm{odd}}(x):=\frac{f(x)-f(-x)}{2},\quad
f_{\mathrm{even}}(x):=\frac{f(x)+f(-x)}{2},$

so that $f=f_{\mathrm{odd}}+f_{\mathrm{even}}$. By the structure of the
equation, $V_{\mathrm{odd}}^{\varepsilon}$ and
$V_{\mathrm{even}}^{\varepsilon}$ solve the equation (1.10) separately.

At first glance, the locally uniform convergence follows from the maximum or
comparison principles for a linear parabolic equation [22]. However, a direct
application of the maximum principle is impossible since the domains of
functions $V^{\varepsilon}$ and $V$ are different. We first show the
convergence where initial data is smooth.

For the odd part, the term $\partial_{y}\\{1_{y>0}\\}V$ does not affect since
$V=0$ at $x=0$. Thus, the equation (1.10) is reduced to

(3.1)
$\left\\{\begin{array}[]{rll}\tau_{1}V_{t}&=V_{yy}-a^{2}(V-1),&|y|<L/\varepsilon\\\
\left.V\right|_{t=0}&=V_{0}^{\varepsilon},&V_{y}(\pm
L/\varepsilon,t)=0\quad\text{for}\quad t>0,\end{array}\right.$

where $V_{0}^{\varepsilon}(y)=V^{\varepsilon}(y,0)$ for $|y|<L/\varepsilon$.
Let $V^{\varepsilon}$ be its solution.

We extend an odd function $V_{0}^{\varepsilon}$ outside
$\left(L/\varepsilon,3L/\varepsilon\right)$ for $x$ to be “even” with respect
to $L/\varepsilon$, i.e.,

$\widetilde{V_{0}^{\varepsilon}}\left(x-L/\varepsilon\right)=\widetilde{V_{0}^{\varepsilon}}\left(L/\varepsilon-x\right),\quad
x\in\left(L/\varepsilon,3L/\varepsilon\right),$

where $\widetilde{V_{0}^{\varepsilon}}$ is its extension. We extend outside
$\left(-L/\varepsilon,3L/\varepsilon\right)$ so that the extension
$\overline{V_{0}^{\varepsilon}}$ is periodic in $\mathbf{R}$ with period
$4L/\varepsilon$. Since $V_{0}^{\varepsilon}$ is even with respect to
$L/\varepsilon$, $\left(\overline{V_{0}^{\varepsilon}}\right)_{y}\left(\pm
L/\varepsilon\right)=0$ if $V_{0y}^{\varepsilon}\left(\pm
L/\varepsilon\right)=0$ and smooth. Solution $V^{\varepsilon}$ is the
restriction on $\left(-L/\varepsilon,L/\varepsilon\right)$ of a solution $W$
of

(3.2) $\left\\{\begin{array}[]{rl}&\tau_{1}W_{t}=W_{yy}-a^{2}(W-1),\quad
y\in\mathbf{R}\\\ &W|_{t=0}=\overline{V_{0}^{\varepsilon}}.\end{array}\right.$

Although the maximum principle implies

$\left\lVert
W-V^{\varepsilon}\right\rVert_{L^{\infty}(\mathbf{R})}(t)\leq\left\lVert\overline{V_{0}^{\varepsilon}}-V_{0}\right\rVert_{L^{\infty}(\mathbf{R})},$

our assumption of the convergence $V_{0}^{\varepsilon}\to V_{0}$ does not
guarantee
$\left\lVert\overline{V_{0}^{\varepsilon}}-V_{0}\right\rVert_{L^{\infty}(\mathbf{R})}\to
0$. We argue differently.

We approximate $V_{0}$ by $V_{0\delta}=V_{0}*\rho_{\delta}$, where
$\rho_{\delta}$ is a symmetric mollifier. We also approximate
$V_{0}^{\varepsilon}$ by

$V_{0\delta}^{\varepsilon}=\overline{V_{0}^{\varepsilon}}*\rho_{\delta}.$

Let $V_{\delta}^{\varepsilon}$ be the restriction of
$W_{\delta}^{\varepsilon}$ of (LABEL:ENoS2) with initial data
$V_{0\delta}^{\varepsilon}$. It follows from the parity and periodic condition
that this $V_{\delta}^{\varepsilon}$ solves (3.1); moreover,
$W_{\delta}^{\varepsilon}=W^{\varepsilon}*\rho_{\delta}$, where
$W^{\varepsilon}$ is the solution of (LABEL:ENoS2) with initial data
$\overline{V_{0}^{\varepsilon}}$. Let $V_{\delta}$ be the bounded solution of
(LABEL:ENoS2) with initial data $V_{0\delta}$ and again
$V_{\delta}=V*\rho_{\delta}$. These properties follow from the fact that our
equation is of constant coefficients. For fixed $\delta>0$, we observe that

$V_{\delta}^{\varepsilon}\to V_{\delta}\quad\text{in}\quad
L^{\infty}\left((-M,M)\times(0,T)\right)$

for $M>0$. Indeed, by the maximum principle

$\displaystyle\left\lVert
V_{\delta}^{\varepsilon}-1\right\rVert_{\infty,\varepsilon}(t)$
$\displaystyle\leq\left\lVert
V_{0\delta}^{\varepsilon}-1\right\rVert_{\infty,\varepsilon},$
$\displaystyle\left\lVert\partial_{y}V_{\delta}^{\varepsilon}\right\rVert_{\infty,\varepsilon}(t)$
$\displaystyle\leq\left\lVert\partial_{y}V_{0\delta}^{\varepsilon}\right\rVert_{\infty,\varepsilon},$
$\displaystyle\left\lVert\partial_{t}V_{\delta}^{\varepsilon}\right\rVert_{\infty,\varepsilon}(t)$
$\displaystyle\leq\left\lVert\partial_{t}V_{0\delta}^{\varepsilon}\right\rVert_{\infty,\varepsilon},$

where we interpret $\tau_{1}\partial_{t}V_{0\delta}^{\varepsilon}=V_{0\delta
yy}^{\varepsilon}-a^{2}\left(V_{0\delta}^{\varepsilon}-1\right)$ and
$\|\cdot\|_{\infty,\varepsilon}$ is the sup norm on
$\left(-L/\varepsilon,L/\varepsilon\right)$. Because of the mollifier, the
right-hand side is bounded by a constant multiple of $\left\lVert
V_{0}^{\varepsilon}\right\rVert_{\infty,\varepsilon}$, which is uniformly
bounded in $\varepsilon<1$. By the Arzelà–Ascoli theorem, $V^{\varepsilon}$
converges (locally uniformly in $\mathbf{R}\times[0,\infty)$) to a bounded
(weak) solution to (LABEL:ENoS2) with initial data $V_{0}$ by taking a
subsequence. Since $V$ is bounded, by the uniqueness of the limit problem, the
convergence is now full (without taking a subsequence). Note that we only
invoke the locally uniform convergence of $V_{0}^{\varepsilon}$ to $V_{0}$
other than uniform bound on derivatives.

We note that

$V-V^{\varepsilon}=V-V_{\delta}+V_{\delta}-V_{\delta}^{\varepsilon}+V_{\delta}^{\varepsilon}-V^{\varepsilon}$

and observe that

$\|V-V^{\varepsilon}\|(t)\leq\|V-V_{\delta}\|(t)+\|V_{\delta}-V_{\delta}^{\varepsilon}\|(t)+\|V_{\delta}^{\varepsilon}-V^{\varepsilon}\|(t)=:I+II+III$

where the norm is $\|\cdot\|$ taken in $L^{\infty}(-M,M)$ for $M>0$. By the
maximum principle,

$\displaystyle\|V_{\delta}^{\varepsilon}-V^{\varepsilon}\|(t)\leq\|V_{0\delta}^{\varepsilon}-V_{0}^{\varepsilon}\|_{L^{\infty}(-L/\varepsilon,L/\varepsilon)}$
$\displaystyle\text{for}\quad L>\varepsilon M,$
$\displaystyle\|V-V_{\delta}\|(t)\leq\|V_{0\delta}-V_{0}\|_{L^{\infty}(\mathbf{R})}$
$\displaystyle\text{for}\quad L>\varepsilon M.$

Since

$\displaystyle
V_{0\delta}^{\varepsilon}=\rho_{\delta}*\left(\overline{V_{0}^{\varepsilon}}-V_{0}\right)+\rho_{\delta}*V_{0},$
$\displaystyle
V_{0}^{\varepsilon}=\overline{V_{0}^{\varepsilon}}-V_{0}+V_{0},$

and $\|\rho_{\delta}*f\|_{\infty}\leq\|f\|_{\infty}$, we see that

$\|V_{0\delta}^{\varepsilon}-V_{0}^{\varepsilon}\|_{L^{\infty}(-L/\varepsilon,L/\varepsilon)}\leq
2\|\overline{V_{0}^{\varepsilon}}-V_{0}\|_{\infty,\varepsilon}+\|\rho_{\delta}*V_{0}-V_{0}\|_{\infty}.$

Thus

$\sup_{0<t<T}III\leq
2\|V_{0}^{\varepsilon}-V_{0}\|_{\infty,\varepsilon}+\|\rho_{\delta}*V_{0}-V_{0}\|_{\infty}.$

Fixing $\delta>0$ and sending $\varepsilon\to 0$, we observe that

$\varlimsup_{\varepsilon\downarrow 0}\sup_{0<t<T}(I+II+III)\leq
2\|\rho_{\delta}*V_{0}-V_{0}\|_{\infty}$

since $V_{\delta}^{\varepsilon}\to V_{\delta}$ in
$L^{\infty}\left((-M,M)\times[0,T]\right)$. Sending $\delta\downarrow 0$, we
obtain

$\varlimsup_{\varepsilon\downarrow 0}\sup_{0<t<T}\|V-V^{\varepsilon}\|(t)=0$

since $V_{0}$ is uniformly continuous.

We next study the even part. The general strategy is the same but more
involved than the odd part. For the even part, we first note that
$V_{\mathrm{even}}^{\varepsilon}$ solves 2 &τ_1 V_t = V_yy - a^2(V-1), y
∈I_ε:= ( -L/ε, 0 ), t¿0,
V_y(0,t) + bV(0,t) = 0, V_y (-L/ε, t ) = 0, t¿0,
V—_t=0(y) = V_0even^ε, where $V_{0}^{\varepsilon}(y)=V^{\varepsilon}(y,0)$ for
$y\in I_{\varepsilon}$. We suppress the word “even” from now on. We shall
approximate $V_{0}^{\varepsilon}$ by a smoother function
$V_{0\delta}^{\varepsilon}$ and approximate $V_{0}=\lim_{\varepsilon\to
0}V_{0}^{\varepsilon}$ by a smoother function uniformly. There are many
possible ways, and we rather like an abstract way. Let $BUC(I_{\varepsilon})$
denote the space of all bounded uniformly continuous functions in
$\overline{I_{\varepsilon}}$. It is a Banach space equipped with the norm

$\lVert f\rVert_{\infty,\varepsilon}=\sup\left\\{\left|f(x)\right|\bigm{|}x\in
I_{\varepsilon}\right\\}.$

If $\varepsilon=0$, then $I_{\varepsilon}$ should be interpreted as
$(-\infty,0]$. Let $A$ be the operator on $BUC(I_{0})$ defined by

$Af:=(-\partial_{y}^{2}+a^{2})f\quad\text{(in the distribution sense)}$

with the

$D(A)=\left\\{f\in BUC(I_{0})\bigm{|}Af\in BUC(I_{0}),\
f_{y}(0)+bf(0)=0\right\\}.$

A standard theory [16] implies that $-A$ generates an analytic semigroup
$e^{-tA}$ in $BUC(I_{0})$. In particular,

$\displaystyle\left\lVert Ae^{-tA}f\right\lVert_{\infty,0}$ $\displaystyle\leq
Ct^{-1}\lVert f\rVert_{\infty,0},$ $\displaystyle\left\lVert
e^{-tA}f\right\lVert_{\infty,0}$ $\displaystyle\leq C\lVert
f\rVert_{\infty,0}$

with some constant $C>0$ independent of time $t\in(0,1)$ and $f\in
BUC(I_{0})$. For a function $h\in BUC(I_{\varepsilon})$, we extend it to
$\widetilde{h}$ so that $\widetilde{h}(x)=h\left(-L/\varepsilon\right)$ for
$x<-L/\varepsilon$. For $V_{0}$, we set $V_{0\delta}=e^{-\delta A}V_{0}$. For
$V_{0}^{\varepsilon}$, we tempt to set $V_{0\delta}^{\varepsilon}=e^{-\delta
A}\widetilde{V_{0}^{\varepsilon}}$. However, unfortunately,
$V_{0\delta}^{\varepsilon}$ does not satisfy the boundary condition at
$-L/\varepsilon$ though it satisfies
$\left(\partial_{y}V_{0\delta}^{\varepsilon}+bV_{0\delta}^{\varepsilon}\right)(0)=0$
and is $C^{2}$ (actually smooth). We set
$\sigma(x)=\rho_{\delta^{\prime}}*\left(1-|x|\right)_{+}$ for a fixed
$\delta^{\prime}>0$ so that $\sigma^{\prime}\left(1/2\right)=-\kappa$ with
$\kappa<0$. We set
$\sigma_{\delta^{\prime\prime}}(y)=\delta^{\prime\prime}\sigma\left(y/\delta^{\prime\prime}\right)$
for small $\delta^{\prime\prime}>0$. For a given $h\in C^{2}(-\infty,0]$, we
modify

$h_{\delta^{\prime\prime}}(y)=h(y)+c\sigma_{\delta^{\prime\prime}}(y-\nu)$

where we take $c$ so that $c\kappa=h_{y}\left(-L/\varepsilon\right)$,
$-L/\varepsilon-\nu=\delta^{\prime\prime}/2$. By this modification,

$\partial_{y}h_{\delta^{\prime\prime}}\left(-L/\varepsilon\right)=0$

and $h_{\delta^{\prime\prime}}\to h$ in $L^{\infty}(-\infty,0)$ as
$\delta^{\prime\prime}\to 0$. We set

$V_{0\delta}^{\varepsilon}=\left(e^{-\delta
A}\widetilde{V_{0}^{\varepsilon}}\right)_{\delta}.$

Since $V_{0\delta}^{\varepsilon}$ satisfies the boundary condition on the
boundary of $I_{\varepsilon}$ and is smooth, we observe that
$\partial_{t}V_{\delta}^{\varepsilon}$ is continuous up to the boundary of
$I_{\varepsilon}\times[0,T)$, where $V_{\delta}^{\varepsilon}$ denotes the
solution of (3) with initial data $V_{0\delta}^{\varepsilon}$. By the maximum
principle,

(3.3) $\displaystyle\begin{aligned} &\lVert
V_{\delta}^{\varepsilon}-1\rVert_{\infty,\varepsilon}(t)\leq\lVert
V_{0\delta}^{\varepsilon}-1\rVert_{\infty,\varepsilon}\\\
&\lVert\partial_{y}V_{\delta}^{\varepsilon}\rVert_{\infty,\varepsilon}(t)\leq\lVert\partial_{y}V_{0\delta}^{\varepsilon}\rVert_{\infty,\varepsilon}+b\lVert
V_{\delta}^{\varepsilon}\rVert_{\infty,\varepsilon}\\\
&\lVert\partial_{t}V_{\delta}^{\varepsilon}\rVert_{\infty,\varepsilon}(t)\leq\lVert\partial_{t}V_{0\delta}^{\varepsilon}\rVert_{\infty,\varepsilon},\end{aligned}$

where $\partial_{t}V_{0\delta}^{\varepsilon}$ for initial data
$V_{0\delta}^{\varepsilon}$ should be interpreted as in the proof for the odd
part. The term involving $b$ appears because of the Robin type boundary
condition. As in the case for (LABEL:ENoS2), by the Arzelà–Ascoli theorem and
the uniqueness of the limit equation, we can prove that
$V_{\delta}^{\varepsilon}$ converges to $V_{\delta}$ locally uniformly in
$\mathbf{R}\times[0,\infty)$. Note that for a fixed $\delta>0$, the right-hand
sides of (LABEL:EMAX) are uniformly bounded as $\varepsilon\to 0$ since
$V_{0}^{\varepsilon}$ converges to $V_{0}$ uniformly. The comparison principle
implies that

$\displaystyle\|V_{\delta}^{\varepsilon}-V^{\varepsilon}\|_{\infty,\varepsilon}(t)\leq\|V_{0\delta}^{\varepsilon}-V_{0}^{\varepsilon}\|_{\infty,\varepsilon},$
$\displaystyle t>0$
$\displaystyle\|V-V_{\delta}\|_{\infty,0}(t)\leq\|V_{0\delta}-V_{0}\|_{\infty,0}$
$\displaystyle t>0.$

Thus

$\displaystyle\|V-V^{\varepsilon}\|(t)$
$\displaystyle\leq\|V-V_{\delta}\|(t)+\|V_{\delta}-V_{\delta}^{\varepsilon}\|(t)+\|V_{\delta}^{\varepsilon}-V^{\varepsilon}\|(t)$
$\displaystyle\leq\|V_{0}-V_{0\delta}\|_{\infty,0}+\|V_{\delta}-V_{\delta}^{\varepsilon}\|(t)+\|V_{0\delta}^{\varepsilon}-V_{0}^{\varepsilon}\|_{\infty,\varepsilon},$

where the norm $\|\cdot\|$ is taken in $L^{\infty}(0,M)$ for $M>0$. Taking the
supremum in $t\in(0,T)$ and sending $\varepsilon\to 0$, we obtain that

$\sup_{0<t<T}\|V-V^{\varepsilon}\|(t)\leq 2\|V_{0\delta}-V_{0}\|_{\infty}$

since we know $\sup_{0<t<T}\|V_{\delta}-V_{\delta}^{\varepsilon}\|(t)\to 0$ as
$\varepsilon\to 0$. Sending $\delta\to 0$, we conclude that $V^{\varepsilon}$
converges to $V$ locally uniformly in $[0,\infty)\times[0,\infty)$.

Since we know that $V^{\varepsilon}$ is continuous up to $y=0$ and $t=0$, this
gives the local uniform convergence of $V^{\varepsilon}(0,t)$ in $[0,\infty)$.
The proof is now complete. ∎

If we only assume that the initial data for (LABEL:ENoS2) and (3) is bounded
and Lipschitz, we have a similar estimate in (LABEL:EMAX) up to the first
derivative of the solution. However, the estimate for the time derivative
should be altered. Since we used such an estimate in Lemma 2.2, we state it in
the case of $\varepsilon=0$ for the reader’s convenience.

######  Lemma 3.1.

Let $V$ be the bounded solution of (1.10) in $\mathbf{R}\times(0,\infty)$ with
a bounded and Lipschitz continuous initial data $V_{0}$. Then for each $T>0$,
there is a constant $C$ depending only on $a$, $b$, and $T$ such that

$t^{1/2}\|\partial_{t}V\|_{L^{\infty}(\mathbf{R})}(t)\leq
C\left(\|\partial_{y}w_{0}\|_{L^{\infty}(\mathbf{R})}+\|w_{0}\|_{L^{\infty}(\mathbf{R})}+1\right)\quad\text{for}\quad
t\in(0,T).$

###### Proof.

We give direct proof. We may assume that $\tau_{1}=1$. We set $w=V-1$ and
$u=e^{a^{2}t}w$ to get

$u_{t}=u_{xx}-2b\partial_{x}\\{1_{x>0}\\}\left(u+e^{a^{2}t}\right),$

where we denote by $x$ instead of $y$. We consider this equation with initial
data $w_{0}=V_{0}-1$. It suffices by simple scaling
$u_{\lambda}(x,t)=u(\lambda x,\lambda^{2}t)$ to prove the desired estimate for
some $T$ independent of $w_{0}$.

Let $E(x,t)$ be the Gauss kernel as before. Then the solution can be
represented as

(3.4)
$u(x,t)=(E\ast_{x}w_{0})(x,t)-\int_{0}^{t}E(x,t-\tau)h(\tau)\,d\tau,\quad
h(t)=2b\left(u(0,t)+e^{a^{2}t}\right),$

where

$(E\ast_{x}w_{0})(x,t)\coloneqq\int_{-\infty}^{\infty}E(x-y,t)w_{0}(y)\,dy.$

Since we can approximate a smooth $w_{0}$, establishing

(3.5) $\|\partial_{t}u\|_{L^{\infty}(\mathbf{R})}(t)\leq
Ct^{-1/2}\left(\|\partial_{x}w_{0}\|_{L^{\infty}(\mathbf{R})}+\|w_{0}\|_{L^{\infty}(\mathbf{R})}+1\right)\quad\text{for}\quad
t\in(0,T)$

with some positive constants $C$ and $T$ independent of $w_{0}$ suffices,
assuming that $\partial_{t}u$ exists and is bounded in $\mathbf{R}\times(0,T)$
for small $T$. By the maximum principle (LABEL:EMAX) and the corresponding
estimate for the odd part, we know that

(3.6) $\|u\|_{L^{\infty}(\mathbf{R})}(t)\leq
c\left(1+\|w_{0}\|_{L^{\infty}(\mathbf{R})}\right)$

with $c$ independent of $w_{0}$ and $t>0$. We estimate $\partial_{t}u$ in
(3.4). Since
$\|f\ast_{x}g\|_{L^{\infty}(\mathbf{R})}\leq\|f\|_{L^{1}(\mathbf{R})}\|g\|_{L^{\infty}(\mathbf{R})}$
and

$\partial_{t}(E\ast_{x}w_{0})=(\partial_{x}E)\ast_{x}\partial_{x}w_{0},$

we easily see (cf. [5, Chapter 1]) that

(3.7)
$\left\|\partial_{t}(E*_{x}w_{0})\right\|_{L^{\infty}(\mathbf{R})}(t)\leq\frac{c^{\prime}}{t^{1/2}}\|\partial_{x}w_{0}\|_{L^{\infty}(\mathbf{R})}\quad\text{for}\quad
t>0$

with $c^{\prime}$ independent of $w_{0}$. The second term of the right-hand
side of (3.4) is more involved than the first term because $h$ contains $u$.
We observe that

$\partial_{t}\int_{0}^{t}E(x,t-\tau)h(\tau)\,d\tau=\int_{0}^{t}E(x,t-\tau)\partial_{\tau}h(\tau)\,d\tau+E(x,t)h(0).$

Since $|E|\leq(4\pi t)^{-1/2}$, it holds that

$\left|\int_{0}^{t}E(x,t-\tau)\partial_{\tau}h(\tau)\,d\tau\right|\leq\int_{0}^{t}\frac{1}{(4\pi(t-\tau))^{1/2}}\frac{1}{\tau^{1/2}}\,d\tau\cdot\sup_{0<t<T}|t^{1/2}\partial_{t}h(t)|$

for $t\in(0,T)$.

Thus,

$\displaystyle\sup_{0<t<T}\left\|\partial_{t}\int_{0}^{t}E(x,t-\tau)h(\tau)\,d\tau\right\|_{\infty}$
$\displaystyle\leq\,C_{1}\sup_{0<t<T}\left|t^{1/2}\partial_{t}h(t)\right|+(4\pi
t)^{-1/2}\left|h(0)\right|$
$\displaystyle\leq\,C_{2}\left(\sup_{0<t<T}\|t^{1/2}\partial_{t}u\|_{L^{\infty}(\mathbf{R})}(t)+1\right)+C_{3}t^{-1/2}\left(\|u\|_{L^{\infty}(\mathbf{R})}(t)+1\right)$

with some constants $C_{j}$ ($j=1,2,3$). By (3.4) and (3.7), we now observe
that

$\sup_{0<t<T}\left\|t^{1/2}\partial_{t}u\right\|_{L^{\infty}(\mathbf{R})}(t)\leq
C_{4}\left\|\partial_{x}w_{0}\right\|_{L^{\infty}(\mathbf{R})}\\\
+C_{2}T^{1/2}\left(\sup_{0<t<T}\|t^{1/2}\partial_{t}u\|_{L^{\infty}(\mathbf{R})}+1\right)+C_{3}\sup_{0<t<T}\left(\|u\|_{L^{\infty}(\mathbf{R})}(t)+1\right)$

with $C_{4}$ independent of $w_{0}$, $u$, and $T$. Applying estimate for
$\|u\|_{\infty}$ in (3.6), we conclude (3.5) for sufficiently small $T$ by
absorbing
$C_{2}T^{1/2}\sup_{0<t<T}\|t^{1/2}\partial_{t}u\|_{L^{\infty}(\mathbf{R})}$ in
the right-hand side to the left. ∎

##  4\. Dirichlet condition for the total variation flow

In this section, we recall a notion of total variation flow for a given $v$
and prove Lemma 1.1. We consider

(4.1) $\partial_{t}u=\operatorname{div}\left(\beta\nabla u/|\nabla
u|\right)\quad\text{in}\quad I\times(0,T),$

where $\beta\in C\left(I\times[0,T]\right)$ is a given nonnegative function;
here, $I=(p_{0},p_{1})$ is an open interval and $T>0$. If we impose the
Dirichlet boundary condition

(4.2) $u=g\quad\text{on}\quad\partial I,$

(4.1) with (4.2) should be interpreted as an $L^{2}$-gradient flow of a time-
dependent total variation type energy

$\Phi^{t}(u)=\int_{I}\beta(x,t)|u_{x}|+\sum_{i=0}^{1}|\gamma
u-g|(p_{i})\beta(p_{i},t)$

when $\int\beta|u_{x}|$ is a weighted total variation of $\beta$ and $\gamma
u$ is a trace of $u$ on $\partial I$. We consider this energy in $L^{2}(I)$ by
$\Phi^{t}(u)=\infty$ when $\Phi^{t}(u)$ is not finite. It is clear that
$\Phi^{t}$ is convex in $L^{2}(I)$. If $\beta$ is spatially constant, it is
well known that $\Phi^{t}$ is also lower semicontinuous; see, e.g. [1]. The
solution of (4.1) with (4.2) should be interpreted as the gradient flow of
form

(4.3) $u_{t}\in-\partial\Phi^{t}(u),$

where $\partial\Phi^{t}$ denotes the subdifferential of $\Phi^{t}$ in
$L^{2}(I)$, i.e.,

$\partial\Phi^{t}(u)=\left\\{f\in
L^{2}(I)\biggm{|}\Phi^{t}(u+h)-\Phi^{t}(u)\geq\int_{I}hf\,dx\ \text{for all}\
h\in L^{2}(I)\right\\}.$

It is standard that (4.3) is uniquely solvable for given initial data
$u_{0}\in L^{2}(I)$ if $\Phi^{t}$ does not depend on time and is lower
semicontinuous and convex on the Hilbert space $L^{2}(I)$ (see, for instance,
[14, 2]). It applies to the total variation flow case when $\beta$ is a
content. In this case, the subdifferential becomes

$\partial\Phi(u)=\left\\{v\in L^{2}(I)\biggm{|}v=-(\beta z)_{x},\
\|z\|_{\infty}\leq 1,\ \Phi(u)=\int_{I}uv\,dx\right\\},$

when $\Phi=\Phi^{t}$; see [1]. The equation (4.3) is

$u_{t}=\operatorname{div}\beta z$

with $|z|\leq 1$ in $I$ and $\Phi(u)=\int_{I}u\operatorname{div}\beta z\,dx$.
We mimic this notion of the solution. A function $u\in
C\left([0,T),L^{2}(I)\right)$ is a solution to (4.1) with (4.2) if there is
$z\in L^{\infty}\left(I\times(0,T)\right)$ such that

(4.4) $\displaystyle u_{t}=(\beta z)_{x}\quad\text{in}\quad I\times(0,T)$
(4.5) $\displaystyle|z|\leq 1\quad\text{in}\quad I\times(0,T)$ (4.6)
$\displaystyle-\int_{I}u(\beta
z)_{x}\,dx=\int_{I}\beta|u_{x}|+\sum_{i=0}^{1}|\gamma u-g|\beta(p_{i},t).$

Under this preparation, we shall prove Lemma 1.1.

###### Proof of Lemma 1.1.

Since $u_{t}^{b}=0$, (4.4) says that $\beta z$ is a constant $c$. The
condition (4.5) is equivalent to saying that $|c|\leq\beta$. Since

$\int_{-L}^{L}u(\beta z)_{x}\,dx=\int_{0}^{L}b(\beta
z)_{x}\,dx=-\beta(0)z(0)b,$

(4.6) is equivalent to

$\beta(0)z(0)b=\beta(0)b.$

In other words, $c$ must be $\beta(0)$. Thus, the existence of $z$ satisfying
$|z|\leq 1$ is guaranteed if and only if

$c=\beta(0)\leq\beta(x)\quad\text{for all}\quad x\in(-L,L).$

The equation (4.6) is fulfilled by taking $g=b$. Thus $u^{b}$ is a stationary
solution to (1) with (1.3). ∎

##  5\. Numerical experiment

In this section, we calculate the solution of (1.7)–(1.9) with
$v_{0}(x)=1-ce^{-a|x|/\varepsilon}$ and compare its value at $x=0$ with an
explicit solution of (2.16) whose explicit form is given in (2.18).

###  5.1. Numerical scheme

Since the initial function, $v_{0}$, is an iven function, the original problem
(1.7)–(1.9) is reduced to

$\displaystyle\begin{dcases*}\frac{\tau_{1}}{\varepsilon}v_{t}=\varepsilon
v_{xx}-\frac{a^{2}(v-1)}{\varepsilon}&in $(0,L)\times(0,\infty)$,\\\
-v_{x}(0,t)+\frac{b}{\varepsilon}v(0,t)=0,&for $t>0$,\\\ v_{x}(L,t)=0,&for
$t>0$,\\\ v(x,0)=v_{0}(x),&for $x\in[0,L]$.\end{dcases*}$

The computational region $[-L,L]$ is divided into unfiorm mesh partitions:

$\displaystyle x_{i}=i\Delta x,\quad i=-1,0,\ldots,N,N+1,\quad\Delta
x=\frac{L}{N}.$

The points $x_{-1}$ and $x_{N+1}$ are needed to handle the Neumann boundary
conditions.

The approximation of $v$ at $x=x_{i}$ is written as $v_{i}$. The central
finite difference approximates the Laplace operator, and the time derivative
is approximated by the backward difference, yielding the folliwng linear
system

$\displaystyle\frac{\tau_{1}}{\varepsilon}\frac{v_{i}-\hat{v}_{i}}{\Delta
t}=\varepsilon\frac{v_{i-1}-2v_{i}+v_{i+1}}{(\Delta
x)^{2}}-\frac{a^{2}(v_{i}-1)}{\varepsilon},\quad i=0,1,\ldots,N,$

where $\hat{v}_{i}$ is the value known at the current time, and $v_{i}$ is the
value to be found at the next time. The Neumann boundary conditions at $x=0$
and $x=L$ are approximated as

$\displaystyle-\frac{v_{1}-v_{-1}}{2\Delta
x}+\frac{b}{\varepsilon}v_{0}=0\quad\text{and}\quad\frac{v_{N+1}-v_{N-1}}{2\Delta
x}=0$

by the central finite differences.

###  5.2. Results

Some results are shown for different values of $c$ for parameters

$\displaystyle L=1,\quad a=1,\quad N=200,\quad\Delta t=(\Delta x)^{2},\quad
b=1.$

The results of the numerical experiments are summarised in Figure 2.

(a)

(b)

(c)

$\varepsilon$ $c$ | $0$ | $1/4$ | $2$
---|---|---|---
$1$ | $6.8\cdot 10^{-2}$ | $6.8\cdot 10^{-2}$ | $6.7\cdot 10^{-2}$
$2^{-1}$ | $9.1\cdot 10^{-3}$ | $9.1\cdot 10^{-3}$ | $9.3\cdot 10^{-3}$
$2^{-2}$ | $1.4\cdot 10^{-4}$ | $1.5\cdot 10^{-4}$ | $2.6\cdot 10^{-4}$
$2^{-3}$ | $4.6\cdot 10^{-5}$ | $4.6\cdot 10^{-5}$ | $4.6\cdot 10^{-5}$
(d)

Figure 2. Results of numerical experiments: (a) $c=0$, (b) $c=1/4$, (c) $c=2$,
(d) table of $L^{\infty}$-errors for different $\varepsilon$ and $c$.

In (a), (b), and (c), the horizontal axis represents time, and the vertical
axis represents the value at the origin. As $\varepsilon$ is decreased, the
numerical solution converges to the exact solution to the extend that the
exact and numerical solutions overlap. Indeed, the table of
$L^{\infty}$-errors for different values of $\varepsilon$ and $c$ is shown in
(d). The errors for $\varepsilon=2^{-3}$ are of order $10^{-5}$, indicating
that the solution for small $\varepsilon$ is an excellent approximation to the
solution of the fractional time differential equation obtained as the singular
limit.

## Acknowledgements

The work of the first author was partly supported by JSPS KAKENHI Grant
Numbers JP19H00639 and JP20K20342, and by Arithmer Inc., Ebara Corporation,
and Daikin Industries, Ltd. through collaborative grants. The work of the
third author was partly supported by JSPS KAKENHI Grant Number JP20K20342. The
work of the fifth author was partly supported by JSPS KAKENHI Grant Numbers
JP22K03425, JP22K18677, 23H00086.

## References

  * [1] F. Andreu-Vaillo, V. Caselles, J. M, Mazón, Parabolic quasilinear equations minimizing linear growth functionals. Progress in Mathematics, 223. _Birkhäuser Verlag, Basel,_ 2004\. xiv+340 pp.
  * [2] H. Brézis, Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert. North-Holland Mathematics Studies, No. 5. Notas de Matemática, No. 50. _North-Holland Publishing Co., Amsterdam-London; American Elsevier Publishing Co., Inc., New York,_ 1973\. vi+183 pp.
  * [3] L. Caffarelli and L. Silvestre, An extension problem related to the fractional Laplacian. _Comm. Partial Differential Equations_ 32 (2007), no. 7–9, 1245–1260.
  * [4] M.-H. Giga, Y. Giga and R. Kobayashi, Very singular diffusion equations. _Taniguchi Conference on Mathematics Nara ’98,_ 93–125, Adv. Stud. Pure Math., 31, _Math. Soc. Japan, Tokyo,_ 2001\.
  * [5] M.-H. Giga, Y. Giga and J. Saal, Nonlinear partial differential equations. Asymptotic behavior of solutions and self-similar solutions. Progress in Nonlinear Differential Equations and their Applications, 79. _Birkhäuser Boston, Ltd., Boston, MA,_ 2010\. xviii+294 pp.
  * [6] Y. Giga, H. Mitake and S. Sato, On the equivalence of viscosity solutions and distributional solutions for the time-fractional diffusion equation. _J. Differential Equations_ 316 (2022), 364–386.
  * [7] Y. Giga and T. Namba, Well-posedness of Hamilton–Jacobi equations with Caputo’s time fractional derivative. _Comm. Partial Differential Equations_ 42 (2017), no. 7, 1088–1120.
  * [8] Y. Giga, J. Okamoto, K. Sakakibara and M. Uesaka, On a singular limit of the Kobayashi–Warren–Carter energy. _Indiana Univ. Math. J.,_ to appear.
  * [9] Y. Giga, J. Okamoto and M. Uesaka, A finer singular limit of a single-well Modica–Mortola functional and its applications to the Kobayashi–Warren–Carter energy. _Adv. Calc. Var._ 16 (2023), 163–182.
  * [10] A. Ito, N. Kenmochi and N. Yamazaki, A phase-field model of grain boundary motion. _Appl. Math._ 53 (2008), no. 5, 433–454.
  * [11] R. Kobayashi, J. A. Warren and W. C. Carter, A continuum model of grain boundaries. _Phys. D_ 140 (2000), 141–150.
  * [12] R. Kobayashi, J. A. Warren and W. C. Carter, Grain boundary model and singular diffusivity. _Free boundary problems: theory and applications, II (Chiba, 1999),_ 283–294, GAKUTO Internat. Ser. Math. Sci. Appl. 14, Gakkōtosho, Tokyo, 2000.
  * [13] R. Kobayashi, J. A. Warren and W. C. Carter, Modeling grain boundaries using a phase-field technique. _J. Cryst. Growth_ 211 (2000), 18–20.
  * [14] Y. Kōmura, Nonlinear semi-groups in Hilbert space. _J. Math. Soc. Japan_ 19 (1967), 493–507.
  * [15] A. Kubica, K. Ryszewska and M. Yamamoto, Time-fractional differential equations – a theoretical introduction. SpringerBriefs in Mathematics. _Springer, Singapore,_ 2020\. x+134 pp.
  * [16] A. Lunardi, Analytic semigroups and optimal regularity in parabolic problems. [2013 reprint of the 1995 original] Modern Birkhäuser Classics. _Birkhäuser/Springer Basel AG, Basel,_ 1995\. xviii+424 pp.
  * [17] S. Moll and K. Shirakawa, Existence of solutions to the Kobayashi–Warren–Carter system. _Calc. Var. Partial Differential Equations_ 51 (2014), no. 3–4, 621–656.
  * [18] S. Moll, K. Shirakawa and H. Watanabe, Energy dissipative solutions to the Kobayashi–Warren–Carter system. _Nonlinearity_ 30 (2017), no. 7, 2752–2784.
  * [19] S. Moll, K. Shirakawa and H. Watanabe, Kobayashi–Warren–Carter type systems with nonhomogeneous Dirichlet boundary data for crystalline orientation, in preparation.
  * [20] T. Namba, On existence and uniqueness of viscosity solutions for second order fully nonlinear PDEs with Caputo time fractional derivatives. _NoDEA Nonlinear Differential Equations Appl._ 25 (2018), no. 3, Paper No. 23, 39 pp.
  * [21] I. Podlubny, Fractional differential equations. An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. Mathematics in Science and Engineering, 198. _Academic Press, Inc., San Diego, CA,_ 1999\. xxiv+340 pp.
  * [22] M. H. Protter and H. F. Weinberger, Maximum principles in differential equations. Corrected reprint of the 1967 original. _Springer-Verlag, New York,_ 1984\. x+261 pp.
  * [23] K, Sakamoto and M. Yamamoto, Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems. _J. Math. Anal. Appl._ 382 (2011), no. 1, 426–447.
  * [24] K. Shirakawa and H. Watanabe, Energy-dissipative solution to a one-dimensional phase field model of grain boundary motion. _Discrete Contin. Dyn. Syst. Ser. S_ 7 (2014), no. 1, 139–159.
  * [25] K. Shirakawa, H. Watanabe and N. Yamazaki, Solvability of one-dimensional phase field systems associated with grain boundary motion. _Math. Ann._ 356 (2013), no. 1, 301–330.
  * [26] E. Topp and M. Yangari, Existence and uniqueness for parabolic problems with Caputo time derivative. _J. Differential Equations_ 262 (2017), no. 12, 6018–6046.
  * [27] H. Watanabe and K. Shirakawa, Qualitative properties of a one-dimensional phase-field system associated with grain boundary, in: _Nonlinear Analysis in Interdisciplinary Sciences – Modellings, Theory and Simulations,_ GAKUTO Internat. Ser. Math. Sci. Appl. 36, Gakkōtosho, Tokyo (2013), 301–328.
  * [28] R. Zacher, Weak solutions of abstract evolutionary integro-differential equations in Hilbert spaces. _Funkcial. Ekvac._ 52 (2009), no. 1, 1–18.
  * [29] R. Zacher, Time fractional diffusion equations: solution concepts, regularity, and long-time behavior. _Handbook of fractional calculus with applications. Vol. 2,_ 159–179, _De Gruyter, Berlin,_ 2019\.