#  Effective Boundary Conditions for the Fisher-KPP Equation on a Domain with
3-dimensional Optimally Aligned Coating  00footnotetext: E-mail addresses:
<EMAIL_ADDRESS>

Xingri Genga,b
a Department of Mathematics, Southern University of Science and Technology,
Shenzhen, P.R. China
b Department of Mathematics, National University of Singapore, Singapore

###### Abstract

We consider the Fisher-KPP equation on a three-dimensional domain surrounded
by a thin layer whose diffusion rates are drastically different from that in
the bulk. The bulk is isotropic, while the layer is considered to be
anisotropic and “optimally aligned”, where the normal direction is always an
eigenvector of the diffusion tensor. To see the effect of the layer, we derive
effective boundary conditions (EBCs) by the limiting solution of the Fisher-
KPP equation as the thickness of the layer shrinks to zero. These EBCs contain
some exotic boundary conditions including the Dirichlet-to-Neumann mapping and
the Fractional Laplacian. Moreover, we emphasize that each EBC keeps effective
indefinitely, even as time approaches infinity.

Keywords. thin layer, asymptotic behavior, eigenvalue problem, effective
boundary conditions, lifespan.

AMS subject classifications. 35K05, 35B40, 35B45,74K35.

##  1 Introduction

Motivated by the scenario of a nature reserve surrounded by a road/buffer
zone. The diffusion rate is isotropic in the nature reserve, while it is
anisotropic and drastically different in the road/buffer zone. Moreover, the
road/buffer zone is thin compared to the scale of the nature reserve,
resulting in multi-scales in the spatial variable. The difference in diffusion
rate and spatial size leads to computational difficulty. To deal with such
situations, we view the road/buffer zone as a thickless surface as its
thickness shrinks to zero, on which “effective boundary conditions” (EBCs) are
imposed. These EBCs not only provide an alternative way for numerical
computation but also give us an analytic interpretation of the effects of the
road/buffer zone.

In this paper, we consider the Fisher-KPP equation on a 3-dimensional domain
surrounded by a thin coating (see Figure 1), which is anisotropic and
“optimally aligned” (a notion proposed in [19]) in the sense that any vector
normal to the interface is always an eigenvector of the diffusion tensor. When
the tangent diffusion rates in the coating are not larger than $o(1/\delta)$,
EBCs and their lifespan were studied by Li, Wang, and Wu [12]. This paper
treats the case when the coating is “optimally aligned” with arbitrary
diffusion rates in the tangent direction, leading to the emergence of exotic
EBCs. Moreover, we also show that the maximal time interval (called lifespan)
of each EBC during which the EBC remains effective is infinite.

$\Omega_{2}$$\partial\Omega$$\Omega_{1}$Figure 1:
$\Omega=\overline{\Omega}_{1}\cup\Omega_{2}$. $\Omega_{2}$ is a thin layer
with thickness $\delta$.$\Gamma$p$\delta$n

Let us introduce our mathematical model as follows: let the domain
$\Omega_{1}$ be surrounded by the coating $\Omega_{2}$ with uniform thickness
$\delta>0$; let the domain
$\Omega=\overline{\Omega}_{1}\cup\Omega_{2}\subset\mathbb{R}^{3}$ as shown in
Figure 1. For any finite $T>0$, consider the initial boundary value problem
with the Dirichlet boundary condition

$\left\\{\begin{array}[]{llr}u_{t}-\nabla\cdot(A(x)\nabla
u)=f(u),&\mbox{$(x,t)\in Q_{T},$}\\\ u=0,&\mbox{$(x,t)\in S_{T},$}&\\\
u=u_{0},&\mbox{$(x,t)\in\Omega\times\\{0\\},$}\end{array}\right.$ (1.1)

where $Q_{T}:=\Omega\times(0,T)$, $S_{T}:=\partial\Omega\times(0,T)$, $0\leq
u_{0}\in L^{\infty}(\Omega)$, and $f(u)=u(1-u)$. Assume that the diffusion
tensor $A(x)$ is given by

$A(x)=\left\\{\begin{aligned} &kI_{3\times 3},&x\in\Omega_{1},\\\
&\left(a_{ij}(x)\right)_{3\times 3},&x\in\Omega_{2},\end{aligned}\right.$
(1.2)

where $k$ is a positive constant independent of $\delta>0$ and the positive-
definite matrix $(a_{ij}(x))$ is anisotropic and “optimally aligned” in the
coating $\Omega_{2}$, and further satisfies

$A(x)\textbf{n}(p)=\sigma\textbf{n}(p)\;\text{ and
}\;A(x)\textbf{s}(p)=\mu\textbf{s}(p),\quad\forall x\in\Omega_{2},$ (1.3)

where $\Omega_{2}$ is thin enough and $\Gamma$ is smooth enough such that the
projection $p$ of $x$ onto $\Gamma$ is unique; $\textbf{n}(p)$ is the unit
outer normal vector of $\Gamma$ at $p$; $\textbf{s}(p)$ is an arbitrary unit
tangent vector of $\Gamma$ at $p$; $\sigma$ and $\mu$ are called the normal
conductivity and the tangent conductivity, respectively.

One main purpose of this paper is to address EBCs on $\Gamma\times(0,T)$ as
the thickness of the layer shrinks to zero.

######  Theorem 1.1.

Suppose that $A(x)$ is given in (1.2) and (1.3). Let $0\leq u_{0}\in
L^{\infty}(\Omega)$ and $f(u)=u(1-u)$. Assume further that $\sigma$ and $\mu$
satisfy

$\lim_{\delta\to 0}\sigma\mu=\gamma\in[0,\infty],\quad\lim_{\delta\to
0}\frac{\sigma}{\delta}=\alpha\in[0,\infty].$

Let $u$ be the weak solution of (1.1), then as $\delta\to 0$, $u\to v$ weakly
in $W_{2}^{1,0}(\Omega_{1}\times(0,T))$ and strongly in
$C([0,T];L^{2}(\Omega_{1}))$, where $v$ is the weak solution of

$\left\\{\begin{array}[]{ll}v_{t}-k\Delta
v=f(v),&(x,t)\in\Omega_{1}\times(0,T),\\\
v=u_{0},&(x,t)\in\Omega_{1}\times\\{0\\},\end{array}\right.$ (1.4)

subject to the effective boundary conditions on $\Gamma\times(0,T)$ listed in
Table 1.

Table 1: Effective boundary conditions on $\Gamma\times(0,T)$ for the Dirichlet problem (1.1). As $\delta\to 0$ | $\frac{\sigma}{\delta}\to 0$ | $\frac{\sigma}{\delta}\to\alpha\in(0,\infty)$ | $\frac{\sigma}{\delta}\to\infty$
---|---|---|---
$\sigma\mu\to 0$ | $\frac{\partial v}{\partial\textbf{n}}=0$ | $k\frac{\partial v}{\partial\textbf{n}}=-\alpha v$ | $v=0$
$\sqrt{\sigma\mu}\to\gamma\in(0,\infty)$ | $k\frac{\partial v}{\partial\textbf{n}}=\gamma\mathcal{J}^{\infty}[v]$ | $k\frac{\partial v}{\partial\textbf{n}}=\gamma\mathcal{J}^{\gamma/\alpha}[v]$ | $v=0$
$\sigma\mu\to\infty$ |  $\nabla_{\Gamma}v=0$, $\int_{\Gamma}\frac{\partial v}{\partial\textbf{n}}=0$ |  $\nabla_{\Gamma}v=0$, $\int_{\Gamma}(k\frac{\partial v}{\partial\textbf{n}}+\alpha v)=0$ | $v=0$

The boundary condition $\nabla_{\Gamma}v=0$ on $\Gamma\times(0,T)$ implies
that $v$ is a constant in $x$, where $\nabla_{\Gamma}$ represents the surface
gradient on $\Gamma$; $k\frac{\partial
v}{\partial\textbf{n}}=\beta\Delta_{\Gamma}v$ can be interpreted as a second-
order partial differential equation on $\Gamma$, where the operator
$\Delta_{\Gamma}$ is the Laplacian-Beltrami defined on $\Gamma$, signifying
that the thermal flux across $\Gamma$ in the outer normal direction results in
heat accumulation that then diffuses with a rate $\beta$; $\mathcal{J}^{h}$ is
the Dirichlet-to-Neumann mapping as defined in the following: for
$h\in(0,\infty)$, and smooth $g$ defined on $\Gamma$, $\mathcal{J}^{h}[g]$ is
defined by

$\mathcal{J}^{h}[g](s):=\Psi_{R}(s,0),$

where $\Psi$ is the bounded solution of

$\left\\{\begin{array}[]{ll}\Psi_{RR}+\Delta_{\Gamma}\Psi=0,&\Gamma\times(0,H),\\\
\Psi(s,0)=g(s),&\Psi(s,H)=0.\end{array}\right.$

The analytic formula for $\mathcal{J}^{h}[g]$ is given in the future, from
which
$\mathcal{J}^{\infty}[g]=\underset{h\to\infty}{\lim}\mathcal{J}^{h}[g]=-\left(-\Delta_{\Gamma}\right)^{1/2}g$
is the fractional Laplacian-Beltrami of order $1/2$.

It is worth mentioning that the convergence of $u$ is valid only in the finite
time interval $[0,T]$. Our interest lies in determining the maximal time
interval on which the uniform convergence persists, referred to as the
lifespan of the EBC. The lifespan represents the duration during which the EBC
remains effective. Notably, the lifespan of an EBC is not necessarily
infinite. In the case of $f(u)\equiv 0$ in (1.1), it turns out that $u\to 0$
as $t\to\infty$, while $v\to\overline{u}_{0}$ when the EBC is Neumann, where
$\overline{u}_{0}$ is the average of $u_{0}$ over $\Omega_{1}$. Thus, the
lifespan of the Neumann EBC is not infinite if $\overline{u}_{0}\neq 0$. In
particular, Pond [18] shows that if $\sigma=\mu=\delta^{m}$ with $m>2$, then
the lifespan of the Neumann EBC is at the order of $o(\delta^{1-m})$ rather
than $O(\delta^{1-m})$.

The remaining part of this paper is to prove the following theorem in which
the lifespan of each EBC in Table 1 is infinite.

######  Theorem 1.2.

Suppose all conditions in Theorem 1.1 hold, then as $\delta\to 0$, the weak
solution $u$ of (1.1) satisfies

$\underset{0\leq
t\leq\infty}{\max}||u(\cdot,t)-v(\cdot,t)||_{L^{2}(\Omega_{1})}\longrightarrow
0,$

where $v$ is the weak solution of (1.4) subject to any effective boundary
condition in Table 1.

The proof of this theorem hinges on a comprehensive understanding of the non-
negative steady states of (1.1) and (1.4). For general $f\neq 0$, the proof
becomes complicated due to the impact of $f$ on altering the steady states of
both $u-$problem and $v-$problem. However, in our specific case where
$f=f(u)=u(1-u)$, this non-linearity facilitates the proximity between the two
heat flows ($u-$flow and $v-$flow), even for large time.

There are lots of profound and intriguing findings in the literature regarding
the study of EBCs. The original derivation of EBCs can be tracked back to the
classic book by Carslaw and Jaeger [2]. Subsequently, the elliptic problem
involving reinforcement/coating was rigorously investigated by Brezis,
Caffarelli, and Friedman [1], followed by Li and Zhang [13, 16] for further
developments. Since then, numerous studies have emerged. See ([4, 5, 6, 14,
12]) for the optimally aligned coating, [8] for the model of wolf movement in
which the diffusion tensor is degenerate along a road, and [10, 9] for the
error estimates between the solution of the original problem and the one of
the effective problem. In the context of ecological problem in our case, see
[15] for the derivation of bulk-surface model, [11] for the road-field model,
and [3] for the effect of an EBC. We also refer interested readers to a review
paper by Wang [20] on the study of EBCs for diffusion equations.

The organization of this paper is as follows. In Section 2, some preliminaries
are presented, including the definition of weak solutions of (1.1) and (1.4),
and basic energy estimates from standard parabolic $L^{2}-$ theory. Section 3
is devoted to deriving EBCs on $\Gamma\times(0,T)$ for the Fisher-KPP
equation, in which an auxiliary function is proposed. In Section 4, based on
the estimates for eigenvalue problems, we investigate the lifespan of EBCs by
proving Theorem 1.2.

Throughout this article, we always assume that $\Omega_{1}$ is fixed and
bounded with $C^{2}$ smooth boundary $\Gamma$; the coating $\Omega_{2}$ is
uniformly thick with $\partial\Omega$ approaching $\Gamma$ as $\delta\to 0$;
$\sigma$ and $\mu$ are positive functions of $\delta$.

##  2 Preliminaries

In the sequel, for notational convenience, let $C(T)$ represent a generic
positive constant depending only on $T$ and $O(1)$ represent a quantity
independent of $\delta$ but may vary from line to line.

###  2.1 Eneregy estimates

Let $W^{1,0}_{2}(Q_{T})$ be the subspace of functions belonging to
$L^{2}(Q_{T})$ with first order weak derivatives in $x$ also being in
$L^{2}(Q_{T})$; $W^{1,1}_{2}(Q_{T})$ is defined similarly with the first order
weak derivative in $t$ belonging to $L^{2}(Q_{T})$; $W^{1,0}_{2,0}(Q_{T})$ is
the closure in $W^{1,0}_{2}(Q_{T})$ of $C^{\infty}$ functions vanishing near
$\overline{S}_{T}$, and $W^{1,1}_{2,0}(Q_{T})$ is defined similarly.
Furthermore, denote $V^{1,0}_{2,0}(Q_{T}):=W^{1,0}_{2,0}(Q_{T})\cap
C\left([0,T];L^{2}(\Omega)\right).$

Let us define one more Sobolev space on
$Q^{1}_{T}=\Omega_{1}\times(0,T):V^{1,0}_{2}(Q^{1}_{T})=W^{1,0}_{2}(Q^{1}_{T})\cap
C\left([0,T];L^{2}(\Omega_{1})\right).$ We endow all these spaces with natural
norms.

Rewrite $\int_{0}^{T}\int_{\Omega}u(x,t)dxdt$ as $\int_{Q_{T}}u(x,t)dxdt$ for
simplicity. We begin with some energy estimates.

######  Definition 2.1.

A function $u$ is said to be a weak solution of the Dirichlet problem (1.1),
if $u\in V^{1,0}_{2,0}(Q_{T})$ and for any $\xi\in
C^{\infty}(\overline{Q}_{T})$ satisfying $\xi=0$ at $t=T$ and also near
$S_{T}$, it holds that

$\mathcal{A}[u,\xi]:=-\int_{\Omega}u_{0}\xi(x,0)dx+\int_{Q_{T}}\left(A(x)\nabla
u\cdot\nabla\xi-u\xi_{t}-f(u)\xi\right)dxdt=0.$ (2.1)

Note that $0$ and $M=\max\\{1,||u_{0}||_{L^{\infty}(\Omega)}\\}$ are the
corresponding lower and upper solution of (1.1). By the method of lower and
upper solutions, for any small $\delta>0$, (1.1) admits a unique weak solution
$u\in V^{1,0}_{2,0}(Q_{T})$, and $0\leq u\leq M$ on $Q_{T}$. Moreover, $u$
satisfies the following “transmission conditions” in the weak sense

$u_{1}=u_{2},\quad k\nabla u_{1}\cdot\textbf{n}=\sigma\nabla
u_{2}\cdot\textbf{n}\quad\text{ on }\Gamma,$ (2.2)

where $u_{1}$ and $u_{2}$ are the restrictions of $u$ on
$\Omega_{1}\times(0,T)$ and $\Omega_{2}\times(0,T)$, respectively.

######  Lemma 2.1.

Suppose $0\leq u_{0}\in L^{\infty}(\Omega)$. Then, any weak solution $u$ of
(1.1) satisfies the following inequalities.

$\begin{split}\text{(i)}&\;0\leq
u(x,t)\leq\max\\{1,||u_{0}||_{L^{\infty}(\Omega)}\\}.\\\
\text{(ii)}&\max_{t\in[0,T]}t\int_{\Omega}\nabla u\cdot A(x)\nabla
udx+\int_{Q_{T}}tu_{t}^{2}dxdt+\int_{Q_{T}}\nabla u\cdot A(x)\nabla udxdt\leq
C\left(T,|\Omega|,||u_{0}||_{L^{\infty}(\Omega)}\right).\\\ \end{split}$

###### Proof.

$(i)$ can be proved by the lower and upper solutions method [17] by choosing
the lower solution 0 and the upper solution
$u=\max\\{1,||u_{0}||_{L^{\infty}(\Omega)}\\}$. $(ii)$ follows from
multiplying (1.1) by $u$ and $tu_{t}$ respectively and performing the
integration by parts in both $x$ and $t$ over $\Omega\times(0,T)$. By applying
the same analysis to the Galerkin approximation of $u$, this formal argument
can be rendered rigorous. Hence, we omit the details. ∎

###  2.2 Weak solutions of effective problems

We define weak solutions of (1.4) together with the boundary conditions in
Table 1.

######  Definition 2.2.

Let the test function $\xi\in C^{\infty}(\overline{Q^{1}_{T}})$ satisfy
$\xi=0$ at $t=T$.

$(1)$ A function $v$ is said to be a weak solution of (1.4) with the Dirichlet
boundary condition $v=0$ if $v\in V^{1,0}_{2,0}(Q_{T}^{1})$, and for any test
function $\xi$, $v$ satisfies

$\begin{split}\mathcal{L}[v,\xi]:=&-\int_{\Omega_{1}}u_{0}(x)\xi(x,0)dx+\int_{0}^{T}\int_{\Omega_{1}}\left(k\nabla
v\cdot\nabla\xi-v\xi_{t}-f(v)\xi\right)dxdt=0.\end{split}$ (2.3)

$(2)$ A function $v$ is said to be a weak solution of (1.4) with the boundary
conditions $\nabla_{\Gamma}v=0$ and $\int_{\Gamma}(k\frac{\partial
v}{\partial\textbf{n}}+\alpha v)=0$ for $\alpha\in[0,\infty)$ if for almost
everywhere fixed $t\in(0,T)$, the trace of $v$ on $\Gamma$ is a constant, and
if $\nabla_{\Gamma}\xi=0$ on $\Gamma$, it holds that $v\in
V^{1,0}_{2}(Q_{T}^{1})$ and $v$ satisfies

$\mathcal{L}[v,\xi]=-\int_{0}^{T}\int_{\Gamma}\alpha v\xi dsdt.$

$(3)$ A function $v$ is said to be a weak solution of (1.4) with the boundary
condition $k\frac{\partial v}{\partial\textbf{n}}=\mathcal{B}[v]$, where
$\mathcal{B}[v]=-\alpha v$, or $\gamma\mathcal{J}^{h}[v]$ for
$h\in(0,\infty]$, if $v\in V^{1,0}_{2}(Q_{T}^{1})$ and if for any test
function $\xi$, $v$ satisfies

$\mathcal{L}[v,\xi]=\int_{0}^{T}\int_{\Gamma}v\mathcal{B}[\xi]dsdt.$

A weak solution of (1.4) satisfies the initial value in the sense that
$v(\cdot,t)\to u_{0}(\cdot)$ in $L^{2}(\Omega_{1})$ as $t\to 0.$

######  Theorem 2.1.

Suppose that $\Gamma\in C^{1}$ and $0\leq u_{0}\in L^{\infty}(\Omega_{1})$.
Then, (1.4) with any boundary condition in Tables 1 has one and only one weak
solution as defined in Definition 2.2.

###### Proof.

The theorem can be proved via the lower and upper solutions method by choosing
the lower solution $0$ and the upper solution
$v=\max\\{1,||u_{0}||_{L^{\infty}(\Omega_{1})}\\}$. ∎

We end this section by introducing the following parametrization of the
coating $\Omega_{2}$. Define $X$ by

$\Gamma\times(0,\delta)\longmapsto
x=X(p,r)=p+r\textbf{n}(p)\in\mathbb{R}^{3},$

where $p$ is the projection of $x$ on $\Gamma$; $\textbf{n}(p)$ is the unit
normal vector of $\Gamma$ pointing out of $\Omega_{1}$ at $p$; $r$ is the
distance from $x$ to $\Gamma$. As is well known ([7, Lemma 14.16]), for a
small $\delta>0$, $X$ is a $C^{1}$ smooth diffeomorphism from
$\Gamma\times(0,\delta)$ to $\Omega_{2}$; $r=r(x)$ is a $C^{2}$ smooth
function of $x$ and is seen as the inverse of the mapping $x=X(p,r)$.

By using local coordinates $s=(s_{1},s_{2})$ in a typical chart on $\Gamma$,
we then have

$x=X(p(s),r)=X(s,r),\quad dx=(1+2Hr+\kappa r^{2})dsdr\quad\text{
in}\quad\overline{\Omega}_{2},$ (2.4)

where $ds$ represents the surface element; $H(s)$and $\kappa(s)$ are the mean
curvature and Gaussian curvature at $p$ on $\Gamma$, respectively. Thus the
Riemannian metric tensor at $x\in\overline{\Omega}_{2}$ induced from
$\mathbb{R}^{3}$ is defined as $G(s,r)$ with elements

$g_{ij}(s,r)=g_{ji}(s,r)=<X_{i},X_{j}>_{\mathbb{R}^{3}},\quad i,j=1,2,3,$

where $X_{i}=X_{s_{i}}$ for $i=1,2$ and $X_{3}=X_{r}$. Let $|G|:=detG$ and
$g^{ij}(s,r)$ be the element of the inverse matrix of $G$, denoted by
$G^{-1}$.

In the curvilinear coordinates $(s,r)$, the derivatives of $u$ are given as

$\nabla u=u_{r}\textbf{n}+\nabla_{s}u,$ (2.5)

where

$\begin{split}\nabla_{s}u=\sum_{i,j=1,2}g^{ij}(s,r)u_{s_{j}}&X_{s_{i}}(s,r).\end{split}$
(2.6)

In particular, if $r=0$, then we have
$\nabla_{\Gamma}u=\sum_{i,j=1,2}g^{ij}(s,0)u_{s_{j}}p_{s_{i}}(s).$ Moreover,
if $A(x)$ satisfies (1.3), we have

$\nabla\cdot\left(A(x)\nabla
u\right)=\frac{\sigma}{\sqrt{|G|}}\left(\sqrt{|G|}u_{r}\right)_{r}+\mu\Delta_{s}u,$
(2.7)

where

$\begin{split}\Delta_{s}u=\nabla_{s}\cdot\nabla_{s}u&=\frac{1}{\sqrt{|G|}}\sum_{ij=1,2}\left(\sqrt{|G|}g^{ij}(s,r)u_{s_{i}}\right)_{s_{j}}.\end{split}$
(2.8)

Consequently, in $\overline{\Omega}_{2}$, we can rewrite $A(x)$ as

$A(x)=\sigma\textbf{n}(p)\otimes\textbf{n}(p)+\mu\sum_{ij}g^{ij}(s,r)X_{s_{i}}(s,r)\otimes
X_{s_{j}}(s,r).$ (2.9)

##  3 EBCs for the Fisher-KPP equation

The aim of this section is to prove Theorem 1.1 in which EBCs are obtained on
$\Gamma\times(0,T)$.

###  3.1 Auxiliary functions

Enlightened by [4], an auxiliary function $\psi$ is constructed in
$\Omega_{2}$. In the curvilinear coordinates $(s,r),$ for any $t\in[0,T]$ and
smooth $g(s)$ defined on $\Gamma$, let $\psi(s,r,t)$ be the bounded solution
of

$\left\\{\begin{array}[]{ll}\sigma\psi_{rr}+\mu\Delta_{\Gamma}\psi=0,&\Gamma\times(0,\delta),\\\
\psi(s,0,t)=g(s),&\psi(s,\delta,t)=0.\end{array}\right.$ (3.1)

It follows from the maximum principle that $\psi$ is unique. Multiplying (3.1)
by $\psi$ and integrating by parts over $\Gamma\times(0,\delta)$, we arrive at

$\int_{0}^{\delta}\int_{\Gamma}\left(\sigma\psi_{r}^{2}+\mu|\nabla_{\Gamma}\psi|^{2}\right)dsdr=-\int_{\Gamma}\sigma\psi_{r}(s,0,t)g(s)ds.$
(3.2)

Multiplying (3.1) by $u$ and performing the integration by parts again, we get

$\begin{split}\int_{0}^{\delta}\int_{\Gamma}\left(\sigma\psi_{r}u_{r}+\mu\nabla_{\Gamma}\psi\cdot\nabla_{\Gamma}u\right)dsdr&=-\int_{\Gamma}\sigma\psi_{r}(s,0,t)u(s,0,t)ds.\end{split}$
(3.3)

We next assume $r=R\sqrt{\sigma/\mu}$ and
$\Psi(s,R)=\psi(s,R\sqrt{\sigma/\mu},t).$ Suppressing $t$ and plugging $r$
into (3.1) lead to

$\left\\{\begin{array}[]{ll}\Psi_{RR}+\Delta_{\Gamma}\Psi=0,&\Gamma\times(0,h),\\\
\Psi(s,0)=g(s),&\Psi(s,h)=0,\end{array}\right.$ (3.4)

where
$h=\delta\sqrt{\frac{\mu}{\sigma}}=\frac{\mu\delta}{\sqrt{\sigma\mu}}=\frac{\sqrt{\sigma\mu}}{\sigma/\delta}.$
Then the Dirichlet-to-Neumann operator is defined by

$\mathcal{J}^{h}[g](s):=\Psi_{R}(s,0),$ (3.5)

from which

$\begin{split}\sigma\psi_{r}(s,0,t)&=\sqrt{\sigma\mu}\Psi_{R}(s,0)=\sqrt{\sigma\mu}\mathcal{J}^{h}[g](s).\end{split}$
(3.6)

By using separation of variables, $\mathcal{J}^{h}[g]$ can be given in
eigenvalues and eigenfunctions of the operator $-\Delta_{\Gamma}$. In
particular, we have

$\begin{split}\Psi(s,R)&=\sum_{n=1}^{\infty}\frac{-g_{n}e^{-\sqrt{\lambda_{n}}h}}{2sinh(\sqrt{\lambda_{n}}h)}\left(e^{\sqrt{\lambda_{n}}R}-e^{\sqrt{\lambda_{n}}(2h-R)}\right)e_{n}(s),\end{split}$
(3.7)

where $g_{n}=<e_{n},g>=\int_{\Gamma}e_{n}gds$; $\lambda_{n}$ and $e_{n}(s)$
are the eigenvalues and the corresponding eigenfunctions of $-\Delta_{\Gamma}$
defined on $\Gamma$. Thus, we obtain

$\begin{split}\mathcal{J}^{h}[g](s)=&-\sum_{n=1}^{\infty}\frac{\sqrt{\lambda_{n}}e_{n}(s)g_{n}}{tanh(\sqrt{\lambda_{n}}h)},\end{split}$
(3.8)

which results from (3.5) and (3.7). Furthermore, if $h\to H\in(0,\infty]$, we
then have

$\begin{split}\mathcal{J}^{h}[g](s)-\mathcal{J}^{H}[g](s)&=\sum_{n=1}^{\infty}\sqrt{\lambda_{n}}e_{n}(s)g_{n}\left(\frac{1}{tanh(\sqrt{\lambda_{n}}H)}-\frac{1}{tanh(\sqrt{\lambda_{n}}h)}\right)\\\
&=|H-h|\sum_{n=1}^{\infty}\lambda_{n}e_{n}(s)g_{n}\frac{-4}{(e^{\sqrt{\lambda_{n}}h^{\prime}}-e^{-\sqrt{\lambda_{n}}h^{\prime}})^{2}}\\\
&=O(|H-h|),\end{split}$

for some $h^{\prime}$ between $h$ and $H$, which means that
$\mathcal{J}^{h}[g]\to\mathcal{J}^{H}[g]$ uniformly as $h\to H$, and
$\mathcal{J}^{\infty}[g]=-(-\Delta_{\Gamma})^{1/2}g.$

In what follows, we will estimate the size of the term $\Psi_{R}(s,0)$ for a
sufficiently small $\delta$. If $h\to 0$ as $\delta\to 0$, then it holds from
(3.8) that

$\begin{split}\left|\Psi_{R}(s,0)+\frac{g(s)}{h}\right|&\leq
h\|g\|_{C^{2}(\Gamma)}.\end{split}$

Combining this with (3.6), we are led to

$\begin{split}\sqrt{\sigma\mu}\Psi_{R}(s,0)&=\frac{\sigma}{\delta}\left(-g(s)+O(h^{2})\right).\end{split}$
(3.9)

On the other hand, if $h\to H\in(0,\infty]$ as $\delta\to 0$, then from the
Taylor expansion for $\Psi(s,R)$, we obtain

$\Psi_{R}(s,0)=\frac{\Psi(s,R)-\Psi(s,0)}{R}-\frac{R}{2}\Psi_{RR}(s,\overline{R}),$

for some $\overline{R}\in[0,R]$. Then, it holds from the maximum principle
that

$\begin{split}\|\Psi_{R}(s,0)\|_{L^{\infty}(\Gamma))}&\leq\frac{2}{R}\|\Psi\|_{L^{\infty}(\Omega_{2})}+R\|\Psi_{RR}\|_{L^{\infty}(\Omega_{2})}\leq\frac{3\|g\|_{C^{2}(\Gamma)}}{R},\end{split}$

from which it turns out that

$\sqrt{\sigma\mu}\|\Psi_{R}\|_{L^{\infty}(\Gamma)}=\frac{O(1)\sqrt{\sigma\mu}}{R}.$
(3.10)

######  Remark 3.1.

For $h\in(0,\infty)$, $\mathcal{J}^{h}[g]$ is defined for smooth $g$. However,
it is worth noting that it is also well-defined for $g\in
H^{\frac{1}{2}}(\Gamma)$ where $H^{\frac{1}{2}}(\Gamma)$ is defined by the
completion of smooth functions under the $H^{\frac{1}{2}}(\Gamma)$ norm.
Moreover, $\mathcal{J}^{H}:H^{\frac{1}{2}}(\Gamma)\to
H^{-\frac{1}{2}}(\Gamma)$ is linear and symmetric with
$H^{-\frac{1}{2}}(\Gamma)$ being the dual space of $H^{\frac{1}{2}}(\Gamma)$.

###  3.2 Main proof of Theorem 1.1

The main proof consists of two steps. Firstly, we establish the weak
convergence of $u\to$ some $v$ in
$W^{1,0}_{2}\left(\Omega_{1}\times(0,T)\right)$, and strong convergence of
that in $C([0,T];H^{1}(\Omega_{1}))$ after passing to a subsequence of
$\delta>0$. Subsequently, we use the auxiliary function $\psi$ to show that
$v$ is a weak solution of (1.4) subject to the effective boundary conditions
in Table 1.

###### Proof of Theorem 1.1.

Step 1. We begin with the proof by considering the compactness of weak
solutions $\\{u\\}_{\delta>0}$ of (1.1).

By Lemma 2.1, $\\{u\\}_{\delta>0}$ is bounded in
$W^{1,0}_{2}\left(\Omega_{1}\times(0,T)\right)$ and
$C([t_{0},T];H^{1}(\Omega_{1}))$ for any small $t_{0}\in(0,T]$. By Banach-
Eberlein theorem, $u\to$ some $v$ weakly in
$W^{1,0}_{2}\left(\Omega_{1}\times(0,T)\right)$ and
$C([t_{0},T];H^{1}_{0}(\Omega_{1}))$ after passing to a subsequence of
$\delta\to 0.$ Since the embedding $H^{1}(\Omega_{1})\hookrightarrow
L^{2}(\Omega_{1})$ is compact, for any fixed $t_{0}$, $\\{u\\}_{\delta>0}$ is
pre-compact in $L^{2}(\Omega_{1}).$ Moreover, the functions
$\\{u\\}_{\delta>0}$ : $t\in[t_{0},T]\mapsto u(\cdot,t)\in L^{2}(\Omega_{1})$
are equicontinuous because it holds from Lemma 2.1 that the term
$\int_{Q_{T}}tu_{t}^{2}dxdt$ is bounded. Consequently, by the generalized
Arzela-Ascoli theorem, after passing to a further subsequence of $\delta\to
0,$ $u\to v$ strongly in $C\left([t_{0},T];L^{2}(\Omega_{1})\right)$.

In the following, it suffices to prove that the strong convergence is in
$C\left([0,T];L^{2}(\Omega_{1})\right)$. We then take a sequence $u_{0}^{n}\in
C_{0}^{\infty}(\Omega_{1})$ that can be constructed by multiplying $u_{0}$ by
cut-off functions in the outer normal direction of $\Gamma$ such that
$u_{0}^{n}\equiv 0$ in $\Omega_{2}$,
$\|u_{0}-u_{0}^{n}\|_{L^{2}(\Omega)}\leq\frac{1}{n}+\|u_{0}\|_{L^{2}(\Omega_{2})}$,
and $\|\nabla u_{0}^{n}\|_{L^{2}(\Omega)}\leq C(n)$. Then, we decompose
$u=u_{1}+u_{2}$, where $u_{1}$ is the weak solution of (1.1) with $f\equiv 0$
and the initial value replaced by $u_{0}-u_{0}^{n}$, and $u_{2}$ is the weak
solution of (1.1) with the initial value replaced by $u_{0}^{n}$.

Since the equation concerning $u_{1}$ is homogeneous, we have

$\|u_{1}(\cdot,t)\|_{L^{2}(\Omega)}\leq\|u_{0}-u_{0}^{n}\|_{L^{2}(\Omega)}\leq\frac{1}{n}+\|u_{0}\|_{L^{2}(\Omega_{2})}.$
(3.11)

Multiplying the equation concerning $u_{2}$ by $(u_{2})_{t}$ and performing
integration by parts, we get

$\begin{split}\int_{0}^{t}\int_{\Omega}(u_{2})^{2}_{t}dxd\tau+\int_{\Omega}\nabla
u_{2}(x,t)\cdot A(x)\nabla
u_{2}(x,t)dx\leq&\int_{0}^{t}\int_{\Omega}f^{2}(u)dxd\tau+\int_{\Omega}\nabla
u_{0}^{n}\cdot A(x)\nabla u_{0}^{n}dx\\\
\leq&\int_{0}^{t}\int_{\Omega}f^{2}(u)dxd\tau+k\int_{\Omega_{1}}|\nabla
u_{0}^{n}|^{2}dx\\\ =:&F(f(u),n).\end{split}$

For any $t\in[0,t_{0}]$, it follows from this and (3.11) that

$\begin{split}\|u_{2}(\cdot,t)-u_{0}^{n}(\cdot)\|^{2}_{L^{2}(\Omega)}&=2\int_{0}^{t}\int_{\Omega}\left(u_{2}(x,t)-u_{0}^{n}(x)\right)(u_{2})_{t}dxd\tau\\\
&\leq
2\left(\int_{0}^{t}\int_{\Omega}\left(u_{2}(x,t)-u_{0}^{n}(x)\right)^{2}\right)^{\frac{1}{2}}\left(\int_{0}^{t}\int_{\Omega}(u_{2})_{t}^{2}\right)^{\frac{1}{2}}\\\
&\leq
2\sqrt{t_{0}}\max_{t\in[0,t_{0}]}\|u_{2}(\cdot,t)-u_{0}^{n}(\cdot)\|_{L^{2}(\Omega)}\left(F(f(u),n)\right)^{\frac{1}{2}},\end{split}$

resulting in

$\max_{t\in[0,t_{0}]}\|u_{2}(\cdot,t)-u_{0}^{n}(\cdot)\|\leq
2\sqrt{t_{0}}\left(F(f(u),n)\right)^{\frac{1}{2}}.$ (3.12)

Finally, it holds from (3.11) and (3.12) that

$\begin{split}\|u(\cdot,t)-u_{0}(\cdot)\|_{L^{2}(\Omega_{1})}&\leq\|u_{1}(\cdot,t)\|_{L^{2}(\Omega_{1})}+\|u_{2}(\cdot,t)-u_{0}(\cdot)\|_{L^{2}(\Omega_{1})}+\|u_{0}-u_{0}(\cdot)\|_{L^{2}(\Omega_{1})}\\\
\
&\leq\frac{2}{n}+2\|u_{0}\|_{L^{2}(\Omega_{2})}+2\sqrt{t_{0}}\left(F(f(u),n)\right)^{\frac{1}{2}},\end{split}$

for any $n$ and sufficiently small $t_{0}$ and $\delta$. Thus,
$\|u(\cdot,t)-u_{0}(\cdot)\|_{L^{2}(\Omega_{1})}$ can be arbitrary small for
$t\in[0,t_{0}]$. Combining the fact that $u\to v$ strongly in
$C\left([t_{0},T];L^{2}(\Omega_{1})\right)$, we conclude that $u\to v$
strongly in $C\left([0,T];L^{2}(\Omega_{1})\right)$ if we define
$v(\cdot,0)=u_{0}.$

Step 2. In what follows, we aim to find all effective boundary conditions in
Table 1.

By the compactness argument in Step 1, given any subsequence of $\delta$, we
can ensure that $u\to$ some $v$ in all above spaces after passing to a further
subsequence. In this step, we will show that $v$ is a weak solution of (1.4)
with effective boundary conditions listed in Table 1. Since $v$ is unique by
Theorem 2.1, $u\to v$ without passing to any subsequence of $\delta>0$.

Let the test function $\xi\in C^{\infty}(\overline{\Omega}_{1}\times[0,T])$
satisfy $\xi=0$ at $t=T$, and extend $\xi$ to $\overline{\Omega}\times[0,T]$
by defining

$\overline{\xi}(x,t)=\left\\{\begin{aligned}
&\xi(x,t),&x\in\overline{\Omega}_{1},\\\
&\psi(p(x),r(x),t),&x\in\Omega_{2},\end{aligned}\right.$

where $\psi$ is introduced in (3.1) with $g(s)=\xi(s,0,t)$. It is easy to see
$\overline{\xi}\in W^{1,1}_{2,0}(Q_{T})$, and $\overline{\xi}$ is called the
harmonic extension of $\xi$. According to Definition 2.1, it holds that

$\begin{split}\mathcal{A}[u,\overline{\xi}]&=-\int_{\Omega}u_{0}(x)\overline{\xi}(x,0)dx+\int_{0}^{T}\int_{\Omega}\left(\nabla\overline{\xi}\cdot
A\nabla u-u\overline{\xi}_{t}-f(u)\overline{\xi}\right)dxdt=0,\end{split}$

which is equivalent to

$\begin{split}&\int_{0}^{T}\int_{\Omega_{1}}k\nabla\xi\cdot\nabla
udxdt-\int_{\Omega}u_{0}(x)\overline{\xi}(x,0)dx-\int_{0}^{T}\int_{\Omega}(u\overline{\xi_{t}}+f(u)\overline{\xi})dxdt\\\
=&-\int_{0}^{T}\int_{\Omega_{2}}\nabla\psi\cdot A(x)\nabla udxdt.\end{split}$
(3.13)

By what we have proven in Step 1, $u\to v$ weakly in
$W^{1,0}_{2}\left(\Omega_{1}\times(0,T)\right)$, and strongly in
$C\left([0,T];L^{2}(\Omega_{1})\right)$ as $\delta\to 0$. Thus, it gives rise
to

$\left\\{\begin{array}[]{ll}&\int_{Q_{T}}u\xi_{t}dxdt\rightarrow\int_{Q^{1}_{T}}v\xi_{t}dxdt,\\\
&\int_{Q^{1}_{T}}\nabla u\cdot\nabla\xi dxdt\to\int_{Q^{1}_{T}}\nabla
v\cdot\nabla\xi dxdt,\\\
&\int_{Q_{T}}f(u)\overline{\xi}dxdt\to\int_{Q^{1}_{T}}f(v)\xi
dxdt.\end{array}\right.$

Then the left-hand side of (3.13) is rewritten as

$\begin{split}\mathcal{L}[v,\xi]:=&\int_{0}^{T}\int_{\Omega_{1}}k\nabla\xi\cdot\nabla
vdxdt-\int_{\Omega_{1}}u_{0}(x)\xi(x,0)dx-\int_{0}^{T}\int_{\Omega_{1}}(v\xi_{t}+f(v)\xi)dxdt.\end{split}$
(3.14)

The remainder of the following focuses on the right-hand side of (3.13). In
the curvilinear coordinates $(s,r)$, due to (2.4), (2.5) and (2.9), we have

$\begin{split}RHS:=&-\int_{0}^{T}\int_{\Omega_{2}}\nabla\psi\cdot A(x)\nabla
udxdt\\\
=&-\int_{0}^{T}\int_{\Gamma}\int_{0}^{\delta}\left(\sigma\psi_{r}u_{r}+\mu\nabla_{s}\psi\nabla_{s}u\right)\left(1+2H(s)r+\kappa(s)r^{2}\right)drdsdt\\\
=&-\int_{0}^{T}\int_{\Gamma}\int_{0}^{\delta}\left(\sigma\psi_{r}u_{r}+\mu\nabla_{\Gamma}\psi\nabla_{\Gamma}u\right)-\int_{0}^{T}\int_{\Gamma}\int_{0}^{\delta}(\sigma\psi_{r}u_{r}+\mu\nabla_{\Gamma}\psi\nabla_{\Gamma}u)(2Hr+\kappa
r^{2})\\\
&-\int_{0}^{T}\int_{\Gamma}\int_{0}^{\delta}\mu(\nabla_{s}\psi\nabla_{s}u-\nabla_{\Gamma}\psi\nabla_{\Gamma}u)(1+2Hr+\kappa
r^{2})\\\ =&I+II+III.\end{split}$ (3.15)

Firstly, it follows from (3.3) and (3.6) that

$I=\sqrt{\sigma\mu}\int_{0}^{T}\int_{\Gamma}u(s,0,t)\Psi_{R}(s,0)dsdt.$ (3.16)

Secondly, in view of (3.2), (3.6) and Lemma 2.1, we get

$\begin{split}|II|=&\left|\int_{0}^{T}\int_{\Gamma}\int_{0}^{\delta}(\sigma\psi_{r}u_{r}+\mu\nabla_{\Gamma}\psi\nabla_{\Gamma}u)(2Hr+\kappa
r^{2})drdsdt\right|\\\
=&O(\delta)\int_{0}^{T}\left(\int_{\Gamma}\int_{0}^{\delta}\sigma\psi_{r}^{2}+\mu|\nabla_{\Gamma}\psi|^{2}\right)^{1/2}\left(\int_{\Omega}\sigma
u_{r}^{2}+\mu|\nabla_{\Gamma}u|^{2}\right)^{1/2}dt\\\
=&O(\delta)\int_{0}^{T}\frac{1}{\sqrt{t}}\left(\int_{\Gamma}\sigma|\psi_{r}(s,0,t)|\right)^{1/2}dt\\\
=&O(\delta)\sqrt{T}(\sigma\mu)^{1/4}||\Psi_{R}(s,0)||^{1/2}_{L^{\infty}(\Gamma)},\end{split}$
(3.17)

where we have used Höder inequality. Thirdly, according to (2.6), (3.2), (3.6)
and Lemma 2.1, it turns out that

$\begin{split}|III|=&\left|\int_{0}^{T}\int_{\Gamma}\int_{0}^{\delta}\mu(\nabla_{s}\psi\nabla_{s}u-\nabla_{\Gamma}\psi\nabla_{\Gamma}u)(1+2Hr+\kappa
r^{2})drdsdt\right|\\\
=&O(\delta)\left|\int_{0}^{T}\int_{\Gamma}\int_{0}^{\delta}\mu\sum_{ij}\psi_{s_{i}}u_{s_{j}}drdsdt\right|\\\
=&O(\delta)\int_{0}^{T}\left(\int_{\Gamma}\int_{0}^{\delta}\sigma\psi_{r}^{2}+\mu|\nabla_{\Gamma}\psi|^{2}\right)^{1/2}\left(\int_{\Omega}\sigma
u_{r}^{2}+\mu|\nabla_{\Gamma}u|^{2}\right)^{1/2}dt\\\
=&O(\delta)\sqrt{T}(\sigma\mu)^{1/4}||\Psi_{R}(s,0)||^{1/2}_{L^{\infty}(\Gamma)},\end{split}$
(3.18)

where the Taylor expansion for $g(s,r)$ and Höder inequality were used.

We next investigate the asymptotic behavior of RHS as $\delta\to 0$. To this
end, we consider the following cases

$(1)\frac{\sigma}{\delta}\to
0,\quad(2)\frac{\sigma}{\delta}\to\alpha\in(0,\infty),\quad(3)\frac{\sigma}{\delta}\to\infty.$

Case $1$. $\frac{\sigma}{\delta}\to 0$ as $\delta\to 0$.

Subcase $(1i)$. $\sigma\mu\to 0$ as $\delta\to 0$. By the boundedness of $u$,
(3.9) and (3.10), (3.16) can be rewritten as

$|I|\leq C(T)\max\\{\frac{\sigma}{\delta},(\sigma\mu)^{1/2}\\}.$

Similarly, it holds from (3.9) and (3.10) that

$|II|+|III|=O(\delta)\sqrt{T}\max\\{\sqrt{\frac{\sigma}{\delta}},(\sigma\mu)^{1/4}\\}.$

Combining these with (3.15), we get

$\begin{split}|RHS|\leq&|I|+|II|+|III|\to 0\;\text{ as }\;\delta\to
0,\end{split}$

from which we have $\mathcal{L}[v,\xi]=0$, indicating that $v$ satisfies
$\frac{\partial v}{\partial\textbf{n}}=0$ on $\Gamma\times(0,T)$.

Subcase $(1ii)$. $\sqrt{\sigma\mu}\to\gamma\in(0,\infty)$ as $\delta\to 0$. In
this case, $h\to\infty$. By the weak convergence of $u$, it follows from (3.8)
and (3.16) that

$\begin{split}I&\longrightarrow\gamma\int_{0}^{T}\int_{\Gamma}v\mathcal{J}^{\infty}[g]dsdt\;\text{
as }\;\delta\to 0.\end{split}$

Due to (3.10), (3.17) and (3.18), it gives

$|II|+|III|\to 0\;\text{ as }\;\delta\to 0,$

from which we are led to

$\mathcal{L}[v,\xi]=\gamma\int_{0}^{T}\int_{\Gamma}v\mathcal{J}^{\infty}[\xi].$

Thus, $v$ satisfies the effective boundary condition $k\frac{\partial
v}{\partial\textbf{n}}=\gamma\mathcal{J}^{\infty}[v]$ on $\Gamma\times(0,T)$.

Subcase $(1iii)$. $\sigma\mu\to\infty$ as $\delta\to 0$. In this case,
$h\to\infty$. Divided both sides of (3.13) by $\sqrt{\sigma\mu}$ and sending
$\delta\to 0$, we get

$\int_{0}^{T}\int_{\Gamma}v\mathcal{J}^{\infty}[g]=0.$

Since the range of $\mathcal{J}^{\infty}[\cdot]$ contains
$\\{e_{n}\\}_{n=1}^{\infty}$, it implies that $\nabla_{\Gamma}v=0$ on $\Gamma$
for almost everywhere $t\in(0,T)$.

Now, we further take the test function $\xi$ satisfying $\xi(s,0,t)=m(t)$ for
some smooth function $m(t)$, and then define a new auxiliary function $\psi$
as

$\psi(s,r,t)=(1-\frac{r}{\delta})m(t).$

Consequently, direct computation leads to

$\begin{split}RHS=&-\int_{0}^{T}\int_{\Omega_{2}}\nabla\psi\cdot A\nabla
udxdt\\\ =&\int_{0}^{T}\frac{\sigma
m(t)}{\delta}\left(\int_{0}^{\delta}\int_{\Gamma}u_{r}(1+2Hr+\kappa
r^{2})dsdr\right)dt\\\ =&\int_{0}^{T}\frac{\sigma
m(t)}{\delta}\left(\int_{\Gamma}u(s,0,t)ds\right)dt+\int_{0}^{T}\frac{\sigma
m(t)}{\delta}\int_{0}^{\delta}\int_{\Gamma}u(2H+2\kappa r)dsdr\\\
=&\frac{\sigma}{\delta}\int_{0}^{T}m(t)dt\left(O(1)+O(\delta)\|u\|_{L^{\infty}(Q_{T})}\right)\to
0\;\text{ as }\;\delta\to 0,\end{split}$ (3.19)

from which we have $\mathcal{L}[v,\xi]=0$ as $\delta\to 0$. Thus, $v$
satisfies $\int_{\Gamma}\frac{\partial v}{\partial\textbf{n}}ds=0$ on
$\Gamma\times(0,T)$ due to the arbitrariness of the test function $\xi$.

Case $2$. $\frac{\sigma}{\delta}\to\alpha\in(0,\infty)$ as $\delta\to 0$.

Subcase $(2i)$. $\sigma\mu\to 0$ as $\delta\to 0$. In this case, $h\to 0$. It
follows from (3.9) and (3.16) that

$I\to-\alpha\int_{0}^{T}\int_{\Gamma}v\xi\;\text{ as }\delta\to 0.$

Subsequently, in view of (3.9), (3.17) and (3.18), we have

$|II+III|\to 0\;\text{ as }\;\delta\to 0,$

from which we are led to
$\mathcal{L}[v,\xi]=-\alpha\int_{0}^{T}\int_{\Gamma}v\xi,$ meaning that $v$
satisfies $k\frac{\partial v}{\partial\textbf{n}}=-\alpha v$ on
$\Gamma\times(0,T)$.

Subcase $(2ii)$. $\sqrt{\sigma\mu}\to\gamma\in(0,\infty)$ as $\delta\to 0$. In
this case, $h\to H=\gamma/\alpha\in(0,\infty)$. By virtue of (3.9) and (3.15)-
(3.18), it holds that

$I\longrightarrow\gamma\int_{0}^{T}\int_{\Gamma}v\mathcal{J}^{\gamma/\alpha}[\xi]\quad\text{
and }\quad|II+III|\longrightarrow 0\text{ as }\delta\to 0,$

from which we get
$\mathcal{L}[v,\xi]=\gamma\int_{0}^{T}\int_{\Gamma}v\mathcal{J}_{D}^{\gamma/\alpha}[\xi],$
implying that $v$ satisfies $k\frac{\partial
v}{\partial\textbf{n}}=\gamma\mathcal{J}^{\gamma/\alpha}[v]$ on
$\Gamma\times(0,T)$.

Subcase $(2iii)$. $\sigma\mu\to\infty$ as $\delta\to 0$. In this case,
$h\to\infty$. Divided both sides of (3.13) by $\sqrt{\sigma\mu}$ and sending
$\delta\to 0$, we obtain

$\int_{0}^{T}\int_{\Gamma}v\mathcal{J}^{\infty}[\xi]=0,$

from which we get $\nabla_{\Gamma}v=0$ on $\Gamma$ by using the same argument
as in Subcase $(1iii)$. Then take the same test function $\xi$ and the same
auxiliary function $\psi$ as in Subcase $(1iii)$, and thus it follows from
(3.19) that

$\mathcal{L}[v,\xi]=-\alpha\int_{0}^{T}\int_{\Gamma}v\xi,$

which means $v$ satisfies $\int_{\Gamma}\left(k\frac{\partial
v}{\partial\textbf{n}}+\alpha v\right)=0$ on $\Gamma\times(0,T)$.

Case $3$. $\frac{\sigma}{\delta}\to\infty$ as $\delta\to 0$.

Subcase $(3i)$. $\sqrt{\sigma\mu}\to\gamma\in[0,\infty)$ as $\delta\to 0$. In
this case, $h\to 0$. Divided both sides of (3.13) by $\sigma/\delta$ and
sending $\delta\to 0$, in view of (3.2) and (3.15)- (3.18), we are led to

$\frac{\delta}{\sigma}I\longrightarrow-\int_{0}^{T}\int_{\Gamma}v\xi=0,$

from which $v=0$ on $\Gamma\times(0,T)$ since $\xi$ is arbitrary.

Subcase $(3ii)$. $\sigma\mu\to\infty$ as $\delta\to 0$. In this case, after
passing to a subsequence, we have $h\to H\in[0,\infty]$. If $H=0$, then
divided both sides of (3.13) by $\sigma/\delta$ and sending $\delta\to 0$, we
have $v=0$ on $\Gamma\times(0,T)$.

On the other hand, if $H\in(0,\infty]$, then divided both sides of (3.13) by
$\sqrt{\sigma\mu}$ and sending $\delta\to 0$, we have

$\frac{I}{\sqrt{\sigma\mu}}\longrightarrow\int_{0}^{T}\int_{\Gamma}v\mathcal{J}^{H}[\xi]=0.$

By using the same argument as in Subcase $(1iii)$, for almost everywhere
$t\in(0,T)$, we have $\nabla_{\Gamma}v=0$, and

$\int_{0}^{T}\int_{\Gamma}vm(t)=0,$

from which we have $v=0$ on $\Gamma\times(0,T)$ for $m(t)$ is an arbitrary
function in $t$.

Therefore, we complete the whole proof of Theorem 1.1. ∎

##  4 Lifespan of EBCs

In this section, we delve into the lifespan of each EBC in Table 1 to
determinate how long these EBCs remain effective. Moreover, we prove that the
lifespan of each EBC is infinite, which means that the convergence of
$u(\cdot,t)\to v(\cdot,t)$ in $L^{2}(\Omega_{1})$ is global in $t$.

###  4.1 Eigenvalue problems

Consider the eigenvalue problem of

$\left\\{\begin{array}[]{llr}-\nabla\cdot(A(x)\nabla u)=\lambda
u,&\mbox{$x\in\Omega,$}\\\ u=0,&\mbox{$x\in\partial\Omega,$}&\\\
\end{array}\right.$ (4.1)

where $A(x)$ is given in (1.2) and (1.3). Let
$(\lambda_{1}^{\delta},e_{1}^{\delta})$ be an eigenpair of (4.1) such that
$e_{1}^{\delta}>0$ in $\Omega$ and $||e_{1}^{\delta}||_{L^{2}(\Omega)}=1$. It
is well-known that

$\lambda_{1}^{\delta}=\underset{0\not\equiv u\in
H^{1}_{0}(\Omega)}{\inf}\frac{\int_{\Omega}A\nabla u\cdot\nabla
udx}{\int_{\Omega}u^{2}dx}.$ (4.2)

######  Lemma 4.1.

Suppose that $\frac{\sigma}{\delta}\to 0$ as $\delta\to 0$. Then, for any
small $\delta>0$, it holds,

$C_{1}\frac{\sigma}{\delta}\leq\lambda_{1}^{\delta}\leq
C_{2}\frac{\sigma}{\delta},$

and

$\int_{\Omega_{1}}\left(e_{1}^{\delta}-\frac{1}{\sqrt{|\Omega_{1}|}}\right)^{2}dx\leq
C\left(\frac{\sigma}{\delta}+\delta^{2}\right),$

where $C_{1},C_{2}$ and $C$ are positive constants independent of $\delta$.

###### Proof.

The can be found in [9, Lemma 2.1], and hence we omit the details. ∎

######  Theorem 4.1.

Suppose $A(x)$ is given in (1.2) and (1.3). Assume that $\sigma$ and $\mu$
satisfy

$\lim_{\delta\to 0}\sigma\mu=\gamma\in[0,\infty],\quad\lim_{\delta\to
0}\frac{\sigma}{\delta}=\alpha\in[0,\infty].$

Then, as $\delta\to 0$,

$\lambda_{1}^{\delta}\to\lambda_{1}\quad\text{ and }\quad e_{1}^{\delta}\to
e_{1}\quad\text{ in }L^{2}(\Omega_{1}),$

where $e_{1}>0$ with $||e_{1}||_{L^{2}(\Omega_{1})}=1$, and
$(\lambda_{1},e_{1})$ is the principal eigenpair of

$-k\Delta v=\lambda v,\quad x\in\Omega_{1}$ (4.3)

subject to the effective boundary conditions listed in Table 1.

###### Proof of Theorem 4.1.

Denote by $(\lambda_{1}^{D},e_{1}^{D})$ the principal eigenpair of (4.3)
subject to the Dirichlet boundary condition. Taking a special test function
$\phi\in H^{1}_{0}(\Omega)$ in (4.2) as

$\phi=\left\\{\begin{array}[]{llr}e_{1}^{D},&\mbox{$x\in\Omega_{1},$}\\\
0,&\mbox{$x\in\Omega_{2},$}&\end{array}\right.$ (4.4)

we are led to

$\lambda_{1}^{\delta}<\frac{k\int_{\Omega_{1}}|\nabla
e_{1}^{D}|^{2}dx}{\int_{\Omega_{1}}(e_{1}^{D})^{2}dx}=\lambda_{1}^{D}.$

Moreover, it holds from (4.2) that $\lambda_{1}^{\delta}>0$. Then after
passing to a subsequence of $\delta>0$, we have

$\lambda_{1}^{\delta}\to\;\text{ some }\;\lambda_{1}^{0}.$

Subsequently, it is easy to see that

$\lambda_{1}^{\delta}=k\int_{\Omega_{1}}|\nabla
e_{1}^{\delta}|^{2}dx+\int_{\Omega_{2}}\sigma(e_{1}^{\delta})_{r}^{2}+\mu|\nabla_{s}e_{1}^{\delta}|^{2}dx\leq\lambda_{1}^{D},$
(4.5)

from which $\\{e_{1}^{\delta}\\}_{\delta>0}$ is bounded in
$H^{1}(\Omega_{1})$. Thus, after passing to a further subsequence of
$\delta\to 0$, we have

$e_{1}^{\delta}\to\;\text{some}\;e_{1}^{0}$

weakly in $H^{1}(\Omega_{1})$ and strongly in $L^{2}(\Omega_{1})$.

Consequently, similar to the parabolic case in the proof of Theorem 1.1, we
can show that $e_{1}^{0}$ is a weak solution of (4.3) subject to the effective
boundary conditions in Table 1 with the corresponding $\lambda_{1}^{0}$. We
next prove that $e_{1}^{0}\not\equiv 0$, after which it holds that $e^{0}_{1}$
is the eigenfunction of $\lambda_{1}^{0}$.

In view of (2.4), we get

$\begin{split}\int_{\Omega_{2}}(e_{1}^{\delta})^{2}dx=&\int_{\Gamma}\int_{0}^{\delta}(e_{1}^{\delta})^{2}(1+H(s)r+\kappa(s)r^{2})drds\\\
=&\int_{\Gamma}\int_{0}^{\delta}\left(\int_{r}^{\delta}(e_{1}^{\delta})_{\tau}d\tau\right)^{2}(1+H(s)r+\kappa(s)r^{2})drds\\\
\leq&C\delta\int_{\Gamma}\int_{0}^{\delta}\left(\int_{0}^{\delta}(e_{1}^{\delta})^{2}_{\tau}d\tau\right)(1+H(s)r+\kappa(s)r^{2})drds\\\
\leq&C\delta^{2}\int_{\Gamma}\int_{0}^{\delta}(e_{1}^{\delta})^{2}_{r}drds,\end{split}$
(4.6)

from which, together with (4.5), we have

$\int_{\Omega_{2}}(e_{1}^{\delta})^{2}dx\leq C\frac{\delta^{2}}{\sigma}\to
0\;\text{ as }\;\delta\to 0,$

in the case of $\frac{\sigma}{\delta}\to\alpha\in(0,\infty]$. Combining (4.5)
and (4.6) with Lemma 4.1, as $\delta\to 0$, we have

$\int_{\Omega_{2}}(e_{1}^{\delta})^{2}dx\leq C\delta\to 0$

in the case of $\frac{\sigma}{\delta}\to 0$. Thus, we obtain
$e_{1}^{0}\not\equiv 0$ and $||e_{1}^{0}||_{L^{2}(\Omega_{1})}=1$, resulting
from

$1=\int_{\Omega}(e_{1}^{\delta})^{2}dx=\int_{\Omega_{1}}(e_{1}^{\delta})^{2}dx+\int_{\Omega_{2}}(e_{1}^{\delta})^{2}dx.$

Consequently, by the uniqueness of the principal eigenpair, it follows that
$(\lambda_{1},e_{1})=(\lambda_{1}^{0},e_{1}^{0})$ and all above convergences
hold without passing to a subsequence. ∎

We next consider the eigenvalue problems of the operator $-k\Delta$ subject to
some nonlocal boundary conditions in Table 1.

######  Lemma 4.2.

(i) Consider the eigenvalue problem of

$\left\\{\begin{array}[]{llr}-k\Delta v=\mu v,&\mbox{$x\in\Omega_{1},$}\\\
k\frac{\partial
v}{\partial\textbf{n}}=\gamma\mathcal{J}^{\infty}[v],&\mbox{$x\in\partial\Omega_{1},$}&\\\
\end{array}\right.$ (4.7)

where $\alpha=0$ and $\gamma\in(0,\infty)$. Then, the principle eigenvalue of
(4.7) is $0$, and the second eigenvalue is positive.

(ii) Consider the eigenvalue problem of

$\left\\{\begin{array}[]{llr}-k\Delta v=\mu v,&\mbox{$x\in\Omega_{1},$}\\\
\nabla_{\Gamma}v=0,\;\int_{\partial\Omega_{1}}\frac{\partial
v}{\partial\textbf{n}}ds=0,&\mbox{$x\in\partial\Omega_{1},$}&\\\
\end{array}\right.$ (4.8)

where $\alpha=0$ and $\gamma=\infty$. Then, the principle eigenvalue of (4.8)
is $0$, and the second eigenvalue is positive.

(iii) Consider the eigenvalue problem of

$\left\\{\begin{array}[]{llr}-k\Delta v=\mu v,&\mbox{$x\in\Omega_{1},$}\\\
k\frac{\partial
v}{\partial\textbf{n}}=\gamma\mathcal{J}^{\gamma/\alpha}[v],&\mbox{$x\in\partial\Omega_{1},$}&\\\
\end{array}\right.$ (4.9)

where $\alpha\in(0,\infty)$ and $\gamma\in(0,\infty)$. Then, the principle
eigenvalue of (4.9) is positive.

(iv) Consider the eigenvalue problem of

$\left\\{\begin{array}[]{llr}-k\Delta v=\mu v,&\mbox{$x\in\Omega_{1},$}\\\
\nabla_{\Gamma}v=0,\;\int_{\partial\Omega_{1}}\left(k\frac{\partial
v}{\partial\textbf{n}}+\alpha
v\right)ds=0,&\mbox{$x\in\partial\Omega_{1},$}&\\\ \end{array}\right.$ (4.10)

where $\alpha\in(0,\infty)$ and $\gamma=\infty$ . Then, the principle
eigenvalue of (4.10) is positive.

###### Proof.

The proof of this Lemma can be found in [9, Lemma 2.2-2.3], and hence we omit
the details. ∎

We now turn to the positive steady state of (1.1) as

$\left\\{\begin{array}[]{llr}-\nabla\cdot(A(x)\nabla
U)=U(1-U),&\mbox{$x\in\Omega,$}\\\ U=0,&\mbox{$x\in\partial\Omega,$}&\\\
\end{array}\right.$ (4.11)

where $U$ is the unique positive solution. It is well-known that $(i)$ the
positive steady state $U$ exists if and only if the principal eigenvalue
$\lambda_{1}^{\delta}<1$; $(ii)$ if $U$ exists, then it is unique; $(iii)$ as
$t\to\infty$, $u\to U$ if $U$ exists ($u\to 0$ if $U$ does not exist).

######  Theorem 4.2.

Suppose that $A(x)$ is given in (1.2) and (1.3). Assume that $\sigma$ and
$\mu$ satisfy

$\lim_{\delta\to 0}\sigma\mu=\gamma\in[0,\infty],\quad\lim_{\delta\to
0}\frac{\sigma}{\delta}=\alpha\in[0,\infty].$

Then, as $\delta\to 0$, $U\to V$ in $L^{2}(\Omega_{1})$, where $V$ is the
unique positive solution of

$-k\Delta V=V(1-V),\;x\in\Omega_{1}$ (4.12)

subject to the effective boundary conditions listed in Table 1 (with $v$
replaced by $V$).

######  Remark 4.1.

As the results about $U$, it holds that the positive steady state $V$ exists
if and only if the principal eigenvalue $\lambda_{1}<1$; if $V$ exists, then
it is unique; as $t\to\infty$, $v\to V$ if $V$ exists ($v\to 0$ if $V$ does
not exist). Moreover, if $\alpha=0$, then $V\equiv 1$.

###### Proof of Theorem 4.2.

By performing the process similar to the parabolic case in the proof of
Theorem 1.1, we can show that after passing to a subsequence of $\delta\to 0$,

$U\to\;\text{ some }\;U_{0}\;\text{ in }\;L^{2}(\Omega_{1}),$

where $U_{0}$ is a weak solution of (4.12) ($U_{0}$ may be 0 in $\Omega_{1}$).
It is easy to see that the weak solution of (4.12) is unique. Thus, all above
convergences hold without passing to a subsequence. We next will show that
$U_{0}\not\equiv 0$ to ensure that $U_{0}$ is positive, namely, $U_{0}>0$ in
$\Omega_{1}$. Thus, we have $U_{0}=V$.

In what follows, it suffices to prove that $U_{0}\not\equiv 0$ due to the fact
$U\geq 0$. In the case of $\frac{\sigma}{\delta}\to 0$ as $\delta\to 0$, we
can construct a lower solution by defining

$\underline{U}(x)=\varepsilon e_{1}^{\delta}(x)$

where $\varepsilon>0$ is sufficiently small. For small $\delta$, it holds from
Lemma 4.1 that

$-\nabla\cdot(A(x)\nabla\underline{U})-\underline{U}(1-\underline{U})=\underline{U}(\lambda_{1}^{\delta}-1+\varepsilon
e_{1}^{\delta})\leq 0.$

Since $U\equiv 1$ is a upper solution, it follows from the lower and upper
solutions method that

$U\geq\underline{U}=\varepsilon e_{1}^{\delta},$

from which we have $U_{0}\geq\underline{U}>0$ in $\Omega_{1}$.

On the other hand, in the case of
$\frac{\sigma}{\delta}\to\alpha\in(0,\infty]$ as $\delta\to 0$, we will prove
$U_{0}>0$ by contradiction. If not, then $U_{0}\equiv 0$. Thus, let
$W=\frac{U}{||U||_{L^{2}(\Omega)}}>0$, and by (4.11), $W$ satisfies the
following elliptic problem

$\left\\{\begin{array}[]{llr}-\nabla\cdot(A(x)\nabla
W)=W(1-U),&\mbox{$x\in\Omega,$}\\\ W=0,&\mbox{$x\in\partial\Omega.$}&\\\
\end{array}\right.$

Multiplying both sides by $W$ and performing integration by parts, we have

$k\int_{\Omega_{1}}|\nabla W|^{2}dx+\int_{\Omega_{2}}A(x)\nabla W\cdot\nabla
Wdx+\int_{\Omega}W^{2}(U-1)dx=0,$ (4.13)

from which $\\{W\\}_{\delta>0}$ is bounded in $H^{1}_{0}(\Omega_{1})$. As a
result, after passing to a subsequence of $\delta$, $W\to$ some $W_{0}$ weakly
in $H^{1}_{0}(\Omega_{1})$ and strongly in $L^{2}(\Omega_{1})$. Moreover,
similar to (4.6), we get

$\begin{split}\int_{\Omega_{2}}W^{2}dx=&\int_{\Gamma}\int_{0}^{\delta}W^{2}(1+H(s)r+\kappa(s)r^{2})drds\\\
\leq&C(\delta^{2})\int_{\Gamma}\int_{0}^{\delta}W^{2}_{r}drds\\\
=&O(\frac{\delta^{2}}{\sigma})\to 0\;\text{ as }\;\delta\to 0.\end{split}$
(4.14)

Consequently, we have $||W_{0}||_{L^{2}(\Omega_{1})}=1$ because of the fact
$||W||_{L^{2}(\Omega)}=1$. Following the process as in the parabolic case in
the proof of Theorem 1.1, we can show that $W_{0}$ is a weak solution of

$-k\Delta W_{0}=W_{0},\quad x\in\Omega_{1},$

subject to the corresponding effective boundary conditions in Table 1 (with
$v$ replaced by $W_{0}$). This means that $\lambda_{1}=1$, which contradicts
the necessary condition for the existence of $V$.

∎

###  4.2 Main proof of Theorem 1.2

In the remaining part, our primary objective is to establish the proof of
Theorem 1.2, which unveils the infinite lifespan of each EBC as in Table 1.

###### Proof of Theorem 1.2.

If the initial value $u_{0}\equiv 0$, then it is easy to see that $u=v=0$ in
$\Omega_{1}\times(0,\infty)$, which implies this theorem trivially. Thus, we
always have the assumption of $u_{0}\not\equiv 0$.

Case 1. $\frac{\sigma}{\delta}\to 0$ as $\delta\to 0$, i.e., $\alpha=0$.

According to Theorem 4.1, as $\delta\to 0$, it holds that
$\lambda_{1}^{\delta}\to\lambda_{1}$, where $\lambda_{1}$ is the principal
eigenvalue of (4.3) with the corresponding effective boundary conditions
(second column in Table 1). By Lemma 4.1, it turns out that $\lambda_{1}=0$.
Thus, for a sufficiently small $\delta>0$, it is easy to see that
$\lambda_{1}^{\delta}<1$, implying that (4.12) has a unique positive steady
state solution $U$. Moreover, it follows from Theorem 4.2 and Remark 4.1 that
as $\delta\to 0$, $U\to V$ in $L^{2}(\Omega_{1})$ with $V\equiv 1$.

Subsequently, we consider the eigenvalue problem of

$\left\\{\begin{array}[]{llr}-\nabla\cdot(A(x)\nabla u)-u(1-U)=\eta
u,&\mbox{$x\in\Omega,$}\\\ u=0,&\mbox{$x\in\partial\Omega,$}&\\\
\end{array}\right.$ (4.15)

from which $(\eta_{1},u)=(0,U)$ is the principal eigenpair with $U>0$ in
$\Omega$. Let $\widetilde{u}=u-U$ and then $\widetilde{u}$ satisfies

$\left\\{\begin{array}[]{llr}\widetilde{u}_{t}-\nabla\cdot(A(x)\nabla\widetilde{u})=\widetilde{u}(1-u-U),&\mbox{$(x,t)\in\Omega\times(0,\infty),$}\\\
\widetilde{u}=u_{0}-U,&\mbox{$(x,t)\in\partial\Omega\times(0,\infty).$}&\\\
\end{array}\right.$ (4.16)

Multiplying (4.16) by $\widetilde{u}$ and integrating by parts in $x$, we have

$\begin{split}\frac{1}{2}\frac{d}{dt}\int_{\Omega}\widetilde{u}^{2}dx=&-\int_{\Omega}\nabla\widetilde{u}\cdot
A\nabla\widetilde{u}+\widetilde{u}^{2}(1-U)dx-\int_{\Omega}\widetilde{u}^{2}udx\\\
\leq&-\int_{\Omega}\widetilde{u}^{2}udx\\\ \leq&0,\end{split}$ (4.17)

which results from Lemma 2.1 and the variational characterization of
$\eta_{1}$ as

$0=\eta_{1}=\underset{0\not\equiv u\in
H^{1}_{0}(\Omega)}{\inf}\frac{\int_{\Omega}A\nabla u\cdot\nabla
u+u^{2}(U-1)dx}{\int_{\Omega}u^{2}dx}.$ (4.18)

∎

Consequently, given a fixed $T>0$, and for any $t\geq T$, (4.17) gives

$\begin{split}||u(\cdot,t)-U||_{L^{2}(\Omega_{1})}\leq&||u(\cdot,t)-U||_{L^{2}(\Omega)}\\\
\leq&||u(\cdot,T)-U||_{L^{2}(\Omega)}\\\
\leq&||u(\cdot,T)-v(\cdot,T)||_{L^{2}(\Omega_{1})}+||v(\cdot,T)-V||_{L^{2}(\Omega_{1})}\\\
&+||U-V||_{L^{2}(\Omega_{1})}+||u(\cdot,T)-U||_{L^{2}(\Omega_{2})}.\end{split}$
(4.19)

Since $\lambda_{1}=0$ and $u_{0}\geq 0,\not\equiv 0$, as $t\to\infty$, it
holds that

$v(x,t)\to V\equiv 1\;\text{ uniformly in }\;\overline{\Omega}_{1},$

which means that for any $\varepsilon>0$, there exists some large
$T_{\varepsilon}>0$ such that

$|v(x,t)-V(x)|<\varepsilon,\quad\forall x\in\overline{\Omega}_{1},$

provided $t\geq T_{\varepsilon}$. Thus, we have

$||v(\cdot,t)-V||_{L^{2}(\Omega_{1})}\leq|\Omega_{1}|^{1/2}\varepsilon$

for $t\geq T_{\varepsilon}$. Furthermore, for a sufficiently small $\delta>0$,
thanks to Theorem 1.1 and Lemma 2.1, we are led to

$\begin{split}||u(\cdot,T_{\varepsilon})-v(\cdot,T_{\varepsilon})||_{L^{2}(\Omega_{1})}\leq\varepsilon,\end{split}$
(4.20)

and

$\begin{split}||u(\cdot,T_{\varepsilon})-U||_{L^{2}(\Omega_{2})}\leq\varepsilon.\end{split}$
(4.21)

As a result, for any $t\geq T_{\varepsilon}$, due to Theorem 4.2 and above
estimates, (4.19) gives

$\begin{split}||u(\cdot,t)-U||_{L^{2}(\Omega_{1})}\leq&||u(\cdot,T_{\varepsilon})-v(\cdot,T_{\varepsilon})||_{L^{2}(\Omega_{1})}+||v(\cdot,T_{\varepsilon})-V||_{L^{2}(\Omega_{1})}\\\
&+||U-V||_{L^{2}(\Omega_{1})}+||u(\cdot,T)-U||_{L^{2}(\Omega_{2})}\\\
\leq&C\varepsilon.\end{split}$ (4.22)

Therefore, in view of (4.19)-(4.21), we have

$\begin{split}&\underset{t\in[T_{\varepsilon},\infty]}{\max}||u(\cdot,t)-v(\cdot,t)||_{L^{2}(\Omega_{1})}\\\
\leq&\underset{t\in[T_{\varepsilon},\infty]}{\max}||u(\cdot,t)-U||_{L^{2}(\Omega_{1})}+\underset{t\in[T_{\varepsilon},\infty]}{\max}||v(\cdot,t)-V||_{L^{2}(\Omega_{1})}+||U-V||_{L^{2}(\Omega_{1})}\\\
\leq&||u(\cdot,T_{\varepsilon})-U||_{L^{2}(\Omega_{1})}+\underset{t\in[T_{\varepsilon},\infty]}{\max}||v(\cdot,t)-V||_{L^{2}(\Omega_{1})}+||U-V||_{L^{2}(\Omega_{1})}\\\
\leq&C\varepsilon.\\\ \end{split}$ (4.23)

Thus, for a small $\delta>0$, it holds that

$\begin{split}\underset{0\leq
t\leq\infty}{\max}||u(\cdot,t)-v(\cdot,t)||_{L^{2}(\Omega_{1})}&\leq\underset{t\in[0,T_{\varepsilon}]}{\max}||u(\cdot,t)-v(\cdot,t)||_{L^{2}(\Omega_{1})}+\underset{t\in[T_{\varepsilon},\infty]}{\max}||u(\cdot,t)-v(\cdot,t)||_{L^{2}(\Omega_{1})}\\\
&\leq C\varepsilon,\end{split}$ (4.24)

which completes the proof of the first case.

Case 2. $\frac{\sigma}{\delta}\to\alpha\in(0,\infty)$ as $\delta\to 0$.

Thanks to Theorem 4.1 and Lemma 4.2, we have $\lambda_{1}>0$, where
$\lambda_{1}$ is the principal eigenvalue of (4.3) with the corresponding
boundary conditions (third column in Table 1).

If $\lambda_{1}<1$, then it holds that $\lambda_{1}^{\delta}<1$ for small
$\delta>0$, from which both $U$ and $V$ exist. Similar to the case as in the
proof of Case 1, we complete the proof of this theorem by deriving (4.24).

If $\lambda_{1}>1$, then it holds that $\lambda_{1}^{\delta}>1$ for small
$\delta>0$, from which neither $U$ nor $V$ exists, namely, $U=V=0$. For small
$\delta>0$, notice

$\begin{split}\frac{1}{2}\frac{d}{dt}\int_{\Omega}u^{2}dx=&-\int_{\Omega}\nabla
u\cdot A\nabla u+u^{2}dx-\int_{\Omega}u^{3}dx\\\
\leq&(-\lambda_{1}^{\delta}+1)\int_{\Omega}u^{2}dx\\\ \leq&0,\end{split}$
(4.25)

where we have used (4.2) and $u\geq 0$. Thus, similar to the case as in the
proof of Case 1, we get (4.24) that finishes this proof.

If $\lambda_{1}=1$, then we have two possible cases for the behavior of
$\lambda_{1}^{\delta}$. Suppose after passing to a subsequence of $\delta>0$,
we have $\lambda_{1}^{\delta}<1$ and $\lambda_{1}^{\delta}\to\lambda_{1}$, in
which case, $U$ exists along this subsequence but $V$ does not (i.e., $V\equiv
0$). Thus, the proof can be accomplished by the same argument as in Case 1.

On the other hand, there exists a subsequence of $\delta$, along which we have
$\lambda_{1}^{\delta}\geq 1$. As a result, neither $U$ nor $V$ exist along
this subsequence, namely, $U=V=0$. Moreover, direct computation yields (4.25),
and similar to the case as in the proof of Case 1, we prove this theorem by
obtaining (4.24).

Case 3. $\frac{\sigma}{\delta}\to\infty$ as $\delta\to 0$.

If $\lambda_{1}^{D}<1$, it follows from Theorem 4.1 that
$\lambda_{1}^{\delta}<1$ for small $\delta>0$, resulting in the existence of
$U$ and $V$. By using the method analogous to Case 1, the proof can be
finished.

If $\lambda_{1}^{D}>1$, we have $\lambda_{1}^{\delta}>1$ for small $\delta>0$,
from which $U=V\equiv 0$. By (4.25), we can prove this theorem by performing
the same argument as in Case 1.

If $\lambda_{1}^{D}=1$, we have $\lambda_{1}^{\delta}<1$ for small $\delta>0$,
leading to $U=V\equiv 0$. The proof is same as in Case 1.

Therefore, we complete the whole proof of this theorem.

## References

  * [1] H. Brezis, L.A. Caffarelli and A. Friedman  _Reinforcement problems for elliptic equations and variational inequalities_ ,  Ann. Mat. Pura Appl., 123 (1980), 219–246.
  * [2] H. Carslaw and J. Jaeger,  _Conduction of heat in solids_ ,  Reprint of the second edition, New York, 1988.
  * [3] X. Chen, J. He and X. Wang,  _Asymptotic propagation speeds of the Fisher-KPP equation with an effective boundary condition on a road_ ,  Arch. Ration. Mech. Anal., 247 (2023), 0003-9527.
  * [4] X. Chen, C. Pond and X. Wang,  _Effective boundary conditions resulting from anisotropic and optimally aligned coatings: the two dimensional case_ ,  Arch. Ration. Mech. Anal., 206 (2012), 911–951.
  * [5] X. Geng,  _Effective Boundary Conditions Arising from the Heat Equation with Three-dimensional Interior Inclusion_ ,  Comm. Pure Appl. Anal., 22 (2023), 1394-1419.
  * [6] ———–,  _Effective Boundary Conditions for Heat Equation Arising from Anisotropic and Optimally Aligned Coatings in Three Dimensions_ ,  preprint, (2023).
  * [7] D. Gilbarg and N. Trudinger,  _Elliptic partial differential equations of second order_ ,  Reprint of the 1998 edition, Springer-Verlag, Berlin, 2001.
  * [8] T. Hillen, K. Painter and M. Winkler,  _Anisotropic diffusion in oriented environments can lead to singularity formation_ ,  European Journal of Applied Mathematics, 24 (2013), 371-413.
  * [9] H. Li, J. Li and X. Wang,  _Error estimates and lifespan of effective boundary conditions for 2-dimensional optimally aligned coatings_ ,  J. Differential Equations, 303 (2022), 1-41.
  * [10] H. Li and J. Li,  _Lifespan of effective boundary conditions for the heat equation_ ,  J. Math. Phys., 63 (2022), 0022-2488.
  * [11] H. Li and X. Wang,  _Using effective boundary conditions to model fast diffusion on a road in a large field_ ,  Nonlinearity, 30 (2017), 3853–3894.
  * [12] H. Li, X. Wang and Y. Wu,  _Effective boundary conditions and lifespan of the logistic diffusion equation on a coated body_ ,  J. Differential Equations, 257 (2014), 3640–3668.
  * [13] J. Li,  _Asymptotic behavior of solutions to elliptic equations in a coated body_ ,  Comm. Pure App. Anal., 8 (2009), 1251–1267.
  * [14] J. Li, S. Rosencrans, X. Wang and K. Zhang,  _Asymptotic analysis of a Dirichlet problem for the heat equation on a coated body_ ,  Proc. Amer. Math. Soc., 137 (2009), 1711-1721.
  * [15] J. Li, L. Su, X. Wang and Y. Wang,  _Bulk-surface coupling: derivation of two models_ ,  J. Differential Equations, 289 (2021), 1-34.
  * [16] J. Li and K. Zhang,  _Reinforcement of the Poisson equation by a thin layer_ ,  Math. Models Methods Appl. Sci., 21 (2011), 1153–1192.
  * [17] C.V. Pao,  _Nonlinear parabolic and elliptic equations_ ,  Springer Science, 2012.
  * [18] C. Pond,  _Effective Boundary Conditions on Coated Bodies_ ,  Ph.D. thesis, Tulane University, 2011.
  * [19] S. Rosencrans and X. Wang,  _Suppression of the Dirichlet eigenvalues of a coated body_ ,  SIAM J. Appl. Math., 66 (2006), 1895–1916; Corrigendum, SIAM J. Appl. Math., 68 (2008), p1202.
  * [20] X. Wang,  _Effective boundary conditions of diffusion equations on domains containing thin layers (in Chinese)_ ,  Sci. Sin. Math., 46 (2016), 709-724.


# Signatures of Quantum Chaos and fermionization in the incoherent transport
of bosonic carriers in the Bose-Hubbard chain

P. S. Muraev  Kirensky Institute of Physics, Federal Research Centre KSC SB
RAS, 660036, Krasnoyarsk, Russia  School of Engineering Physics and Radio
Electronics, Siberian Federal University, 660041, Krasnoyarsk, Russia  IRC
SQC, Siberian Federal University, 660041, Krasnoyarsk, Russia  D. N. Maksimov
Kirensky Institute of Physics, Federal Research Centre KSC SB RAS, 660036,
Krasnoyarsk, Russia  IRC SQC, Siberian Federal University, 660041,
Krasnoyarsk, Russia  A. R. Kolovsky  Kirensky Institute of Physics, Federal
Research Centre KSC SB RAS, 660036, Krasnoyarsk, Russia  School of Engineering
Physics and Radio Electronics, Siberian Federal University, 660041,
Krasnoyarsk, Russia

###### Abstract

We analyse the stationary current of Bose particles across the Bose-Hubbard
chain connected to a battery, focusing on the effect of inter-particle
interactions. It is shown that the current magnitude drastically decreases as
the strength of inter-particle interactions exceeds the critical value which
marks the transition to quantum chaos in the Bose-Hubbard Hamiltonian. We
found that this transition is well reflected in the non-equilibrium many-body
density matrix of the system. Namely, the level-spacing distribution for
eigenvalues of the density matrix changes from Poisson to Wigner-Dyson
distributions. With the further increase of the interaction strength, the
Wigner-Dyson spectrum statistics changes back to the Poisson statistics which
now marks fermionization of the bosonic particles. With respect to the
stationary current, this leads to the counter-intuitive dependence of the
current magnitude on the particle number.

1. Recently, we have seen a surge of interest to quantum transport in one-dimensional systems coupled at their edges to particle reservoirs. Following [1] we refer to these systems as boundary driven systems. Particularly, one very important example of the boundary driven system is the so-called open Bose-Hubbard (BH) model or BH chain [2, 3, 4, 5]. This system can be realised experimentally by using different physical platforms including superconducting circuits [6, 7], photonic crystals [8, 9], and cold Bose atoms in optical lattices [10]. The central question to be addressed with the open BH chain, both theoretically and experimentally, is the stationary current of Bose particles across the chain and its dependence on the strength of inter-particle interactions. It is known that properties of the closed/conservative BH system crucially depend on the ratio of the hopping matrix element $J$ and the interaction constant $U$ which are two of the four parameters of the BH Hamiltonian, the other being the chain length $L$ and the particle number $N$. For example, for integer $N/L$ the ground state of the system shows the quantum phase transition between the super-fluid state for $J\gg U$ and the Mott-insulator state in the opposite limit [11]. As for the excited states they show a qualitative change from the regular to the chaotic [12, 13]. Thus, one may expect that the stationary current in the open BH chain may also crucially depend on the interaction constant.

Up to now the above problem has been analysed only by using the
pseudoclassical approach which approximates the quantum dynamics by that of
the classical counterpart of the BH model [4]. This is because the exact
quantum treatment, both numerical and analytic, is complicated by the fact
that the open BH model does not conserve the particle number. Thus, the system
dynamics takes place in the extended Hilbert space spanned by the direct sum
of the subspaces with given number of particles. Taking into account that the
transition to quantum chaos in the closed BH model depends not only on the
ratio $U/J$ but also on the ratio $N/L$, it is very problematic to track this
transition to in the transport properties of the open BH model.

In the present work we introduce a boundary driven BH model which shows a
behaviour similar to the standard open BH model but conserves the number of
particles. This allows us to separate the dependencies on the interaction
constant $U$ and the particle number $N$, thus, relating the obtained results
to the transition to chaos in the conservative BH model.

2. We consider the BH chain of the length $L$ with incoherent coupling between the first and the $L_{{\rm th}}$ sites. The coupling is described by the the following Lindblad operators,

$\begin{array}[]{cc}{\cal
L}_{1}(\widehat{R})&=\widehat{V}^{\dagger}\widehat{V}\widehat{R}+\widehat{R}\widehat{V}^{\dagger}\widehat{V}-2\widehat{V}\widehat{R}\widehat{V}^{\dagger}\;,\\\
\\\ {\cal
L}_{2}(\widehat{R})&=\widehat{V}\widehat{V}^{\dagger}\widehat{R}+\widehat{R}\widehat{V}\hat{V}^{\dagger}-2\widehat{V}^{\dagger}\widehat{R}\hat{V}\;,\end{array}$
(1)

where $\widehat{V}=\hat{a}_{1}^{\dagger}\hat{a}_{L}$. Thus, the master
equation for the carriers density matrix $\widehat{R}$ reads

$\frac{\partial\widehat{R}}{\partial
t}=-i[\widehat{H},\widehat{R}]-\Gamma_{1}{\cal
L}_{1}(\widehat{R})-\Gamma_{2}{\cal L}_{2}(\widehat{R})\;,$ (2)

where

$\widehat{H}=-\dfrac{J}{2}\sum^{L-1}_{\ell=1}\left(\hat{a}^{\dagger}_{\ell+1}\hat{a}_{\ell}+\mathrm{h.c.}\right)+\dfrac{U}{2}\sum^{L}_{\ell=1}\hat{n}_{\ell}(\hat{n}_{\ell}-1)$
(3)

is the Bose-Hubbard Hamiltonian. It is easy to see that the Lindblad operator
${\cal L}_{1}(\widehat{R})$ induces the incoherent transport of the carriers
from the last to the first sites, while the operator ${\cal
L}_{2}(\widehat{R})$ is responsible for the incoherent transport in the
reverse direction. If the rates $\Gamma_{1}\neq\Gamma_{2}$, there is a non-
zero current in the clockwise or counterclockwise direction depending on the
inequality relationship between the two relaxation constants. We notice that,
by an analogy with electronic devices, the introduced Lindblad operators mimic
the effect of a battery which induces direct current in the electric circuits.

3. We are interested in the stationary current $I={\rm Tr}[\widehat{I}\widehat{R}]$ where $\widehat{R}=\widehat{R}(t\rightarrow\infty)$ is now the steady-state density matrix and $\widehat{I}$ is the current operator,

$\widehat{I}=\dfrac{J}{2i}\sum^{L-1}_{\ell=1}\left(\hat{a}^{\dagger}_{\ell+1}\hat{a}_{\ell}-h.c.\right)\;.$
(4)

First of all, we notice that if $\Gamma_{1}=\Gamma_{2}$ the steady-state
density matrix is proportional to the identity matrix, namely,
$\widehat{R}=\widehat{1}/{\cal N}$, where ${\cal N}$ is the dimension of the
Hilbert space. In what follows we focus on the liner response regime where
$\Gamma_{1}=\Gamma+\Delta\Gamma/2$, $\Gamma_{2}=\Gamma-\Delta\Gamma/2$, and
$\Delta\Gamma\ll\Gamma$. Thus, we have

$\widehat{R}=\frac{\widehat{1}}{{\cal N}}+\Delta\Gamma\widetilde{R}\;,$ (5)

where ${\rm Tr}[\widetilde{R}]=0$. Substituting the Ansatz (5) into the master
equation we obtain

$-i[\widehat{H},\widetilde{R}]-\Gamma\left[{\cal L}_{1}(\widetilde{R})+{\cal
L}_{2}(\widetilde{R})\right]-\frac{2(\hat{n}_{L}-\hat{n}_{1})}{{\cal
N}}=O(\Delta\Gamma)\;.$ (6)

In the limit $\Delta\Gamma\rightarrow 0$ Eq. (6) transforms into the algebraic
equation for the elements of the unknown matrix $\widetilde{R}$. In our
numerical approach, however, we do not solve this algebraic equation but
evolve the density matrix $\widehat{R}(t)$ according to the master equation
(2) and use Eq. (6) to check that we reached the true steady state. We found
this method to be more efficient than the straightforward solution of the
algebraic equation.

Figure 1: Absolute values of the matrix elements of the matrix $\widetilde{R}$
in the basis of the current operator. The system parameters are $L=6$, $N=3$
(the Hilbert space dimension ${\cal N}=56$), $J=1$, $\Gamma=0.04$, and $U=0$
(left) and $U=1$ (right). The upper limit of the colour axis is $0.2$.

Since our primary goal is the stationary current across the chain, we shall
analyse the matrix $\widetilde{R}$ in the basis of the eigenstates of the
current operator,

$\widehat{I}=\sum_{j=1}^{\cal N}\sigma_{j}|\Phi_{j}\rangle\langle\Phi_{j}|\;.$
(7)

Two examples of the matrix $\widetilde{R}$ in this basis are given in Fig. 1
for $U=0$, left panel, and $U=J$, right panel. A qualitative difference
between these two cases is clearly visible from the plot. In the next
paragraph we quantify this difference. We conclude the present paragraph by
displaying the commutation relation between the current operator and the BH
Hamiltonian for $U=0$,

$-i[\widehat{H},\widehat{I}]-\frac{(\hat{n}_{L}-\hat{n}_{1})}{2}=0\;,$ (8)

which we shall use later on.

Figure 2: The scaled eigenvalues of the matrix $\widetilde{R}$ for $U=0$, red
line, and $U=1$, magenta line, compared to the eigenvalues of the current
operator, blue line. Parameters are $L=8$, $N=4$ (${\cal N}=330$), $J=1$, and
$\Gamma=0.04$.

4. Knowing that the BH Hamiltonian (3) exhibits transition to quantum chaos as $U$ is increased, we expect a similar transition for the non-equilibrium density matrix. This expectation is supported by the visual analysis of the matrices depicted in Fig. 1 and results of the relevant studies on the boundary driven spin chains [14, 15]. Following Ref. [14, 15] we consider the spectrum and eigenstates of the non-equilibrium density matrix,

$\widetilde{R}=\sum_{j=1}^{\cal
N}\lambda_{j}|\Psi_{j}\rangle\langle\Psi_{j}|\;.$ (9)

In what follows we restrict ourselves by the case $\Gamma\ll J$. Then, if
$U=0$, the states $|\Psi_{j}\rangle$ practically coincide with the eigenstates
of the current operator $|\Phi_{j}\rangle$, while the eigenvalues are related
to each other as

$\lambda_{j}\approx\sigma_{j}/4{\cal N}\;.$ (10)

To see that let us scale the density matrix $\widetilde{R}\rightarrow 4{\cal
N}\widetilde{R}$ and set $\Delta\Gamma=0$ in Eq. (6). Then the obtained
algebraic equation differs from the commutation relation Eq. (8) by a small
term $\sim\Gamma$ which can be taken into account perturbatively. As expected,
application of the perturbation theory removes degeneracies of the eigenvalues
of the current operator, see the lower inset in Fig. 2.

Figure 3: The integrated level spacing distribution compared to the integrated
Poisson and GUE Wigner-Dyson distributions, $L=10$, $N=5$ (${\cal N}=2002$),
and $\Gamma=0.04$. For the statistical analysis we took 60 percents of
eigenvalues from the central part of the spectrum.

The magenta staircase curve in Fig. 2 depicts the case $U\neq 0$. Here we see
that the width of the spectrum increases as $U$ is increased. The main
difference is, however, in the spectrum statistics. Figure 3 shows the
distribution of the the scaled spacings
$s=(\lambda_{j+1}-\lambda_{j})f(\lambda_{j})$, $f(\lambda)$ being the mean
density of states, as compared to the Poisson and Wigner-Dyson distributions
[16]. For $U=J$ an excellent agreement with the Wigner-Dyson distribution,
which is the hallmark of quantum chaos [17, 18], is noticed.

5. Next we analyse the stationary current $I$ across the chain,

$I={\rm Tr}[\widehat{I}\widehat{R}]=\sum_{j=1}^{\cal
N}\lambda_{j}\langle\Psi_{j}|\widehat{I}|\Psi_{j}\rangle\equiv\sum_{j=1}^{\cal
N}\lambda_{j}I_{j}\;.$ (11)

In the case of vanishing inter-particle interactions one derives by using Eq.
(10) the following semi-analytic equation,

$I=4JN^{2}\Delta\Gamma\int_{0}^{1}\sigma^{2}(x)dx\;,$ (12)

where $x=j/{\cal N}$ and $\sigma(x)$ is the inverse function to the integrated
density of states of the current operator, which interpolates the blue and red
lines in Fig. 2. Thus, as it is intuitively expected, for $U=0$ the total
current increases with the number of particles in the system.

The case $U\neq 0$ is more subtle. Figure 4 shows numerically obtained
dependence of the stationary current on the interaction constant $U$ for $L=6$
and different number of particles $N$. One clearly identifies in this figure
the critical $U_{cr}=U_{cr}(\bar{n})$, $\bar{n}=N/L$, above which the current
drastically decreases. This critical interaction marks the crossover from the
Poisson to Wigner-Dyson spectrum statistics for the non-equilubrium density
matrix $\widetilde{R}$. An unexpected result is that for $U\gg U_{cr}$ the
current decreases with the number of particles. Furthermore, we find that for
these large $U$ the spectrum statistics is again Poissonian.

Figure 4: The total current for $U=0$ as function of the interaction constant
$U$ for three different values of the particle number $N=2,3,4$ in the linear
and logarithmic scales, $L=10$ and $\Gamma=0.04$

6. We relate the observed change of the spectrum statistics and the counter-intuitive dependence of the current on the particle number to the interaction-induced localization of the eigenstates $|\Psi_{j}\rangle$ and the fermionization of the strongly interacting Bose particles [19]. Indeed, the obvious consequence of the eigenstate localization is that the mean $I_{j}=\langle\Psi_{j}|\widehat{I}|\Psi_{j}\rangle$ tends to zero and is strictly zero if all bosons occupies a single site of the chain. Figure 5 shows the quantiles $I_{j}$ for $(L,N)=(6,3)$ and $U=0,0.5,10$. It is seen that the fraction of the delocalized states which support the current decreases in favor of the localised states for which $I_{j}\approx 0$. For example, in Fig. 5 (c) the states corresponding to the minimal and maximal eigenvalues are the Fock states $|0,0,0,0,0,3\rangle$ and $|3,0,0,0,0,0\rangle$, respectively. Along with the localised and partially localised states one can see in Fig. 5 (c) a number of the delocalized states. A closer inspection of these states shows that they are a superposition of the Fock states where occupation numbers of the chain sites are either zero or unity. Since this subspace of the Hilbert space is the Hilbert space of the hard-core bosons, we conclude that the residual conductivity of the system at large $U$ is mainly due to the hard-core bosons. As known, the spectral and transport properties of the hard-core bosons are similar to those of the non-interacting fermions and, hence, they can support ballistic transport for arbitrary large $U$ if $N/L<1$. We also mention that the appearance of the localised states and the integrable ‘hard-core boson states’ is consistent with the observed change of the level-spacing distribution from the Wigner-Dyson distribution back to the Poisson distribution in the limit of large $U$.

Figure 5: The quantities $I_{j}=\langle\Psi_{j}|\widehat{I}|\Psi_{j}\rangle$
for $U=0$ (a) 0.5 (b), and 10 (c). The system parameters are $L=6$, $N=3$, and
$\Gamma=0.04$.

7. In summary, we introduced the model for quantum transport of Bose particles across the Bose-Hubbard chain which conserves the number of particles in the chain. Similar to the standard transport model where Bose-Hubbard chain connects two particle reservoirs with different chemical potentials and where the number of particles is not conserved, the introduced model shows different transport regimes depending on the ratio between the tunnelling and interaction constants in the Bose-Hubbard Hamiltonian (3). Namely, for $U\sim J$ the stationary current of the Bose particles is drastically suppressed as compared to the case $U=0$. In our previous publication [4] we explain this effect by transition to chaotic dynamics of the classical counterpart of the system. In this work we use the genuine quantum approach where the object of interest is the non-equilibrium many-body density matrix of the bosonic carriers in the chain. Up to the best of our knowledge this is the first example where the non-equilibrium density matrix is calculated/analysed for the non-integrable bosonic system. (For transport properties of the fermionic and spin systems we refer the reader to the recent review [20].) We found the spectrum of this matrix exhibits a transition from a regular spectrum for $U\ll J$, which obeys the Poisson statistics, to an irregular spectrum for $U\sim J$, which obeys the Wigner-Dyson statistics. In this sense we confirm the conjecture of Ref. [4] that the drastic reduction of the current for $U\sim J$ is due to the transition to quantum chaos.

Within the framework of the introduced model we also observed a new effect, –
the residual conductivity due to fermionization of the Bose particles. We
notice that this is a pure quantum effect which cannot be addressed by using
the classical (mean-field) or pseudoclassical (truncated Wigner function)
approaches.

Acknowledgment. This work has been supported by Russian Science Foundation
through Grant No. N19-12-00167.

## References

  * Landi  _et al._ [2022] G. T. Landi, D. Poletti, and G. Schaller, Nonequilibrium boundary-driven quantum systems: Models, methods, and properties, Reviews of Modern Physics 94, 045006 (2022).
  * Ivanov  _et al._ [2013] A. Ivanov, G. Kordas, A. Komnik, and S. Wimberger, Bosonic transport through a chain of quantum dots, The European Physical Journal B 86, 10.1140/epjb/e2013-40417-4 (2013).
  * Kolovsky  _et al._ [2018] A. R. Kolovsky, Z. Denis, and S. Wimberger, Landauer-Büttiker equation for bosonic carriers, Physical Review A 98, 043623 (2018).
  * Bychek  _et al._ [2020] A. A. Bychek, P. S. Muraev, D. N. Maksimov, and A. R. Kolovsky, Open Bose-Hubbard chain: Pseudoclassical approach, Physical Review E 101, 012208 (2020).
  * Muraev  _et al._ [2022] P. S. Muraev, D. N. Maksimov, and A. R. Kolovsky, Resonant transport of bosonic carriers through a quantum device, Physical Review A 105, 013307 (2022).
  * Fitzpatrick  _et al._ [2017] M. Fitzpatrick, N. M. Sundaresan, A. C. Li, J. Koch, and A. A. Houck, Observation of a Dissipative Phase Transition in a One-Dimensional Circuit QED Lattice, Physical Review X 7, 011016 (2017).
  * Fedorov  _et al._ [2021] G. Fedorov, S. Remizov, D. Shapiro, W. Pogosov, E. Egorova, I. Tsitsilin, M. Andronik, A. Dobronosova, I. Rodionov, O. Astafiev, and A. Ustinov, Photon Transport in a Bose-Hubbard Chain of Superconducting Artificial Atoms, Physical Review Letters 126, 180503 (2021).
  * Szameit  _et al._ [2005] A. Szameit, D. Blömer, J. Burghoff, T. Schreiber, T. Pertsch, S. Nolte, A. Tünnermann, and F. Lederer, Discrete nonlinear localization in femtosecond laser written waveguides in fused silica, Optics Express 13, 10552 (2005).
  * Cáceres-Aravena  _et al._ [2022] G. Cáceres-Aravena, D. Guzmán-Silva, I. Salinas, and R. A. Vicencio, Controlled Transport Based on Multiorbital Aharonov-Bohm Photonic Caging, Physical Review Letters 128, 256602 (2022).
  * Labouvie  _et al._ [2016] R. Labouvie, B. Santra, S. Heun, and H. Ott, Bistability in a Driven-Dissipative Superfluid, Physical Review Letters 116, 235302 (2016).
  * Greiner  _et al._ [2002] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms, Nature 415, 39 (2002).
  * Kolovsky and Buchleitner [2004] A. R. Kolovsky and A. Buchleitner, Quantum chaos in the Bose-Hubbard model, Europhysics Letters (EPL) 68, 632 (2004).
  * Kolovsky [2016] A. R. Kolovsky, Bose–Hubbard Hamiltonian: Quantum chaos approach, International Journal of Modern Physics B 30, 1630009 (2016).
  * Prosen and Žnidarič [2013] T. Prosen and M. Žnidarič, Eigenvalue Statistics as an Indicator of Integrability of Nonequilibrium Density Operators, Physical Review Letters 111, 124101 (2013).
  * Žnidarič [2013] M. Žnidarič, Magnetization transport in spin ladders and next-nearest-neighbor chains, Physical Review B 88, 205135 (2013).
  * Wigner [1951] E. P. Wigner, On the statistical distribution of the widths and spacings of nuclear resonance levels, Mathematical Proceedings of the Cambridge Philosophical Society 47, 790 (1951).
  * Haake [1991] F. Haake, Quantum Signatures of Chaos, in  _Quantum Coherence in Mesoscopic Systems_, edited by B. Kramer (Springer US, Boston, MA, 1991) pp. 583–595.
  * Stöckmann [1999] H.-J. Stöckmann,  _Quantum Chaos: An Introduction_ (Cambridge University Press, Cambridge, 1999).
  * Paredes  _et al._ [2004] B. Paredes, A. Widera, V. Murg, O. Mandel, S. Fölling, I. Cirac, G. V. Shlyapnikov, T. W. Hänsch, and I. Bloch, Tonks–Girardeau gas of ultracold atoms in an optical lattice, Nature 429, 277 (2004).
  * Bertini  _et al._ [2021] B. Bertini, F. Heidrich-Meisner, C. Karrasch, T. Prosen, R. Steinigeweg, and M. Žnidarič, Finite-temperature transport in one-dimensional quantum lattice models, Reviews of Modern Physics 93, 025003 (2021).