1. Introduction
We consider the Kobayashi–Warren–Carter system, introduced in [Reference Kobayashi, Warren and Carter14, Reference Kobayashi, Warren and Carter15, Reference Warren, Kobayashi and Carter29], to model evolutions of structures in a multi-grain problem. It is a kind of phase-field system in a domain Ω in Rn, formally a gradient flow of the energy
Here, $\alpha_0(v)\geq0$ is a given function, typically $\alpha_0(v)=rv^2$, with a constant r > 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
where $\alpha_w\geq0$ and τ > 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_0v^2$, where τ 0 is a positive constant. We consider the gradient flow of $E_\mathrm{KWC}^\varepsilon$ under this metric, and its explicit form is
An explicit form in [Reference Kobayashi, Warren and Carter14] corresponds to the case when $s=\varepsilon$, $\nu=\varepsilon^2$, $\tau_1=\tau\varepsilon$, $F(v)=(v-1)^2$, $\alpha_0(v)=v^2$, $\alpha_w(v)=\tau_0 v^2/\varepsilon$ with $\tau_0 \gt 0$. In other words,
The function v represents an order parameter, where v = 1 corresponds to a grain region and where v 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 Eqs. (1.3) and (1.4) as $\varepsilon\downarrow0$. It turns out that the correct scaling of time should be $\tau=\tau_1/\varepsilon$, while τ 1 is independent of ɛ. Since the system Eqs. (1.3)–(1.4) is regarded as a gradient flow of $E_\mathrm{KWC}^\varepsilon(u,v)$ of Eq. (1.1), we are tempted to expect that the limit flow is the gradient flow of its limit energy $E_\mathrm{KWC}^0$ which was obtained in our papers [Reference Giga, Okamoto, Sakakibara and Uesaka8, Reference Giga, Okamoto and Uesaka9]. 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 Eq. (1.3) because handing Eq. (1.4) 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 Eqs. (1.3)–(1.4), where Ω is an interval $\Omega=(-L,L)$, and impose the Dirichlet boundary condition for u and the Neumann boundary condition for v. More precisely,
while
We set
We expect that the function
with b > 0, solves Eq. (1.4). Since Eq. (1.4) is of total variation flow type, the definition of a solution is not obvious. Fortunately, under a suitable assumption of v, say $v(0,t)\leq v(x,t)$ for all $x\in(-L,L)$, t > 0, the function ub solves Eq. (1.4) under Eq. (1.7), 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
Then, ub solves
under Eq. (1.7).
We stress that the notion of a solution of the equation for u in lemma 1.1 and also Eq. (1.4) under the Dirichlet condition Eq. (1.7) is not obvious and will be discussed in §4.
Problem Eqs. (1.3)–(1.4) are reduced to
under the boundary condition
and the initial condition
where $1_{x \gt 0}$ is a characteristic function of $(0,+\infty)$, i.e., the Heaviside function, so that its distributional derivative equals the Dirac δ function. We note that if we assume that $v_0^\varepsilon$ is even and non-decreasing for x > 0, i.e., $v_0^\varepsilon(x)=v_0^\varepsilon(-x)$ and $v_{0x}^\varepsilon(x)x\geq0$ for $x\in(-L,L)$, then ub is indeed a stationary solution of Eq. (1.4) with Eq. (1.7). This follows from lemma 1.1 since we have $v^\varepsilon(0,t)\leq v^\varepsilon(x,t)$ for all $x\in(-L,L)$, t > 0 by the maximum principle. By a scaling transformation $y=x/\varepsilon$, Eq. (1.10) becomes
in $(-L/\varepsilon,L/\varepsilon)\times(0,\infty)$ for $V=V^\varepsilon(y,t)=v^\varepsilon(\varepsilon y,t)$, where v ɛ is a solution of Eq. (1.10). Thus, we expect that this limit $V=\lim_{\varepsilon\to0}V^\varepsilon$ solves Eq. (1.13) on R and is bounded. Since the solution v ɛ of Eq. (1.10) 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\to0}\xi^\varepsilon$ solves.
We set
Since $F(v)=a^2(v-1)^2$ to get $G(\zeta)=a(\zeta-1)^2/2$, we have
As proved in [Reference Giga, Okamoto and Uesaka9], this energy is obtained as a value at Ξ of Γ-limit of $E_\mathrm{sMM}^\varepsilon(v)$ in Eq. (1.2) as ɛ → 0 if v converges to a set-valued function Ξ (under the graph convergence) of the form
for $\zeta\in(0,1)$.
A key observation is to derive an equation for ξ. For $\alpha\in\mathbf{R}$, β > 0, we set
where Γ denotes the gamma function. We consider well-prepared initial data in the sense that it is continuous and solves Eq. (1.13) outside y = 0.
Lemma 1.2. (Limit equation)
Assume that $\tau_1=1$. Let V be the bounded solution of Eq. (1.13) in $\mathbf{R}\times(0,\infty)$ with initial data $V(y,0)=1-ce^{-a|y|}$ with some $c\in\mathbf{R}$. Then, $\xi(t)=V(0,t)$ solves
with
Moreover, $m_a^\prime(t) \lt 0$ and $m_a(t) \gt 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 τ 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 Eq. (1.14) becomes the Caputo derivative $2\partial_t^{1/2}$. For even initial data, the solution V of Eq. (1.13) is even. Since Vy is odd in y, we see that the delta part of Vyy equals
where $V(\pm0,t)=\lim_{\gamma\to\pm0}V(\gamma,t)$ and δ denotes the Dirac delta function. Since V is smooth outside y = 0 and Vt has no delta part, Eq. (1.13) deduces that
Thus, Eq. (1.13) in $\mathbf{R}\times(0,\infty)$ is reduced to the Robin boundary problem in $(-\infty,0)\times(0,\infty)$ with
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
with the Robin boundary condition
Assume that $w(x,0)=-c$ for some $c\in\mathbf{R}$; the boundary value $\xi(t)=w(0,t)$ solves
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 [Reference Caffarelli and Silvestre3] show that $(-\Delta)^\gamma$ ($0 \lt \gamma \lt 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 [Reference Caffarelli and Silvestre3]; see remark 2.4 at the end of §2. We do not pursue this problem in this paper.
Since Eq. (1.13) is linear and of constant coefficients, we can use the Laplace transform to obtain the desired equation. As one expects, we are able to prove the convergence of V ɛ.
Lemma 1.4. (Convergence)
Let V be the bounded solution of Eq. (1.13) in $\mathbf{R}\times(0,\infty)$ with initial data V 0, which is bounded and uniformly continuous on R. Let vɛ be the solution of Eq. (1.10) under Eqs. (1.11) and (1.12). Assume that $V_0^\varepsilon=\left.V^\varepsilon\right|_{t=0}\in C\left[-L/\varepsilon,L/\varepsilon\right]$ converges to V0 in the sense that
Then, Vɛ converges to V locally uniformly in $\mathbf{R}\times[0,\infty)$. In particular, $\xi^\varepsilon(t)=v^\varepsilon(0,t)=V^\varepsilon(0,t)$ converges to $\xi(t)=V(0,t)$ locally uniformly in $[0,\infty)$. If we assume that
then
for any T > 0. Moreover,
Although the statement is for linear equations, we have to be a little bit careful since the domains where solutions are defined depend on ɛ. We extend the initial data $V_0^\varepsilon$ to a whole space. If the solution W ɛ to the whole space problem with the extended data $\overline{V_0^\varepsilon}$ solves the Neumann problem in $(-L/\varepsilon,L/\varepsilon)$, it is very convenient. We take this strategy. However, $\overline{V_0^\varepsilon}$ may not be close to V 0 outside $(-L/\varepsilon,L/\varepsilon)$, so the maximum principle does not imply the convergence of the solution W ɛ starting from $\overline{V_0^\varepsilon}$ to V even locally. We regularize initial data and prove that $\{W^\varepsilon\}$ is bounded and equi-continuous in $\mathbf{R}\times[0,T]$, $0 \lt T \lt \infty$. By Arzelà–Ascoli theorem and a diagonal argument, it converges to some function by taking a subsequence locally uniformly in $\mathbf{R}\times[0,T]$. The limit is identified with a solution V, so it is a full convergence. This eventually yields that V ɛ converges to V locally uniformly in $\mathbf{R}\times[0,\infty)$, so it yields the convergence of ξ ɛ. Because of the term $V\partial_y(1_{y \gt 0})$, the actual proof is more involved. We divide the initial data into even part and odd part. For the odd part, the problem is reduced to the problem without the $V\partial_y(1_{y \gt 0})$ term. For the even part, the problem is reduced to the Robin boundary problem and the proof is more involved. Note that we do not assume behaviour of $V_0(y)$ so $|y|\to\infty$ to derive locally uniform convergence.
To obtain the uniform convergence of V ɛ to V in $(-L/\varepsilon,L/\varepsilon)\times(0,T)$, it seems that the behaviour of $V_0(y)$ as $|y|\to\infty$ should be controlled. We observe that our decay assumption for $V_0(y)-1$ is preserved not only for V but also for $\partial_yV$ for $t \gt \delta \gt 0$. This allows us to compare V ɛ with V globally in space. All proofs are quite elementary, but we give full proofs for the reader’s convenience.
Applying lemmas 1.2 and 1.4, we can obtain a characterization of the limit equation.
Theorem 1.5 Assume that $\tau_1=1$. Let vɛ be the solution of Eq. (1.10) under Eqs. (1.11) and (1.12). Assume that $v_0^\varepsilon$ is well-prepared in the sense
as $\varepsilon \to 0$, with some c independent of ɛ. Then, $\xi^\varepsilon(t)$ converges to ξ locally uniformly in $[0,\infty)$, and ξ solves Eq. (1.14). Moreover, the graph of vɛ converges to a set-valued function Ξ of the form
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, Eq. (1.14) is altered because there are a few lower-order terms; see Eq. (2.9). If $V(y,0)=1-ce^{-\mu|x|}$ with $\mu\neq a$, we still get an explicit form corresponding to Eq. (1.14); see Eq. (2.16). In these cases, the solution ξ 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 §2. We also calculate numerically how the solution of the ɛ-problem converges by comparing it with the explicit solution of Eq. (1.14) and more general Eq. (2.16).
We now come back to the singular limit problem of the Kobayashi–Warren–Carter system Eqs. (1.3)–(1.4), with $\tau=\tau_1/\varepsilon$, $F(v)=a^2(v-1)^2$, and $\alpha_0(v)=v^2$ as ɛ → 0. We give a purely formal argument. Eq. (1.4) 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) \gt 0$; in this case, Eq. (1.4) 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 the Neumann boundary condition. As in [Reference Giga, Giga and Kobayashi4], we consider Eq. (1.4) for piecewise constant functions. For an interval $I=(p,q)$ and its given division $p=p_0 \lt p_1 \lt \cdots \lt p_m=q$, we consider a piecewise constant function of form
We interpret a solution by mimicking the notion of a solution when v is independent of time and $\inf\alpha_w(v) \gt 0$ (under the periodic boundary condition for simplicity). It is of the form
where $\chi_j=\operatorname{sgn}(h_{j+1}-h_j)$; we identify $p_m=p_0$ by periodicity. See figure 1. The function z is called a Cahn–Hoffman vector field.
For a given v, it is unclear whether the solution remains in a class of piecewise constant functions [Reference Giga, Giga and Kobayashi4], known as a facet-splitting problem. If v is constant, it is well known that the solution remains in a class of piecewise constant functions.
We postulate that our system Eq. (1.15) has a (spatially) piecewise constant solution. Integrating the first equation of Eq. (1.15) in $(p_{j-1},p_j)$ yields
If $v_\varepsilon\to\Xi$ in the graph topology as ɛ → 0, we arrive at
where the set-valued function Ξ is defined as
Since jump $b_j=|h_{j+1}-h_j|$ determines the ξj-equation as indicated in theorem 1.5, the equation for ξj (assuming $\tau_1=1$) is expected to be
where $M_a f=\int_0^t m_a(t-s)f(s)\,ds$ for well-prepared initial data at least when $h_{j+1}-h_j$ is independent of time t. We postulate that Eq. (1.17) is still valid when hj depends on time. Thus, the singular limit equation of the Kobayashi–Warren–Carter equations Eqs. (1.3)–(1.4) as ɛ → 0 (under the periodic boundary condition) is expected to be the system Eqs. (1.16)–(1.17) for hj and ξj. If we consider other boundary conditions, we always impose the homogeneous Neumann boundary condition for v, like Eq. (1.8). If we consider the Dirichlet boundary condition for u with time-independent data, the values of h 1 and hm are prescribed. We can impose the Neumann condition for u; in this case, imposing $z(p_0)=z(p_m)=0$ in Eq. (1.15) is natural. It is rather standard [Reference Podlubny23] 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 Eqs. (1.5)–(1.6) is still an open problem, even in a one-dimensional setting. If we write it in the form of Eqs. (1.3)–(1.4), the difficulty stems from the fact that the weights αw and α 0 can vanish somewhere. In the literature, αw is assumed to be away from zero. If α 0 is allowed to vanish, the $\Delta u$ term is added in the right-hand side of Eq. (1.6). For example, instead of considering Eq. (1.6), we consider
with δ > 0, $\delta^\prime\geq0$ and $\mu\geq0$ such that $\delta^\prime+\mu \gt 0$. The existence of a solution to Eqs. (1.5) and (1.18) with its large time behaviour is established in [Reference Ito, Kenmochi and Yamazaki10–Reference Kenmochi and Yamazaki13, Reference Moll and Shirakawa19, Reference Moll, Shirakawa and Watanabe20, Reference Shirakawa and Watanabe26, Reference Shirakawa, Watanabe and Yamazaki27, Reference Watanabe and Shirakawa30] under several homogeneous boundary conditions. (The case δ = 0 is included in [Reference Ito, Kenmochi and Yamazaki11, Reference Watanabe and Shirakawa30].) Unfortunately, the uniqueness of their solution is only known in a one-dimensional setting under the relaxation term µ > 0 [Reference Ito, Kenmochi and Yamazaki10, Theorem 2.2]. The extension of these results to inhomogeneous boundary condition is not difficult. In [Reference Moll, Shirakawa and Watanabe21], under non-homogeneous Dirichlet boundary conditions, structured patterns of stationary (i.e., time-independent) solutions were studied. In a one-dimensional setting, they thoroughly characterized all stationary solutions. Our Γ-convergence result in [Reference Giga, Okamoto and Uesaka9] gives the convergence of minimizer of $E_\mathrm{KWC}^\varepsilon$ to that of the limit energy as ɛ → 0. We do not know the convergence of stationary solutions.
This paper is organized as follows. In §2, we derive the limit equation for both well-prepared and non-well-prepared initial data. We also give the large time behaviour of a solution. Most of the calculations are very explicit. In §3, we prove lemma 1.4. Section 4 gives a rigorous definition of the Dirichlet problem for Eq. (1.4), assuming that v is given and $\alpha_w(v)\equiv1$. In §5, we give several numerical tests.
2. Derivation of equations with the fractional time derivative
The first goal of this section is to prove lemma 1.2. We then consider more general initial data. We in particular consider the case $V(y,0)=1-ce^{-\mu|y|}$ with µ > 0 not necessarily equal to a and derive the equation corresponding to Eq. (1.14); see Eq. (2.16). The solution ξ can be written explicitly (see Eq. (2.18)), and it is of the form
with a lower-order term
This asymptotic result is shown for general µ; see lemma 2.3 for a detailed statement. At the end of this section, we give remarks related to corollary 1.3. It turns out that corollary 1.3 follows from an expression of the Dirichlet–Neumann map of the heat equation. The Caputo derivatives can be naturally derived using the Dirichlet–Neumann map.
We begin with recalling several elemental properties of the Laplace transform
for a locally integrable function g in $[0,\infty)$. By definition,
and by definition of the Gamma function, we see
We now arrive at a well-known formula
for $\alpha,\beta \gt 0$.
Lemma 2.1. For a > 0, let $m_a(t)$ be
Then
Moreover, $m_a^\prime(t) \lt 0$ and $m_a(t) \gt 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
If we take $g_1\equiv1$, we see
since $\mathcal{L} [g_1](\lambda)=\lambda^{-1}$ by Eq. (2.2). Setting $\alpha=a^2$, we apply Eqs. (2.2) and (2.4) to get
which yields the desired formula for $\mathcal{L}[m_a]$.
We differentiate $m_a/2$ to get
We observe that
Thus, $\lim_{t\to\infty}m_a(t)=0$ since $\lim_{t\to\infty}f_{1/2}^\alpha(t)=0$, which implies that $m_a(t) \gt 0$ for t > 0 since $m_a^\prime(t) \lt 0$.
Proof of lemma 1.2
Since the property of ma has been proved in lemma 2.1, it suffices to derive Eq. (1.14). Studying the equation for $w=V-1$ instead of Eq. (1.13) is more convenient. Eq. (1.13) for w (with $\tau_1=1$) becomes
The initial data $w(x,0)$ equal
with $w_0(x)=-ce^{-a|x|}$. Let $\widehat{w}(x,\lambda)$ be the Laplace transform of w in the t variable, i.e.,
We note that
Taking the Laplace transform of Eq. (2.5), we arrive at
a second-order linear ordinary differential equation with a jump of the derivatives $\widehat{w}_x$. Since the coefficients are constants, we can solve Eq. (2.6) explicitly.
A bounded solution satisfying
is of form
with $A\in\mathbf{R}$. One should determine A so that Eq. (2.6) holds. We observe
so that
since $(\operatorname{sgn}x)^\prime=2\delta$. Eq. (2.6) imposes that
which implies
Since $\widehat{w}(0,\lambda)=A-c\lambda^{-1}$, we end up with
We set $\eta(t)=w(0,t)$ and observe that
We next extract the time derivative of η. We note that
so that $\widehat{\eta_t}=\lambda\widehat{\eta}+c$ since $\eta(0)=-c$. Then,
Multiplying $\sqrt{\lambda+a^2}+b$ and subtracting $(a+b)\widehat{\eta_t}$ from both sides, we get
The energy for ξ equals
Since $\eta=\xi-1$, we see $\operatorname{grad}E_\mathrm{sMM}^{0,b}(\xi)=2(b+a)\eta+2b$. Thus,
By lemma 2.1 and Eq. (2.3), we now conclude that $\xi=\eta+1$ solves Eq. (1.14) by taking the inverse Laplace transform of Eq. (2.7).
In lemma 1.2, we only considered well-prepared initial data, meaning that $V(y,t)$ is a stationary solution of Eq. (1.13) in $\mathbf{R}\backslash\{0\}$. We take this opportunity to write a general equation corresponding to Eq. (1.14) starting from general initial data. We set
and Eq. (1.14) is
For general initial data, our equation corresponding to Eq. (2.8) becomes more complicated than Eq. (2.8).
We set
which is the Green function of $-\partial_y^2+\sigma$ with $\sigma=\lambda+a^2$, i.e.,
For $w_0\in L^\infty(\mathbf{R})$, we set
where $*_x$ denotes the convolution in space, i.e.,
Lemma 2.2. Let $\tau_1=1$. Let V be the bounded solution of Eq. (1.13) in $\mathbf{R}\times(0,\infty)$ with initial data $V(y,0)=1+w_0$, where w 0 is bounded and Lipschitz continuous in R. Then, $\xi(t)=V(0,t)$ solves
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
is integrable for all $x\in\mathbf{R}$ near s = 0; see lemma 3.4.
We argue in the same way to prove lemma 1.2. We begin with Eq. (2.6), where $w=V-1$. Its solution outside x = 0 is given as
with $A\in\mathbf{R}$. As before, we determine A and obtain that $\eta=\xi-1$ satisfies
Since $\widehat{\eta_t}=\lambda\widehat{\eta}-\eta(0)$, we proceed with
Thus,
Taking the inverse Laplace transform, we obtain
the same as Eq. (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.,
We know that its Laplace transform (as a function of t) is
Let $E^a=e^{-a^2 t}E$, then
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
Thus,
Since $\sqrt{\lambda+a^2}g^a=\frac12\left(e^{-\sqrt{\lambda+a^2}|x|} *_x w_0 \right)(0,t)$, we conclude that
If $w_0=-ce^{-\mu|x|}$, we observe that
where $\operatorname{erfc}$ denotes the complementary error function, i.e.,
Indeed, we proceed with
By a direct calculation, we observe that
Thus,
where we invoked Eqs. (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 Eq. (2.13) becomes
Using Eq. (2.1) again, we now obtain Eq. (2.11).
The formula (2.13) gives another representation of ma. Indeed, since
we see
in particular, implies that $m_a(t)\leq 2f_{1/2}^{a^2}(t)$ and
Moreover, by Eq. (2.15), we see
Therefore, the positivity of ma (lemma 2.1) implies $q^a(t) \gt 0$ for t > 0. The formula Eq. (2.9) in lemma 2.2 becomes
if $w_0=-ce^{-\mu|x|}$. If $\mu=a$, this is reduced to Eq. (1.14) (or Eq. (2.8)).
One can give an explicit form of a solution of Eq. (2.16) starting from $\xi(0)=1-c$. We substitute Eq. (2.12) into Eq. (2.10) to get
Since
by Eqs. (2.13) and (2.1), we have
However, the calculation of $q^b*q^\mu$ is quite involved, and it is easier to calculate $\widehat{\eta}$ in Eq. (2.17) more. We proceed
It is easy to see that
As we already observed,
For the first term,
By definition,
We thus conclude that
Thus, $\xi=\eta+1$ is the solution of Eq. (2.16) with $\xi(0)=1-c$.
From this solution formula, we can establish the solution’s large-time behaviour. We set
(i) The function $\bar{\eta}$ is negative, monotonically decreasing, and converging to $-b/(b+a)$ as $t\to\infty$.
(ii) The estimate
\begin{align*} \left| \frac{\eta_e}{c(\mu+b)} \right| \leq e^{-a^2t} \int_t^\infty q^b(s)\,ds \leq e^{-a^2t} \frac1b \end{align*}holds for t > 0. In particular,
\begin{align*} \lim_{t\to\infty}\eta_e e^{a^2t} = 0. \end{align*}
(i) Since $q^b\geq0$, the monotonicity is clear. We observe that
\begin{align*} \int_0^\infty e^{-a^2s}q^b(s)\,ds = \mathcal{L} \left[q^b\right] (0+a^2) = \frac{1}{\sqrt{a^2}+b} \end{align*}by Eq. (2.13). Thus, $\lim_{t\to\infty} \bar{\eta}(t)=-b/(b+a)$.
(ii) Since
\begin{align*} \int_0^\infty e^{-\mu^2s}q^b(s)\,ds = \frac{1}{\mu+b}, \end{align*}we observe that
\begin{align*} \eta_e &= -c(\mu+b) \int_t^\infty e^{-\mu^2s} q^b(s)\,ds\,e^{-(a^2-\mu^2)t} \\ &= -c(\mu+b) \int_t^\infty e^{\mu^2(t-s)} q^b(s)\,ds\,e^{-a^2t}. \end{align*}Since $e^{\mu^2(t-s)} \leq 1$ for $s\geq t$, this implies $\left| \eta_e/c(\mu+b)\right| \leq e^{-a^2t} \int_t^\infty q^b(s)\,ds \leq e^{-a^2t} \int_0^\infty q^b(s)\,ds = e^{-a^2t}/b$. The proof is now complete.
We conclude this section by giving remarks on corollary 1.3. It turns out that corollary 1.3 can be derived using the expression of the Dirichlet–Neumann map.
Remark 2.4. We consider the initial-boundary value problem for
Then, as in [Reference Caffarelli and Silvestre3], we obtain that
where $\gamma=1/(2-\alpha)$ with some constant $c_\gamma \gt 0$, provided that α < 1. Indeed, as in [Reference Caffarelli and Silvestre3], let ψ be a solution of
Since the Laplace transform $\widehat{w}$ of w satisfies
we see, by scaling, that
Thus,
since $\eta(0)=0$. Thus,
If γ < 1, then $f_{1-\gamma}^0$ is integrable. As noted in [Reference Caffarelli and Silvestre3], $\psi^\prime(0)$ exists (even for the degenerate case, i.e., α > 0) and $\psi^\prime(0) \gt 0$. Thus,
at least for $\gamma=1/(2-\alpha) \lt 1$, i.e., α < 1.
Corollary 1.3 is easily derived by this result since $w_x=\partial_t^{1/2}w$ and $w_x+bw=0$; note that $c_{1/2}=1$.
Remark 2.5. The reader might be interested in how fractional partial differential equations like fractional diffusion equations are derived. We consider
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
This type of equation is a kind of fractional diffusion that has been well-studied; see [Reference Kubica, Ryszewska and Yamamoto17, Reference Zacher32]. Here, we briefly recall only the well-posedness of its initial boundary value problem for Eq. (2.19) in a domain. In the framework of distributions, the well-posedness of its initial boundary value problems has been established in [Reference Sakamoto and Yamamoto25, Reference Zacher31] by using the Galerkin method. The unique existence of viscosity solutions for Eq. (2.19), including general nonlinear problems, has been established in [Reference Giga and Namba7, Reference Namba22] and also in [Reference Topp and Yangari28] for the whole space $\mathbf{R}^{n-1}$. The scope of equations these theories apply is different. However, it has been proved in [Reference Giga, Mitake and Sato6] that two notions of solutions (viscosity solution and distributional solution) agree for Eq. (2.19) when we consider the Dirichlet problem in a smooth bounded domain.
3. Convergence
The goal of this section is to prove lemma 1.4. For this purpose, we shall study a homogeneous version of Eq. (1.13) in a whole space $\mathbf{R}\times(0,\infty)$ and estimate the derivative of a solution as $|y|\to\infty$. We begin with several estimates related to the heat semigroup in R. Let $E(x,t)$ be the Gauss kernel, and write $E_t(x)=E(x,t)$, i.e.,
and $E_t^a=e^{-a^2t}E_t$. Then,
solves
with initial data $f\in L^\infty(\mathbf{R})\cap C(\mathbf{R})$, i.e., f is bounded and continuous on R.
Proposition 3.1. Assume that $f\in C(\mathbf{R})\cap L^\infty(\mathbf{R})$ satisfies a decay condition
Then, for $T \gt \delta \gt 0$,
and
Proof. We may assume a = 0. We notice that
with some C independent of t and x, since $\partial_x E_t=-(x/2t)E_t$ and $\sup_{y \gt 0}ye^{-y^2} \lt \infty$. Thus, it suffices to prove that
We divide
We notice that
with some constant C ʹ independent of $x,y\in\mathbf{R}$ and t < T. Since $2|y|\leq|x|$ implies $|x-y|\geq|x|-|y|\geq|x|/2$, we obtain that
For $II$, we proceed
We thus obtain Eq. (3.1). The proof is now complete.
We also need to estimate decay as $|x|\to\infty$ for
when $g(s)=h(s)\partial_x(1_{x \gt 0})$, $h\in L^\infty(0,T)$. This quantity equals
Proposition 3.2. For $h\in L^\infty(0,T)$ and m > 0,
Proof. Again, we may assume that a = 0. Since
for $|x|\geq M$ and $t\leq T$. Since $t^{-1}e^{-M^2/(16t)}$ is integrable near t = 0, we see that
with some constants C 0, $C^\prime_0$ depending only on M and T. The proof is now complete.
We consider a homogeneous version Eq. (1.13) in $\mathbf{R}\times(0,\infty)$. Namely, we consider
Proposition 3.3. Let w be a bounded solution of Eq. (3.2) with initial data $w_0\in C(\mathbf{R})\cap L^\infty(\mathbf{R})$. Assume that
Then,
for any $\delta,T$ satisfying $T \gt \delta \gt 0$.
Proof. We may assume $\tau_1=1$. By Duhamel’s formula, w is of the form
Applying propositions 3.1 and 3.2, we obtain desired results.
Proof. Proof of lemma 1.4
For a function f on $(-L,L)$, we decompose it into its odd and even parts, i.e.,
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 Eq. (1.13) separately.
At first glance, the locally uniform convergence follows from the maximum or comparison principles for a linear parabolic equation [Reference Protter and Weinberger24]. However, a direct application of the maximum principle is impossible since the domains of functions V ɛ and V are different. We first show the convergence where initial data are smooth.
For the odd part, the term $V\partial_y(1_{y \gt 0})$ does not affect since V = 0 at y = 0. Thus, Eq. (1.13) is reduced to
where $V_0^\varepsilon(y)=V^\varepsilon (y,0)$ for $|y| \lt L/\varepsilon$. Let V ɛ be its solution.
We extend an odd function $V_0^\varepsilon$ to $\left(L/\varepsilon,3L/\varepsilon\right)$ for x to be ‘even’ with respect to $L/\varepsilon$, i.e.,
where $\widetilde{V_0^\varepsilon}$ is its extension. We extend $\widetilde{V_0^\varepsilon}$ outside $\left(-L/\varepsilon,3L/\varepsilon\right)$ so that the extension $\overline{V_0^\varepsilon}$ is periodic in 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 ɛ is the restriction on $\left(-L/\varepsilon,L/\varepsilon\right)$ of a solution W ɛ of
Although the maximum principle implies
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})}\to0$. We argue differently.
We approximate V 0 by $V_{0\delta}=V_0*\rho_\delta$, where ρδ is a symmetric mollifier. We also approximate $V_0^\varepsilon$ by
We set $W_\delta^\varepsilon=W^\varepsilon*\rho_\delta$, where W ɛ is the solution of Eq. (3.4). Since Eq. (3.4) is of constant coefficients, this $W_\delta^\varepsilon$ solves Eq. (3.4)1 with initial data $V_{0\delta}^\varepsilon$. Let $V_\delta^\varepsilon$ be the restriction of $W_\delta^\varepsilon$ on $(-L/\varepsilon,L/\varepsilon)$. It follows from the parity and periodic condition that this $V_\delta^\varepsilon$ solves Eq. (3.3). We set $V_\delta=V*\rho_\delta$ and observe that Vδ is the bounded solution of Eq. (3.4)1 with initial data $V_{0\delta}$ since Eq. (3.4) is of constant coefficients. For fixed δ > 0, we observe that
for M > 0. Indeed, by the maximum principle
where $\|\cdot\|_{\infty,\varepsilon}$ is the sup norm on $\left(-L/\varepsilon, L/\varepsilon\right)$. Here, the notation $\partial_tV_{0\delta}^\varepsilon$ for a function $V_{0\delta}^\varepsilon$ of y should be interpreted as
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 for ɛ < 1. By the Arzelà–Ascoli theorem and a diagonal argument, $V_\delta^\varepsilon$ converges (locally uniformly in $\mathbf{R}\times[0,\infty)$) to a bounded (weak) solution to Eq. (3.4) with initial data $V_{0\delta}$ 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 the uniform bound on derivatives.
We note that
and observe that
where the norm is $\|\cdot\|$ taken in $L^\infty(-M,M)$ for M > 0. By the maximum principle,
Since
and $\|\rho_\delta*f\|_\infty\leq\|f\|_\infty$, we see that
Thus,
Fixing δ > 0 and sending ɛ → 0, we observe that
since $V_\delta^\varepsilon\to V_\delta$ in $L^\infty\left((-M,M)\times[0,T]\right)$. Sending $\delta\downarrow0$, we obtain
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
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\to0}V_0^\varepsilon$ by a smoother function uniformly. There are many possible ways, and we rather like an abstract way. Let $\text{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
If ɛ = 0, then Iɛ should be interpreted as $(-\infty,0]$. Let A be the operator on $\text{BUC}(I_0)$ defined by
with
A standard theory [Reference Lunardi18] implies that −A generates an analytic semigroup e −tA in $\text{BUC}(I_0)$. In particular,
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 \lt -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$ although 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 \gt 0$ so that $\sigma^\prime\left(1/2\right)=-\kappa$ with κ < 0. We set $\sigma_{\delta^{\prime\prime}}(y)=\delta^{\prime\prime}\sigma\left(y/\delta^{\prime\prime}\right)$ for small $\delta^{\prime\prime} \gt 0$. For a given $h\in C^2(-\infty,0]$, we modify
where we take c so that $c\kappa=h_y\left(-L/\varepsilon\right)$, $-L/\varepsilon-\nu=\delta^{\prime\prime}/2$. By this modification,
and $h_{\delta^{\prime\prime}}\to h$ in $L^\infty(-\infty,0)$ as $\delta^{\prime\prime}\to0$. We set
Since $V_{0\delta}^\varepsilon$ satisfies the boundary condition on the boundary of Iɛ 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 Eq. (3.5) with initial data $V_{0\delta}^\varepsilon$. By the maximum principle,
where $\partial_tV_{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 Eq. (3.4), by the Arzelà–Ascoli theorem and the uniqueness of the limit equation, we can prove that $V_\delta^\varepsilon$ converges to Vδ locally uniformly in $(-\infty,0]\times[0,\infty)$. Note that for a fixed δ > 0, the right-hand sides of Eq. (3.6) are uniformly bounded as ɛ → 0 since $V_0^\varepsilon$ converges to V 0 uniformly. The comparison principle implies that
Thus,
where the norm $\|\cdot\|$ is taken in $L^\infty(0,M)$ for M > 0. Taking the supremum in $t\in(0,T)$ and sending ɛ → 0, we obtain that
since we know $\sup_{0 \lt t \lt T}\|V_\delta-V_\delta^\varepsilon\|(t)\to0$ as ɛ → 0. Sending δ → 0, we conclude that V ɛ converges to V locally uniformly in $[0,\infty)\times[0,\infty)$.
Since we know that V ɛ is continuous up to y = 0 and t = 0, this gives the local uniform convergence of $V^\varepsilon(0,t)$ in $[0,\infty)$.
So far, we do not invoke the assumption that $\left(V_0(x)-1\right)|x|\to0$ as $|x|\to\infty$. We shall prove
We only discuss the case of even initial data since the odd case is easier. Since
it suffices to prove that
By our construction, the function $V_\delta(y,t)$ is
for w solving Eq. (3.2) with $w_0=V_0-1$. By proposition 3.3,
The difference $V_\delta-V_\delta^\varepsilon=u^\varepsilon$ satisfies
Let $\bar{u}^\varepsilon$ be the solution of
By the maximum principle,
Since we know that $\sup_{0 \lt t \lt T}\left|u_b^\varepsilon(x)\right|\to0$ as ɛ → 0, sending ɛ → 0 yields
We set $r=u^\varepsilon-\bar{u}^\varepsilon$ and observe that r satisfies
We set $\tilde{r}=e^{a^2t}r$ and observe that $\tilde{r}$ satisfies
We take $M^\varepsilon=\sup_{0\leq t\leq T}e^{a^2t}\left|(\partial_y V_\delta)(-L/\varepsilon,0)\right|$ and observe that −yM and yM are the super- and subsolution, respectively. By the comparison principle,
This implies that
By Eq. (3.8), $M^\varepsilon L/\varepsilon\to0$ as ɛ → 0. We thus conclude that
Sending $\delta\downarrow0$, we obtain Eq. (3.7). The statement that $\sup_{0\leq t\leq T}|V-1|(y)|y|\to0$ as $|y|\to\infty$ follows from proposition 3.3. The proof is now complete.
If we only assume that the initial data for Eqs. (3.4) and (3.5) are bounded and Lipschitz, we have a similar estimate in Eq. (3.6) 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 ɛ = 0 for the reader’s convenience.
Lemma 3.4. Let V be the bounded solution of Eq. (1.13) in $\mathbf{R}\times(0,\infty)$ with the 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
Proof. We give direct proof. We may assume that $\tau_1=1$. We set $w=V-1$ and $u=e^{a^2t}w$ to get
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^2t)$ 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
where
Since we can approximate a smooth w 0, establishing
with some positive constants C and T independent of w 0 suffices, assuming that $\partial_tu$ exists and is bounded in $\mathbf{R}\times(0,T)$ for small T. By the maximum principle Eq. (3.6) and the corresponding estimate for the odd part, we know that
with c independent of w 0 and t > 0. We estimate $\partial_tu$ in Eq. (3.9). Since $\|f\ast_xg\|_{L^\infty(\mathbf{R})}\le\|f\|_{L^1(\mathbf{R})}\|g\|_{L^\infty(\mathbf{R})}$ and
we easily see (cf. [Reference Giga, Giga and Saal5, Chapter 1]) that
with c ʹ independent of w 0. The second term of the right-hand side of Eq. (3.9) is more involved than the first term because h contains u. We observe that
Since $|E|\le(4\pi t)^{-1/2}$, it holds that
for $t\in(0,T)$.
Thus,
with some constants Cj ($j=1,2,3$). By Eqs. (3.9) and (3.12), we now observe that
with C 4 independent of w 0, u, and T. Applying estimate for $\|u\|_\infty$ in Eq. (3.11), we conclude Eq. (3.10) for sufficiently small T by absorbing $C_2T^{1/2}\sup_{0 \lt t \lt T}\|t^{1/2}\partial_tu\|_{L^\infty(\mathbf{R})}$ on 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
where $\beta\in C\left(I\times[0,T]\right)$ is a given non-negative function; here, $I=(p_0,p_1)$ is an open interval and T > 0. If we impose the Dirichlet boundary condition
Eq. (4.1) with Eq. (4.2) should be interpreted as an L 2-gradient flow of a time-dependent total variation type energy
when $\int\beta|u_x|$ is a weighted total variation of β 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 β is spatially constant, it is well known that $\Phi^t$ is also lower semicontinuous; see e.g. [Reference Andreu-Vaillo, Caselles and Mazón1]. The solution of Eq. (4.1) with Eq. (4.2) should be interpreted as the gradient flow of form
where $\partial\Phi^t$ denotes the subdifferential of $\Phi^t$ in $L^2(I)$, i.e.,
It is standard that Eq. (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, [Reference Brézis2, Reference Kōmura16]). It applies to the total variation flow case when β is a constant. In one-dimensional case, $u\in BV(I)$ implies $u\in L^\infty(I)$, so the subdifferential becomes
when $\Phi=\Phi^t$; see [Reference Andreu-Vaillo, Caselles and Mazón1, Proposition 5.10 and Lemma 5.13]Footnote 1. Here, $(\beta z,u_x)$ denotes the Anzellotti pair [Reference Andreu-Vaillo, Caselles and Mazón1] and
Eq. (4.3) is
with $|z|\leq1$ in I and
with $-(\beta z)(p_i)(-1)^i(g-\gamma u)(p_i)=\beta|\gamma u-g|(p_i)$ for $i=0,1$. We mimic this notion of the solution. A function $u\in C\left([0,T),L^2(I)\right)$ is a solution to Eq. (4.1) with Eq. (4.2) if there is $z\in L^\infty\left(I\times(0,T)\right)$ such that
Under this preparation, we shall prove lemma 1.1.
Proof. Proof of lemma 1.1
We set $p_0=-L$, $p_1=L$ so that $I=(-L,L)$. Since $u_t^b=0$, Eq. (4.4) says that $\beta z$ is a constant c. The condition Eq. (4.5) is equivalent to saying that $|c|\leq\min\beta=\beta(0)$. Since
with $g(L)=b$, $g(-L)=0$, Eq. (4.6) is equivalent to
In other words, c must be $\beta(0)$. Thus, the existence of z satisfying $|z|\leq1$ is guaranteed if and only if
Eq. (4.6) is fulfilled with $u=u^b$ by taking $g(-L)=0$ and $g(L)=b$. Thus, ub is a stationary solution to Eq. (1.4) with Eq. (1.7).
5. Numerical experiment
In this section, we calculate the solution of Eqs. (1.10)–(1.12) with $v_0(x)=1-ce^{-a|x|/\varepsilon}$ and compare its value at x = 0 with an explicit solution of Eq. (2.16) whose explicit form is given in Eq. (2.18).
5.1. Numerical scheme
Since the initial function, v 0, is an even function, the original problem Eqs. (1.10)–(1.12) is reduced to
The computational region $[0,L]$ is divided into uniform mesh partitions:
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 vi. The central finite difference approximates the Laplace operator, and the time derivative is approximated by the backward difference, yielding the following linear system:
where $\hat{v}_i$ is the value known at the current time and vi is the value to be found at the next time. The Neumann boundary conditions at x = 0 and x = L are approximated as
by the central finite differences.
5.2. Results
Some results are shown for different values of c for parameters
The results of the numerical experiments are summarized in figure 2. In (a), (b), and (c), the horizontal axis represents time, and the vertical axis represents the value at the origin. As ɛ is decreased, the numerical solution converges to the exact solution to the extent that the exact and numerical solutions overlap. Indeed, the table of $L^\infty$-errors for different values of ɛ and c is shown in (d). The errors for $\varepsilon=2^{-3}$ are of order 10−5, indicating that the solution for small ɛ is an excellent approximation to the solution of the fractional time differential equation obtained as the singular limit.
Acknowledgements
The authors are grateful to the anonymous reviewer for his/her valuable expository suggestions. The work of the first author was partly supported by JSPS KAKENHI Grant Numbers JP24K00531, JP19H00639, and JP20K20342, and by Arithmer Inc., Daikin Industries, Ltd. and Ebara Corporation 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, and 23H00086.