# Stochastic completeness and $L^{1}$-Liouville property for second-order
elliptic operators

Debdip Ganguly  Department of Mathematics,
Indian Institute of Technology Delhi,
IIT Campus, Hauz Khas, Delhi,
New Delhi 110016, India.
E-mail address<EMAIL_ADDRESS>,  Yehuda Pinchover  Department of
Mathematics,
Technion - Israel Institute of Technology,
Haifa 32000, Israel.
E-mail address<EMAIL_ADDRESS>and  Prasun Roychowdhury  Department of
Mathematics,
Indian Institute of Science Education and Research,
Dr. Homi Bhabha Road, Pashan,
Pune 411008, India.
E-mail address<EMAIL_ADDRESS>
###### Abstract.

Let $P$ be a linear, second-order, elliptic operator with real coefficients
defined on a noncompact Riemannian manifold $M$ and satisfies $P1=0$ in $M$.
Assume further that $P$ admits a minimal positive Green function in $M$. We
prove that there exists a smooth positive function $\rho$ defined on $M$ such
that $M$ is stochastically incomplete with respect to the operator
$P_{\rho}:=\rho\,P$, that is,

$\int_{M}k_{P_{\rho}}^{M}(x,y,t)\ {\rm d}y<1\qquad\forall(x,t)\in
M\times(0,\infty),$

where $k_{P_{\rho}}^{M}$ denotes the minimal positive heat kernel associated
with $P_{\rho}$. Moreover, $M$ is $L^{1}$-Liouville with respect to $P_{\rho}$
if and only if $M$ is $L^{1}$-Liouville with respect to $P$. In addition, we
study the interplay between stochastic completeness and the $L^{1}$-Liouville
property of the skew product of two second-order elliptic operators.

###### Key words and phrases:

Green function, $L^{1}$-Liouville, optimal Hardy, stochastically incomplete.

###### 2010 Mathematics Subject Classification:

Primary 31C35, ; Secondary 35J08, 35A08, 35B09, 58G03.

##  1\. Introduction

Let $M$ be a smooth, noncompact, connected Riemannian manifold of dimension
$n$. Let $P$ be a linear, second-order, elliptic operator with real
coefficients defined on $M$, and satisfying $P1=0$ in $M$.

Denote by $\mathcal{C}_{P}(M)$ the cone of positive solutions of the equation
$Pu=0$ in $M$. The generalized principal eigenvalue of the operator $P$ is
defined by

$\lambda_{0}=\lambda_{0}(P,M):=\sup\\{\lambda\in\mathbb{R}\;\mid\;\mathcal{C}_{P-\lambda}(M)\neq\emptyset\\}.$

We say that $P$ is nonnegative in $M$ (and denote it by $P\geq 0$ in $M$) if
$\lambda_{0}:=\lambda_{0}(P,M)\geq 0.$ Since $P1=0$ in $M$, it follows that
$\lambda_{0}\geq 0$, that is, $P\geq 0$ in $M.$

Consider the parabolic operator

$Lu:=\frac{\partial u}{\partial t}+Pu\qquad(x,t)\in M\times(0,\infty),$ (1.1)

and let $k^{M}_{P}(x,y,t)$ be the _minimal positive heat kernel_ of the
parabolic operator $L$ on the manifold $M.$ By definition, for a fixed $y\in
M,$ the function $(x,t)\mapsto k^{M}_{P}(x,y,t)$ is a _minimal_ positive
solution of the equation

$Lu=0\quad\mbox{in}\ M\times(0,\infty),$ (1.2)

subject to the initial data $\delta_{y},$ the Dirac distribution at $y\in M$.

######  Definition 1.1.

Suppose that $\lambda_{0}=\lambda_{0}(P,M)\geq 0$, and let $k_{P}^{M}$ be the
corresponding heat kernel. We say that the operator $P$ is _subcritical_
(respectively, _critical_) in $M$ if for some $x\not=y$, (and therefore for
any $x\not=y$), $x,y\in M$, we have

$\int_{0}^{\infty}k_{P}^{M}(x,y,t)\,\mathrm{d}t<\infty,\quad\left(\mbox{respectively,
}\int_{0}^{\infty}k_{P}^{M}(x,y,t)\,\mathrm{d}t=\infty\right).$ (1.3)

If $P$ is subcritical in $M$, then

$G_{P}^{M}(x,y):=\int_{0}^{\infty}k_{P}^{M}(x,y,t)\,\mathrm{d}t\qquad x,y\in
M$ (1.4)

is called the minimal positive Green function of the operator $P$ in $M$.

Next we introduce the notion of stochastically completeness.

######  Definition 1.2.

Let $(M,g)$ be a connected and noncompact Riemannian manifold of dimension
$N.$ Let $P$ be an elliptic operator satisfying $P1=0$ in $M$. Then $M$ is
said to be _stochastically complete_ (respectively, _stochastically
incomplete_) with respect to $P$ if for some $(x,t)\in M\times(0,\infty)$ (and
therefore for all $(x,t)\in M\times(0,\infty)$) we have

$\int_{M}k_{P}^{M}(x,y,t)\,\mathrm{d}y=1,\quad\left(\mbox{respectively,
}\int_{M}k_{P}^{M}(x,y,t)\,\mathrm{d}y<1\right).$

######  Definition 1.3.

Suppose that the manifold $M$ and the operator $P$ satisfy the above
assumptions. We say that the $L^{1}$-Liouville property holds on $M$ with
respect to $P$ (in short, $M$ is $L^{1}$-Liouville with respect to $P$) if
every nonnegative $L^{1}$-supersolution of the operator $P$ in $M$ is the
constant function.

######  Remark 1.4.

Recall that if $u$ is a positive solution of the equation $Pu=0$ in $M$, then
the generalized maximum principle implies that for any $t>0$, either

$\int_{M}k_{P}^{M}(x,y,t)u(y)\,\mathrm{d}y=u(x)\quad\mbox{or
}\int_{M}k_{P}^{M}(x,y,t)u(y)\,\mathrm{d}y<u(x)\quad\forall x\in M.$

Equivalently, for any $\lambda>0$ either

$\lambda\int_{M}G_{P+\lambda}^{M}(x,y)u(y)\,\mathrm{d}y=u(x)\quad\mbox{or
}\lambda\int_{M}G_{P+\lambda}^{M}(x,y)u(y)\,\mathrm{d}y<u(x)\quad\forall x\in
M.$

Moreover, $P$ is critical in $M$ if and only if $P$ admits a unique (up to a
multiplicative constant) positive supersolution (see [18] and references
therein). Therefore, if $P$ is critical in $M$ and $P1=0$, then $M$ is
stochastically complete and $L^{1}$-Liouville with respect to $P.$

Throughout the paper, unless otherwise stated, we assume that $P$ is
subcritical in $M$. In other words, we assume that $P$ admits a (unique)
minimal positive Green function $G_{P}^{M}$ on $M$. In this paper we are
mainly interested in the following question:

######  Problem 1.5.

Given an operator $P$ of the form (2.1) on $M$, construct an operator
$P_{\rho}:=\rho\,P$ for some positive smooth function $\rho$ on $M$ such that
$M$ is stochastically incomplete with respect to $P_{\rho}$, and such that $M$
is $L^{1}$-Liouville with respect to $P_{\rho}$ if and only if $M$ is
$L^{1}$-Liouville with respect to $P$.

We provide an affirmative answer to the above question. Our proof rests on
certain constructions of positive supersolutions and related Hardy weights [3,
4], and a generalization of the Omori-Yau maximum principle (see [1, 14]).

Stochastically completeness and the $L^{1}$-Liouville property on manifolds
have been studied extensively in the recent years. We mention a few of them
without a claim of completeness (see [1, 6, 9, 13, 14, 16]), a more
comprehensive reference can be found in [5] and references therein. It has
been shown that certain natural geometric assumptions on manifolds ensure the
validity of these properties on a complete manifold, e.g., volume growth,
bound on the curvature, etc. These questions were well understood in case of
the Laplace-Beltrami operator $\Delta_{g}$ using potential theoretic arguments
and exploiting some underlying geometric conditions. We mention also [15],
where the authors study properties like parabolicity, stochastic completeness
and $L^{1}$-Liouville property in a manifold with boundary, under the
Dirichlet boundary condition. In [8] the authors studied the parabolicity of a
manifold with respect to the Neumann boundary conditions. Moreover, recently
stochastic completeness was characterized in [7] in terms of certain
properties of nonlinear evolution equations of fast diffusion. Furthermore, in
[17] the author constructed an example of an operator $P$ defined on a
complete Riemannian manifold $M$, which does not admit invariant positive
solution at the bottom of the spectrum, providing a counter example to
Stroock’s conjecture (see also [5, 20]). It turns out that the constructed
manifold $M$ is stochastically incomplete and not $L^{1}$-Liouville with
respect to $P$. The present article is partly motivated also from the
construction of the example in [17], aiming to find a unified approach to
construct stochastically incomplete manifold with respect to a second-order
elliptic operator.

All the aforementioned results are the primary motivation of the study in this
paper. Indeed, we study the two notions defined in definitions 1.2 and 1.3 for
general second-order elliptic operators $P$ satisfying $P1=0$ in $M$. The rest
of the paper is organized as follows. In Section 2, we state the assumptions
on the operator $P$, and formulate our main result in Theorem 2.3. Section 3
is devoted to the proof of a few key lemmas which provide equivalent
conditions for stochastic completeness and the validity of the
$L^{1}$-Liouville property for second-order elliptic operator $P$. In Section
4, we prove our main result (Theorem 2.3). Section 5 is devoted to the study
of the $L^{1}$-Liouville property and the stochastic (in)completeness for skew
product operators. Finally, in Section 6 we give an alternative proof of
Theorem 2.3.

##  2\. Preliminaries and the main result

This section is devoted to the statement of the main result of the paper.
Before going further we must introduce some technical assumptions and few
definitions.

Let $M$ be a smooth, noncompact, and connected manifold of dimension $n$,
where $n\geq 2$. We consider a second-order elliptic operator $P$ defined on
$M$ with real coefficients which (in any coordinate system
$(U;x_{1},\ldots,x_{n})$) is of the form

$Pu=-\sum_{i,j=1}^{n}a^{ij}(x)\partial_{i}\partial_{j}u+b(x)\cdot\nabla u.$
(2.1)

We assume that for every $x\in M$ the matrix $A(x):=\big{[}a^{ij}(x)\big{]}$
is symmetric and that the real quadratic form

$\xi\cdot A(x)\xi:=\sum_{i,j=1}^{n}\xi_{i}a^{ij}(x)\xi_{j}\qquad\text{for
}\xi\in\mathbb{R}^{n}$

is positive definite. Throughout the paper it is assumed that the coefficients
of $P$ are Hölder continuous in $M$. By a solution $v$ of the equation $Pu=0$
in $M$, we mean that $v\in C^{2}(M)$ and satisfies the equation pointwise. If
$P$ is of the form (2.1), we obviously have

$P1=0.$ (2.2)

Whenever we consider the adjoint operator we assume that $P^{*}$, the formal
adjoint operator of $P$, also has Hölder continuous coefficients. We say that
$P$ is symmetric if $P=P^{*}$.

######  Remark 2.1.

It is well known that any smooth manifold may be considered as an embedded
manifold in $\mathbb{R}^{N}$, for some $N\geq n$. Hence, the Euclidean metric
and the Lebesgue measure on $\mathbb{R}^{N}$ induce a Riemannian metric and a
measure on $M$, which will be considered throughout the paper in case the
given manifold $M$ is not a priori Riemannian. Another possible structure is
induced by the principal part of the operator $P$. Namely, the matrix
$A=\left(a^{ij}\right)$ induces a Riemannian metric $g_{A}$ on the manifold
$M$ by

$\left(g(\partial/\partial x_{i},\partial/\partial
x_{j})\right):=\left(a_{ij}\right)=A^{-1}.$

We recall two fundamental properties of the heat kernel $k_{P}^{M}$ associated
with the operator $P$.

######  Lemma 2.2.

Let $P$ be a nonnegative second-order elliptic operator defined on $M$. Then
the minimal positive heat kernel $k_{P}^{M}(x,y,t)$ satisfies the following
properties:

  1. (1)

$k_{P}^{M}(x,y,t)$ satisfies the Chapman-Kolmogorov equation (the semigroup
property):

$k_{P}^{M}(x,y,s+t)=\int_{M}k_{P}^{M}(x,z,s)k_{P}^{M}(z,y,t)\,\mathrm{d}z\qquad\forall
s,t>0\mbox{ and }\forall\ x,y\in M.$

  2. (2)

$k_{P}^{M}(x,y,t)\geq 0\qquad\forall t>0\mbox{ and }\ \forall x,y\in M.$

Now we are in a position to state the main theorem of the paper.

######  Theorem 2.3.

Let $(M,g)$ be a smooth connected and noncompact manifold of dimension $n\geq
2.$ Let $P$ be a subcritical operator in $M$ of the form (2.1). Assume that
$G(x):=G_{P}^{M}(x,o)$, the minimal positive Green function with a singularity
at $o\in M$, satisfies

$\underset{x\rightarrow\bar{\infty}}{\lim}G(x)\,=\,0,$ (2.3)

where $\bar{\infty}$ is the ideal point of the one-point compactification of
$M$.

Then there exists a positive smooth function $\rho$ defined on $M$ such that
the operator $P_{\rho}:=\rho P$ satisfies the following properties:

$(1)$ $M$ is stochastically incomplete with respect to $P_{\rho}.$

$(2)$ $M$ is $L^{1}$-Liouville with respect to $P_{\rho}$ if and only if $M$
is $L^{1}$-Liouville with respect to $P.$

######  Remark 2.4.

Consider a 2-dimensional Riemannian manifold $(M,g)$, and consider a conformal
change of the metric $g\mapsto\hat{g}:=(\rho(x))^{-1}g$, where $\rho$ is a
positive smooth function on $M$. Then the Laplace-Beltrami operator on
$(M,\hat{g})$ is given by $\Delta_{\hat{g}}=\rho(x)\Delta_{g}$.

##  3\. Key Lemmas

In this section we prove some of the main ingredients used in the proof of
Theorem 2.3. Recently there have been extensive research on the equivalent
conditions of stochastic completeness/incompleteness (see [1, 14, 15, 16]). In
these papers it has been shown that stochastic completeness is equivalent to
Omori-Yau maximum principle. However, all these results are proved in the
particular case of $P:=-\Delta_{g}$, where $\Delta_{g}$ is the Laplace-
Beltrami operator in $M$. It turns out that one can easily extend with some
modifications these results to the case of second-order elliptic operators
satisfying appropriate assumptions. The next lemma is devoted to address this
characterization.

######  Lemma 3.1.

Let $M$ be a smooth, connected, noncompact manifold, and let $P$ be a
subcritical operator (not necessarily symmetric) of the form (2.1). Then the
following assertions are equivalent.

  1. (1)

$M$ is stochastically complete with respect to $P$.

  2. (2)

For every $\lambda<0,$ the only nonnegative, bounded classical solution $u$ of
$Pu=\lambda u$ is $u=0.$

  3. (3)

If for any $T\in(0,\infty),$ the Cauchy problem

$\displaystyle\begin{cases}&\frac{\partial v}{\partial t}+Pv=0,\\\
&v|_{t=0^{+}}=0,\end{cases}$

has a bounded solution in $M\times(0,T)$, then necessarily $v\,=\,0.$

  4. (4)

For $u\in C^{2}(M)$ with $\sup_{M}u<\infty$ and for $\alpha>0$, define

$\Omega_{\alpha}:=\\{x\in M\mid u(x)>\sup_{M}u-\alpha\\}.$ (3.1)

Then for every such $u$ and $\alpha>0$, we have
$\inf_{\Omega_{\alpha}}(-P)u\leq 0.$

  5. (5)

For every $u\in C^{2}(M)$ with $\sup_{M}u<\infty,$ there exists a sequence
$\\{x_{n}\\}$ such that for every $n\in\mathbb{N}$,

$u(x_{n})\geq\sup_{M}u-\frac{1}{n}\,,\quad\emph{and}\quad(-P)u(x_{n})\leq\frac{1}{n}\,.$

The above lemma is known for the case when $P:=-\Delta_{g}.$ A complete proof
can be found, for instance in ([5, Theorem 6.2], [16, Theorem 1.1]), see also
[6]. For the sake of completeness we present the proof below, with the
necessary modifications needed for a general operator $P$ of the form (2.1).

###### Proof.

$(2)\implies(1)$. Suppose that $(1)$ is not true. Then $M$ is stochastically
incomplete with respect to $P$. Define

$u(x,t)=\mathcal{P}_{t}1(x,t):=\int_{M}k^{M}_{P}(x,y,t)\,\mathrm{d}y.$

By our assumption, $0\leq u(x,t)<1$. For any $\lambda>0$ set,

$v(x):=\int_{0}^{\infty}e^{-\lambda
t}u(x,t)\mbox{d}t=\int_{M}G_{P+\lambda}^{M}(x,y)\,\mathrm{d}y.$

It can be easily checked (see [18]) that

$\displaystyle(P+\lambda)v(x)=1.$

Clearly, $0<v<1/\lambda$. Let $w:=1-\lambda v$. Then $(P+\lambda)w=0$ and
$0<w<1$, which contradicts $(2)$.

$(3)\implies(2)$. Suppose $(2)$ is not true, and let $v$ be a bounded nonzero
nonnegative function satisfying $Pv=\lambda v$ with $\lambda<0$. We may assume
that $0<v\leq 1$. Then $u(x,t):=e^{-\lambda t}v(x)$ solves the Cauchy problem

$\begin{cases}&\frac{\partial u}{\partial t}+Pu=0,\\\
&u|_{t=0^{+}}=v(x).\end{cases}$ (3.2)

Also, $w(x,t):=\int_{M}k^{M}_{P}(x,y,t)v(y)\,\mathrm{d}y$ solves the above
Cauchy problem. Using the minimality of $w$, we conclude that $0\leq w\leq 1$.
Moreover, we claim $u\neq w$ for $t>0$. Indeed, since
$\int_{M}k^{M}_{P}(x,y,t)\,\mathrm{d}y\leq 1$, it follows that

$\sup_{x\in M}w(x,t)\leq\sup_{x\in M}v(x).$

On the other hand, for $t>0$

$\sup_{x\in M}u(x,t)=e^{-\lambda t}\sup_{x\in M}v(x)>\sup_{x\in
M}v(x)\geq\sup_{x\in M}w(x,t),$

and the claim is proved. Hence, the nonzero nonnegative and bounded function
$z:=u-w$ solves the Cauchy problem

$\displaystyle\begin{cases}&\frac{\partial z}{\partial t}+Pz=0,\\\
&z|_{t=0^{+}}=0,\end{cases}$ (3.3)

which contradicts $(3).$

$(1)\implies(3).$ We also show this by contraposition. Assume $u(x,t)$ is a
nonzero bounded solution of (3.3) in $M\times(0,T)$ for some $T>0$. Without
loss of generality, we may assume that $\sup_{x\in M}u(x,t)>0$ and $\sup_{x\in
M}|u(x,t)|<1$. Then the function $w:=1-u$ is positive and $\inf_{x\in
M}w(x,t)<1$. Since the function $w$ is a solution to the Cauchy problem

$\displaystyle\begin{cases}&\frac{\partial w}{\partial t}+Pw=0,\\\
&w|_{t=0^{+}}=1,\end{cases}$ (3.4)

and $\mathcal{P}_{t}1:=\int_{M}K_{P}^{M}(x,y,t)\,\mathrm{d}y$ is the minimal
positive solution to (3.4), we conclude that $\mathcal{P}_{t}1\leq w$.
Therefore, for some $x\in M$ and $t\in(0,T)$,

$\mathcal{P}_{t}1=\int_{M}K_{P}^{M}(x,y,t)\,\mathrm{d}y<1,$

and M is stochastically incomplete.

$(4)\implies(5)$ and $(5)\implies(2)$ are trivial. Now to complete the proof
it remains to show $(2)\implies(4).$ We argue by a contradiction. Assume that
there exists a function $u\in C^{2}(M)$ with $\sup_{M}u<\infty$ satisfying for
some $\alpha>0$

$\inf_{\Omega_{\alpha}}(-Pu)\geq C>0,$

where $\Omega_{\alpha}$ is defined in (3.1). Then following [16, Theorem 1.1],
we set

$\Omega^{*}:=\left\\{x\in M:(-P)u>\frac{C}{2}\right\\}.$

Obviously, $\overline{\Omega_{\alpha}}\subset\Omega^{*}$. Moreover,
$u+\alpha-\sup_{M}u$ is a subsolution of $(P+\lambda)$ on $\Omega^{*}$ with
$\lambda=C/(2\alpha)$. Hence
$u_{\alpha}:=\mbox{max}\\{u+\alpha-\sup_{M}u,0\\}$ is a subsolution of
$(P+\lambda)$ in $M$. Clearly, $0<u_{\alpha}\leq\alpha.$ Furthermore, any
positive constant is a supersolution of $P+\lambda$ in $M$, and choosing a
constant strictly greater than $\alpha$ we have a supersolution
$u^{+}>u_{\alpha}.$

Using the “sub/supersolution method” and Perron’s method for $P+\lambda$, we
conclude that there exists a solution $v$ of the equation $(P+\lambda)u=0$ in
$M$ that satisfies $0<u_{\alpha}\leq v\leq u^{+}$ which does not vanish
identically, but this contradicts $(2)$. ∎

######  Remark 3.2.

In view of Lemma 3.1, stochastically incomplete manifold with respect to $P$
is characterized by the existence of a function $u\in C^{2}(M)$ with
$u^{*}:=\sup_{M}u<\infty$ such that for any sequence $\\{x_{n}\\}\in M$
satisfying

$\lim_{n\rightarrow\infty}u(x_{n})=u^{*}$ (3.5)

we have

$\limsup_{n\rightarrow\infty}(-P)u(x_{n})>0.$

Next, we address the $L^{1}$-Liouville property. Grigor’yan proved in [6, 5]
that the $L^{1}$-Liouville property with respect to the Laplace-Beltrami
operator $\Delta_{g}$ is equivalent to the non-integrability of the Green
function of $\Delta_{g}.$ Using Grigor’yan’s approach (see also [2]), we
extend in the following lemma the aforementioned result to the case of a
subcritical operator of the form (2.1). In fact, Corollary 3.4, demonstrates
that the $L^{1}$-Liouville property is equivalent to the non-integrability of
the heat kernel of $P$.

######  Lemma 3.3.

Let $(M,g)$ be a smooth, connected, noncompact manifold, and let $P$ be a
subcritical operator of the form (2.1). Then the following assertions are
equivalent.

  1. (1)

$M$ is not $L^{1}$-Liouville with respect to $P$.

  2. (2)

For any $y\in M$, the minimal Green function $G^{M}_{P}(\cdot,y)$ is in
$L^{1}(M)$.

###### Proof.

$(2)\implies(1).$ Define

$u:=\min\\{G^{M}_{P}(\cdot,y),1\\}.$

Since $P1=0,$ it follows that $u$ is a nonconstant positive
$L^{1}(M)$-supersolution. Hence, $M$ is not $L^{1}$-Liouville with respect to
$P$.

$(1)\implies(2).$ Since $P1=0$, the operator $P$ satisfies the weak and strong
maximum principle. Suppose that $u\in L^{1}(M)$ is a nontrivial nonnegative
supersolution of $P$ in $M$. Then by the strong maximum principle, $u>0$ in
$M$.

Fix $y\in M$ and a compact smooth exhaustion $\\{M_{j}\\}_{j=1}^{\infty}$ of
$M$ such that $K=\overline{B(y,1)}\Subset M_{1}$. Using Harnack’s inequality,

$G^{M}_{P}(x,y)\leq Cu(x),\quad\forall x\in\partial K,$

where $C:=C(K)$ is a positive constant. Then using the _generalized maximum
principle_ , we obtain $G^{M_{j}}_{P}(x,y)\leq Cu(x),\quad\forall x\in
M_{j}\setminus K$, where $G^{M_{j}}_{P}(y,x)$ denotes the Dirichlet Green
function in $M_{j}.$ Now letting $j\rightarrow\infty$, we obtain

$G^{M}_{P}(x,y)\leq Cu(x),\quad\forall x\in M\setminus K.$

Therefore,

$\displaystyle\int_{M}G^{M}_{P}(x,y)\,\mathrm{d}x$
$\displaystyle=\int_{K}G^{M}_{P}(x,y)\,\mathrm{d}x+\int_{M\setminus
K}G^{M}_{P}(x,y)\,\mathrm{d}x$
$\displaystyle\leq\int_{K}G^{M}_{P}(x,y)\,\mathrm{d}x+C\int_{M\setminus
K}u(x)\,\mathrm{d}x<\infty$

since $G^{M}_{P}$ is locally integrable. Hence, $G^{M}_{P}(\cdot,y)$ is in
$L^{1}(M)$. ∎

In view of the above lemma, we have.

######  Corollary 3.4.

Let $M$ be a smooth, connected, noncompact manifold, and let $P$ be an
elliptic operator of the form (2.1). Consider its heat kernel
$k_{P}^{M}(x,y_{,}t)$. Then $M$ is $L^{1}$-Liouville with respect to $P$ if
and only if there exists $y_{0}\in M$ such that

$\int_{M}\int_{0}^{\infty}k_{P}^{M}(x,y_{0},t)\,\mathrm{d}t\,\mathrm{d}x=\infty.$

######  Corollary 3.5.

Let $P$ be a symmetric operator in $M$ of the form (2.1). If $P$ is
stochastically complete, then $M$ is $L^{1}$-Liouville with respect to $P$.

###### Proof.

First consider the case when $P$ is subcritical, fix $y\in M$. By Lemma 3.3,
it suffices to show that $G(\cdot,y)$ is not integrable. Using Tonelli’s
theorem we get

$\displaystyle\int_{M}G_{P}^{M}(x,y)\,\mathrm{d}x=\\!\int_{M}\int_{0}^{\infty}k_{P}^{M}(x,y,t)\,\mathrm{d}t\,\mathrm{d}x=\\!\int_{0}^{\infty}\int_{M}k_{P}^{M}(x,y,t)\,\mathrm{d}x\,\mathrm{d}t=\\!\int_{0}^{\infty}\\!\\!\,\mathrm{d}t=\infty.$

On the other hand, by Remark 1.4, $P$ is critical in $M$ if and only if $P$
admits (up to a multiplicative constant) a unique positive supersolution.
Therefore, since $P1=0$, it follows that in the critical case, $M$ is
$L^{1}$-Liouville with respect to $P$. ∎

We recall an elementary lemma from [17, Lemma 2.2].

######  Lemma 3.6.

Let $P$ be an operator of the form (2.1) and $\rho$ a positive smooth
function. Assume that $P$ is subcritical in $M$. Consider the operator
$P_{\rho}=\rho P$ defined on $M$. Then $P_{\rho}$ is subcritical in $M$ and
the minimal positive Green function $G^{M}_{P_{\rho}}$ of $P_{\rho}$ in $M$ is
given by

$G^{M}_{P_{\rho}}(x,y)=\frac{G^{M}_{P}(x,y)}{\rho(y)}\,.$ (3.6)

##  4\. Proof of the main theorem

This section is devoted to the proof of Theorem 2.3. The proof is essentially
based on finding a “test function” $u\in C^{2}(M)$ which satisfies the
hypothesis of Remark 3.2. An alternative method involving optimal Hardy
weights is presented in Section 6.

###### Proof of Theorem 2.3.

(1) Fix $o\in M$ and a positive Radon measure $\mu$ on $M$ with a smooth
positive density $\mu$. We assume that the corresponding Green potential
$G_{\mu}$ is finite. That is, we assume that for some $x\in M$ (and therefore,
for any $x\in M$)

$G_{\mu}(x):=\int_{M}\mbox{$G^{M}_{P}(x,y)$}\mu(y)\,\mathrm{d}y<\infty.$ (4.1)

Hence, $PG_{\mu}=\mu$. Denote by $G(x)=\mbox{$G^{M}_{P}(x,o)$}$. By the
minimality of $G$, it follows that for $\varepsilon>0$ small enough there
exists $C_{\varepsilon}>0$ such that $G(x)\leq C_{\varepsilon}G_{\mu}(x)$ in
$M\setminus B(o,\varepsilon)$.

We further assume that $G_{\mu}(x)\leq C_{1}G(x)$ in $M$. In other words, we
have

$G_{\mu}\asymp G\qquad\mbox{in }M\setminus B(o,\varepsilon).$ (4.2)

In particular, (2.3) and (4.2) imply that $G_{\mu}$ is bounded in $M$ and

$\lim_{x\rightarrow\bar{\infty}}G_{\mu}(x)=0.$

For necessary and sufficient conditions for the validity of (4.2) see [4,
Lemma 6.1] and references therein.

Set $\rho:=1/\mu$, and consider the operator $P_{\rho}=\rho P$. Lemma 3.6
implies that $P_{\rho}$ is subcritical in $M$ with the minimal positive Green
function

$\mbox{$G^{M}_{P_{\rho}}(x,y)$}=\frac{\mbox{$G^{M}_{P}(x,y)$}}{\rho(y)}\,.$

In view of Remark 3.2, in order to prove stochastic incompleteness of
$P_{\rho}$ in $M$, it is sufficient to construct a function $u\in C^{2}(M)$
with $u^{*}:=\sup_{M}u<\infty$ such that for any sequence $\\{x_{n}\\}\in M$
satisfying $\lim_{n\rightarrow\infty}u(x_{n})=u^{*}$, we have

$\limsup_{n\rightarrow\infty}(-P_{\rho})u(x_{n})>0.$

Let $u:=-G_{\mu}$. Then clearly $u<0$ and $u$ satisfies $u(x)\rightarrow 0$ as
$x\rightarrow\bar{\infty}$, where $\bar{\infty}$ is the ideal point of the
one-point compactification of $M$. Hence, $u^{*}:=\sup_{M}u=0<\infty$ and a
sequence $\\{x_{n}\\}\subset M$ satisfies
$\lim_{n\to\infty}u^{*}(x_{n})=u^{*}=0$ if and only if $x_{n}\to\bar{\infty}$.
On the other hand, since $(-P_{\rho})u=1$, it follows that for any sequence
$\\{x_{n}\\}$ we have

$\limsup_{n\rightarrow\infty}(-P_{\rho})u(x_{n})=1>0.$

Therefore, $M$ is stochastic incomplete with respect to $P_{\rho}$ .

(2) We need to prove that $M$ is $L^{1}$-Liouville with respect to the
operator $P_{\rho}:=\rho P$ if and only if $P$ is $L^{1}$-Liouville with
respect to $P$. In view of Lemma 3.3 this is equivalent to the non-
integrability of the Green function $G^{M}_{P_{\rho}}(x,y)$ as a function of
$x$. Indeed,

$\hskip 56.9055pt\int_{M}\mbox{$G^{M}_{P_{\rho}}(x,y)$}\,\mathrm{d}x=\int_{M}\
\frac{\mbox{$G^{M}_{P}(x,y)$}}{\rho(y)}\,\mathrm{d}x=\mu(y)\int_{M}G^{M}_{P}(x,y)\,{\rm
d}x.$

For fixed $y,$ the conclusion follows immediately. ∎

######  Remark 4.1.

In light of the above proof, it follows that Theorem 2.3 holds true for any
$\rho=1/\mu$ such that $\lim_{x\rightarrow\bar{\infty}}G_{\mu}(x)=0$ (without
assuming (4.2)).

######  Remark 4.2.

Let $W_{\mu}:=\mu/G_{\mu}$. By [4, Lemma 6.2], if (4.2) holds true, then the
operator $P-W_{\mu}$ is positive-critical in $M$ with respect to the Hardy-
weight $W_{\mu}$ and $G_{\mu}$ is its ground state.

##  5\. Skew product operators and stochastically incompleteness

This section is devoted to the study of the skew product of operators $P_{i}$
of the form (2.1), defined on smooth, noncompact, and connected manifolds
$M_{i}$, $i=1,2$, and its relations to the stochastic completeness and
$L^{1}$-Liouville properties of the operators $P_{i}$. Let us recall the
definition of skew product (see [10, 11, 12]).

######  Definition 5.1.

Let $M=M_{1}\times M_{2}$ be the Cartesian product of two manifolds $M_{1}$
and $M_{2}$, and denote a point $x$ in $M$ by $x=(x_{1},x_{2})\in
M=M_{1}\times M_{2}.$ The skew product operator $P$ of operators $P_{i}$ of
the form (2.1) defined on $M_{i}$, $i=1,2$, is given by

$P:=P_{1}\otimes I_{2}+I_{1}\otimes P_{2},$ (5.1)

where $I_{i}$ is the identity map on $M_{i}$.

Next, we recall a well known lemma on the heat kernel of a skew product
operator. For completeness, we provide the proof.

######  Lemma 5.2.

Let $P_{1}$ and $P_{2}$ be nonnegative operators of the form (2.1) defined on
$M_{1}$ and $M_{2}$, respectively. Let $k_{P_{1}}^{M_{1}}$ and
$k_{P_{2}}^{M_{2}}$ be the minimal positive heat kernels of $P_{1}$ and
$P_{2}$, respectively. Then the minimal positive heat kernel $k^{M}_{P}$ of
the skew product operator $P$ on $M=M_{1}\times M_{2}$ satisfies

$k_{P}^{M}(x,y,t)=k_{P_{1}}^{M_{1}}(x_{1},y_{1},t)\>k_{P_{2}}^{M_{2}}(x_{2},y_{2},t),$
(5.2)

where $x=(x_{1},x_{2}),\>y=(y_{1},y_{2})\in M.$ Furthermore, if at least one
of the above operators is subcritical, then $P$ is subcritical in $M$.

###### Proof.

The product formula (5.2) follows immediately from the definition of skew
product. In order to prove subcriticality, we need to prove that for some
$x\neq y$, (and therefore for any $x\neq y$), $x,y\in M$,

$\int_{0}^{\infty}k_{P}^{M}(x,y,t)\,\mathrm{d}t=\int_{0}^{\infty}k_{P_{1}}^{M_{1}}(x_{1},y_{1},t)\>k_{P_{2}}^{M_{2}}(x_{2},y_{2},t)\,\mathrm{d}t<\infty.$

Without loss of generality, we may assume that $P_{2}$ is subcritical in
$M_{2}$. By Fubini’s theorem it is enough to prove that

$\int_{M_{1}}\int_{0}^{\infty}k_{P_{1}}^{M_{1}}(x_{1},y_{1},t)\>k_{P_{2}}^{M_{2}}(x_{2},y_{2},t)\,\mathrm{d}t\>\,\mathrm{d}y_{1}<\infty.$

By the minimality of the heat kernel, it follows that

$\int_{M_{1}}k_{P_{1}}^{M_{1}}(x_{1},y_{1},t)\,\mathrm{d}y_{1}\leq 1\>\text{
for some (and hence for all) }\ (x_{1},t)\in M_{1}\times(0,\infty).$

Therefore, by Tonelli’s theorem and the subcriticality of $P_{2}$ in $M_{2}$,
we obtain,

$\displaystyle\int_{M_{1}}\int_{0}^{\infty}k_{P_{1}}^{M_{1}}(x_{1},y_{1},t)\>k_{P_{2}}^{M_{2}}(x_{2},y_{2},t)\,\mathrm{d}t\>\,\mathrm{d}y_{1}\leq
1\cdot\int_{0}^{\infty}k_{P_{2}}^{M_{2}}(x_{2},y_{2},t)\,\mathrm{d}t<\infty.\qquad\qed$

The next result leads us to the question of stochastic (in)completeness of a
skew product operator.

######  Lemma 5.3.

Let $P_{1}$ and $P_{2}$ be operators of the form (2.1) defined on $M_{1}$ and
$M_{2}$, respectively. Then the manifold $M=M_{1}\times M_{2}$ is
stochastically incomplete with respect to the skew product operator $P$ if at
least one of the manifolds $M_{i}$ is stochastically incomplete with respect
to $P_{i}$.

###### Proof.

Without loss of generality, we may assume that $M_{1}$ is stochastically
incomplete manifold with respect to the operator $P_{1}$, and hence for all
$(x_{1},t)\in M_{1}\times(0,\infty)$ there holds

$\int_{M_{1}}k_{P_{1}}^{M_{1}}(x_{1},y_{1},t)\,\mathrm{d}y_{1}<1.$ (5.3)

Now for $x=(x_{1},x_{2}),\>y=(y_{1},y_{2})\in M$, by Tonelli’s theorem we have

$\displaystyle\int_{M}k_{P}^{M}(x,y,t)\,\mathrm{d}y$
$\displaystyle=\left(\int_{M_{1}}k_{P_{1}}^{M_{1}}(x_{1},y_{1},t)\,\mathrm{d}y_{1}\right)\\!\\!\left(\int_{M_{2}}\>k_{P_{2}}^{M_{2}}(x_{2},y_{2},t)\,\mathrm{d}y_{2}\right)$
$\displaystyle\leq\int_{M_{1}}k_{P_{1}}^{M_{1}}(x_{1},y_{1},t)\,\mathrm{d}y_{1}<1.$

Another simple alternative proof can be easily derived using Remark 3.2. ∎

The following result concerns the interplay between stochastic
(in)completeness and the $L^{1}$-Liouville property of a skew product operator
of the form (5.1).

######  Theorem 5.4.

Let $P_{1}$ and $P_{2}$ be operators of the form (2.1) defined on $M_{1}$ and
$M_{2}$, respectively. Let $P:=P_{1}\otimes I_{2}+I_{1}\otimes P_{2}$ defined
on $M:=M_{1}\times M_{2}$ be the corresponding skew product operator. Then the
following assertions hold true:

  1. (1)

If at least one of $M_{i},$ $i=1,2$ is not $L^{1}$-Liouville with respect to
$P_{i}$, and the other operator is symmetric, then the skew product operator
$P$ on $M$ is not $L^{1}$-Liouville.

  2. (2)

If one of the $M_{i},$ $i=1,2$ is $L^{1}$-Liouville with respect to $P_{i}$
and the other is symmetric, then the product manifold $M$ is $L^{1}$-Liouville
with respect to $P$.

  3. (3)

Assume that both $M_{i}$, $i=1,2$, are stochastically complete with respect to
$P_{i}$, then the product manifold $M$ is stochastically complete with respect
to $P$. Moreover, if $P$ is symmetric in $M$, then $M$ is $L^{1}$-Liouville
with respect to $P$.

  4. (4)

If one of the $M_{i},$ $i=1,2$ is stochastically incomplete with respect to
$P_{i}$, then $M$ is stochastically incomplete with respect to $P$.

###### Proof.

$(1)$ Without loss of generality, we assume that $M_{1}$ is not
$L^{1}$-Liouville with respect to $P_{1}$ and $P_{2}$ is symmetric. Then

$\displaystyle\int_{M}G_{P}^{M}(x_{1},x_{2},y_{1},y_{2})\,\mathrm{d}x_{1}\,\mathrm{d}x_{2}=\int_{0}^{\infty}\int_{M_{1}}\int_{M_{2}}k_{P_{1}}^{M_{1}}(x_{1},y_{1},t)\>k_{P_{2}}^{M_{2}}(x_{2},y_{2},t)\,\mathrm{d}x_{2}\,\mathrm{d}x_{1}\,\mathrm{d}t$
$\displaystyle=\int_{0}^{\infty}\int_{M_{1}}k_{P_{1}}^{M_{1}}(x_{1},y_{1},t)\,\mathrm{d}x_{1}\left(\int_{M_{2}}k_{P_{2}}^{M_{2}}(x_{2},y_{2},t)\,\mathrm{d}x_{2}\right)\,\mathrm{d}t$
$\displaystyle\leq\int_{M_{1}}\int_{0}^{\infty}k_{P_{1}}^{M_{1}}(x_{1},y_{1},t)\,\mathrm{d}t\,\mathrm{d}x_{1}<\infty.$

Hence, $M$ is not $L^{1}$-Liouville with respect to $P$.

$(2)$ Without loss of generality, we assume that $M_{1}$ is $L^{1}$-Liouville
with respect to $P_{1}$ and $P_{2}$ is symmetric. Let $(y_{1},y_{2})$ be any
point on $M_{1}\times M_{2}$ and assume $P$ is not $L^{1}$-Liouville, i.e.,

$\int_{M_{1}}\int_{M_{2}}\int_{0}^{\infty}k_{P_{1}}^{M_{1}}(x_{1},y_{1},t)\>k_{P_{2}}^{M_{2}}(x_{2},y_{2},t)\,\mathrm{d}t\,\mathrm{d}x_{1}\,\mathrm{d}x_{2}<\infty.$

Then using Fubini’s theorem and the symmetricity of $P_{2}$, we conclude,

$\int_{M_{1}}\int_{0}^{\infty}k_{P_{1}}^{M_{1}}(x_{1},y_{1},t)\,\mathrm{d}t\,\mathrm{d}x_{1}<\infty,$

and hence $M_{1}$ is not $L^{1}$-Liouville with respect to $P_{1}$ which is a
contradiction. Therefore, $M$ is $L^{1}$-Liouville with respect to $P$.

$(3)$ Since both $M_{i}$, $i=1,2$ are stochastically complete with respect to
$P_{i}$, therefore, we have for $i=1,2,$

$\int_{M_{i}}k_{P_{i}}^{M_{i}}(x_{i},y_{i},t)\,\mathrm{d}y_{i}=1\text{ for any
}(x_{i},t)\in M_{i}\times(0,\infty).$

Using Tonelli’s theorem we obtain

$\displaystyle\int_{M}k_{P}^{M}(x_{1},x_{2},y_{1},y_{2},t)\,\mathrm{d}y_{1}\,\mathrm{d}y_{2}=\int_{M_{1}}\int_{M_{2}}k_{P_{1}}^{M_{1}}(x_{1},y_{1},t)\>k_{P_{2}}^{M_{2}}(x_{2},y_{2},t)\,\mathrm{d}y_{2}\,\mathrm{d}y_{1}=1.$

So, $M$ is stochastically complete with respect to $P$. Furthermore, under the
symmetricity assumption, Corollary 3.5 implies that $P$ is $L^{1}$-Liouville
as well.

$(4)$ By Lemma 5.3, $M$ is stochastically incomplete with respect to $P$. ∎

The following result concerns a case where the $L^{1}$-Liouville property
holds for the skew product operator $P$.

######  Corollary 5.5.

If one of the $P_{i}$, $i=1,2$ is symmetric in $M_{i}$ and the other is
critical, then the product manifold $M$ is $L^{1}$-Liouville with respect to
$P$.

###### Proof.

Since a critical operator of the form (2.1) is $L^{1}$-Liouville, the
corollary follows from part $(2)$ of Theorem 5.4. ∎

######  Remark 5.6.

We recall that in [5, 17] the authors present a counterexample to a conjecture
of A. Grigor’yan of a skew product manifold $M\times K$ that does not satisfy
the $L^{\infty}$-Liouville property, where $(M,g)$ is a complete Riemannian
manifold (with $\lambda_{0}(-\Delta_{g},M)=0$) satisfying the
$L^{\infty}$-Liouville property in $(M,g)$, and $K$ is any compact Riemannian
manifold.

##  6\. Optimal Hardy weight and Stochastically incomplete

In the present section, we give an alternative proof of Theorem 2.3 under an
additional assumption.

###### Alternative proof of Theorem 2.3.

Step 1 : First we recall the main result of Frass, Devyver and Pinchover [3,
Theorem 4.12] (also see [19]): Let $P$ be a subcritical operator of the form
(2.1) in $M,$ and let $G^{M}_{P}$ be its minimal positive Green function.
Assume that $P$ and $G^{M}_{P}$ satisfy all the hypotheses of Theorem 2.3. Let
$0\leq\varphi\in C_{c}^{\infty}(M),\varphi\neq 0$, and consider the Green
potential given by

$G_{\varphi}^{M}(x)=\int_{M}G^{M}_{P}(x,y)\varphi(y)\,{\rm d}y.$

Then $G_{\varphi}^{M}$ satisfies

$P(\sqrt{G_{\varphi}^{M}})=W\sqrt{G_{\varphi}^{M}}\geq 0,$

where $W:=\dfrac{P(\sqrt{G_{\varphi}^{M}})}{\sqrt{G_{\varphi}^{M}}}$, and
$W=\dfrac{|\nabla(G^{M}_{\varphi})|^{2}_{A}}{4|G^{M}_{\varphi}|^{2}}\quad\mbox{in}\
M\setminus\mathrm{supp}\,\varphi.$ Moreover, $W$ is an optimal Hardy-weight
for $P$ in $M$ (That is, $P-W$ is null-critical in $M$, and in particular, for
any $\lambda>1$ the operator $P-\lambda W\not\geq 0$ outside any compact set
in $M$).

We assume that $W\\!>\\!0$ outside a compact set $K$ satisfying
$\mathrm{supp}\,\varphi\Subset K\subset M$.

Step 2 : Let $\rho$ (to be chosen later) be a positive smooth function on $M.$
Consider the operator $P_{\rho}=\rho P.$ In view of Remark 3.2, in order to
prove stochastic incompleteness, it is sufficient to construct a function
satisfying certain conditions. Let $u$ be a strictly negative $C^{2}$-smooth
function such that

$u(x):=-(G_{\varphi}^{M}(x))^{1/2}\qquad\forall\ x\in M\setminus K.$

Then

$(-P_{\rho})u=\rho W(G^{M}_{\varphi})^{1/2}\qquad\forall\ x\in M\setminus K.$

Moreover, since $G_{\varphi}^{M}\asymp G^{M}_{P}$ in a neighbourhood of
$\bar{\infty}$, it follows from our assumption (2.3) that $u(x)\rightarrow 0$
as $x\rightarrow\bar{\infty}.$ Hence, $u^{*}:=\sup_{M}u=0$, and
$\lim_{n\to\infty}u(x_{n})=u^{*}$ if and only if
$x_{n}\rightarrow\bar{\infty}$. For any such a sequence, we eventually have,

$(-P_{\rho})u(x_{n})=\rho(x_{n})W(x_{n})(G^{M}_{\varphi}(x_{n}))^{1/2}.$

Choose a strictly positive smooth function $\rho$ on $M$ such that

$\rho(x):=\left(W(x)(G^{M}_{\varphi}(x))^{1/2}\right)^{-1}\qquad\forall x\in
M\setminus K,$

letting $n\rightarrow\infty,$ we conclude

$\limsup_{n\rightarrow\infty}(-P_{\rho})u(x_{n})=1>0.$

Hence, by Remark 3.2, $M$ is stochastically incomplete with respect to
$P_{\rho}:=\rho P$. ∎

######  Remark 6.1.

We note that the extra assumption that $W>0$ in a neighbourhood of
$\bar{\infty}$ is satisfied in many cases. For example, if $M$ is a smooth
bounded domain in $\mathbb{R}^{n}$, and the coefficients of $P$ are up to the
boundary Hölder continuous, then by the Hopf lemma, $W>0$ in a neighbourhood
of $\partial M$.

Acknowledgments. D. Ganguly is partially supported by INSPIRE faculty
fellowship (IFA17-MA98). Y.P. acknowledges the support of the Israel Science
Foundation (grant 637/19) founded by the Israel Academy of Sciences and
Humanities. P. Roychowdhury is supported by Council of Scientific & Industrial
Research (File no. 09/936(0182)/2017-EMR-I).

## References

  * [1] C. Bär, and G. P. Bessa, _Stochastic completeness and volume growth_ , Proc. Amer. Math. Soc. 138 (2010), 2629–2640.
  * [2] C. Constantinescu, and A. Cornea, _“Potential Theory on Harmonic Spaces”_ , with a preface by H. Bauer. Die Grundlehren der mathematischen Wissenschaften, Band 158. Springer-Verlag, New York-Heidelberg, 1972.
  * [3] B. Devyver, M. Fraas, and Y. Pinchover, _Optimal Hardy weight for second-order elliptic operator: an answer to a problem of Agmon_ , J. Funct. Anal. 266 (2014), 4422–4489.
  * [4] D. Ganguly, and Y. Pinchover, _Some new aspects of perturbation theory of positive solutions of second-order linear elliptic equations_ , in a special issue dedicated to the memory of A. I. Volpert, Pure Appl. Funct. Anal. 5 (2020), 295–319.
  * [5] A. Grigory’an, _Analytic and geometric background of recurrence and non-explosion of the Brownian motion on Riemannian manifolds_ , Bull. Amer. Math. Soc. (N.S.) 36 (1999), 135–249.
  * [6] A. Grigory’an, X. Huang, and J. Masamune, _On stochastic completeness of jump processes_ , Math. Z. 271 (2012), 1211–1239.
  * [7] G. Grillo, K. Ishige, and M. Muratori, _Nonlinear characterizations of stochastic completeness_ , J. Math. Pures Appl. (9) 139 (2020), 63–82.
  * [8] D. Impera, S. Pigola, and A. G. Setti, _Potential theory for manifolds with boundary and applications to controlled mean curvature graphs_ , J. Reine Angew. Math. 733 (2017), 121–159.
  * [9] L. Mari, and D. Valtorta, _On the equivalence of stochastic completeness and Liouville and Khas’minskii conditions in linear and nonlinear settings_ , Trans. Amer. Math. Soc. 365 (2013), 4699–4727.
  * [10] M. Murata, _On construction of Martin boundaries for second order elliptic equations_ , Publ. Res. Inst. Math. Sci. 26 (1990), 585–627.
  * [11] M. Murata, and T. Tsuchida, _Positive solutions of Schrödinger equations and Martin boundaries for skew product elliptic operators_ , Adv. Differential Equations 22 (2017), 621–692.
  * [12] M. Murata, and T. Tsuchida, _Monotonicity of non-Liouville property for positive solutions of skew product elliptic equations_ , Proc. Roy. Soc. Edinburgh Sect. A 150 (2020), 1429–1449.
  * [13] G. Nenciu, and I. Nenciu, _Drift-diffusion equations on domains in $\mathbb{R}^{d}$: essential self-adjointness and stochastic completeness_, J. Funct. Anal. 273 (2017), 2619–2654.
  * [14] G. Pacelli Bessa, S. Pigola, and A.G. Setti, _On the $L^{1}$-Liouville property of stochastically incomplete manifolds_, Potential Anal. 39 (2013), 313–324.
  * [15] L. F. Pessoa, S. Pigola, and A. G. Setti, _Dirichlet parabolicity and $L^{1}$-Liouville property under localized geometric conditions_, J. Funct. Anal. 273 (2017), 652–693.
  * [16] S. Pigola, M. Rigoli, and A. G. Setti, _A remark on the maximum principle and stochastic completeness_ , Proc. Amer. Math. Soc. 131 (2003), 1283–1288.
  * [17] Y. Pinchover, _On nonexistence of any $\lambda_{0}$-invariant positive harmonic function, a counterexample to Stroock’s conjecture_, Comm. Partial Differential Equations 20 (1995), 1831–1846.
  * [18] Y, Pinchover, _Maximum and anti-maximum principles and eigenfunctions estimates via perturbation theory of positive solutions of elliptic equations_ , Math. Ann. 314 (1999), 555–590.
  * [19] Y. Pinchover, and I. Versano, _On families of optimal Hardy-weights for linear second-order elliptic operators_ , J. Funct. Anal. 278 (2020), 108428, 27 pp..
  * [20] R. G. Pinsky, _“Positive Harmonic Functions and Diffusion”_ , Cambridge Studies in Advanced Mathematics 45, Cambridge University Press, Cambridge, 1995.


# Localization-delocalization effects of a delocalizing dissipation on
disordered XXZ spin chains

Xiansong Xu  Science, Mathematics and Technology Cluster, Singapore University
of Technology and Design, 8 Somapah Road, 487372 Singapore  Dario Poletti
<EMAIL_ADDRESS>Science, Mathematics and Technology Cluster,
Singapore University of Technology and Design, 8 Somapah Road, 487372
Singapore  EPD Pillar, Singapore University of Technology and Design, 8
Somapah Road, 487372 Singapore

###### Abstract

The interplay between interaction, disorder, and dissipation has shown a rich
phenomenology. Here we investigate a disordered XXZ spin chain in contact with
a bath which, alone, would drive the system towards a highly delocalized and
coherent Dicke state. We show that there exist regimes for which the natural
orbitals of the single-particle density matrix of the steady state are all
localized in the presence of strong disorders, either for weak interaction or
strong interaction. We show that the averaged steady-state occupation in the
eigenbasis of the open system Hamiltonian could follow an exponential decay
for intermediate disorder strength in the presence of weak interactions, while
it is more evenly spread for strong disorder or for stronger interactions.
Last, we show that strong dissipation increases the coherence of the steady
states, thus reducing the signatures of localization. We capture such
signatures of localization also with a concatenated inverse participation
ratio which simultaneously takes into account how localized are the
eigenstates of the Hamiltonian, and how close is the steady state to an
incoherent mixture of different energy eigenstates.

> As discovered in 1958 by P. W. Anderson Anderson (1958), in the presence of
> disorder, quantum systems can become insulators, or in other words the
> systems become localized. When the particles in the quantum system interact
> with each other, the effects of disorders can be very different from the
> non-interacting case, known as the many-body localization Basko, Aleiner,
> and Altshuler (2006), and the study of this topic is still attracting
> significant attention. In recent years, researchers have also studied how
> the environment surrounding such disordered and interacting systems can
> affect the localization properties of quantum systems. It was found that
> depending on the magnitude of disorder and on the environment, one can still
> observe possible signatures of localization in these quantum systems. Here
> we study the effects of a particular type of environment which tries to
> drive the system towards a highly delocalized state, while the presence of
> disorder, together with the interaction between the particles, could result
> in a localized system. We find that even when the system shows signatures of
> localization, there can be interesting differences in its long-time
> properties.

##  I Introduction

The effects of disorder in quantum systems took a prominent place in 1958 with
the discovery of Anderson localization Anderson (1958) in non-interacting
systems. The characterization of the effects of disorder in interacting
systems occurred much later Basko, Aleiner, and Altshuler (2006); Oganesyan
and Huse (2007); Pal and Huse (2010); Giamarchi and Schulz (1987, 1988) with
the study of many-body localized (MBL) systems. Some reviews on this phase of
matter which does not thermalize are Refs. Nandkishore and Huse, 2015; Abanin
and Papić, 2017; Abanin  _et al._ , 2019. MBL systems have also been studied
experimentally with ultracold atoms Schreiber  _et al._ (2015a, b); Bordia
_et al._ (2016); Choi  _et al._ (2016); Bordia  _et al._ (2017a, b); Lüschen
_et al._ (2017), trapped ions Smith  _et al._ (2016) and superconducting
qubits Roushan  _et al._ (2017); Guo  _et al._ (2020).

A natural question to ask is how MBL systems respond when coupled to a bath
Nandkishore, Gopalakrishnan, and Huse (2014); Johri, Nandkishore, and Bhatt
(2015); Nandkishore (2015); Žnidarič, Scardicchio, and Varma (2016);
Medvedyeva, Prosen, and Žnidarič (2016); Levi  _et al._ (2016); Fischer,
Maksymenko, and Altman (2016); Lüschen  _et al._ (2017); Hyatt  _et al._
(2017); Everest  _et al._ (2017); Luitz, Huveneers, and De Roeck (2017);
Nandkishore and Gopalakrishnan (2017); Žnidarič  _et al._ (2017); Žnidarič and
Ljubotina (2018); Karlsson, Hopjan, and Verdozzi (2018); Vakulchyk  _et al._
(2018); Xu, Guo, and Poletti (2018); Wu  _et al._ (2019); Rubio-Abadal  _et
al._ (2019); Wybo, Knap, and Pollmann (2020). Studies of strongly interacting
and strongly disordered systems coupled to small systems that act as baths
have shown that the overall system can become localized or thermalized
Nandkishore (2015); Hyatt  _et al._ (2017); Luitz, Huveneers, and De Roeck
(2017). Analyzing the spectral features of an MBL system coupled to a bath has
shown that there can exist regimes, e.g., with large enough interaction and
weak enough coupling to the baths, such that the spectral response to a local
spin flip is gapped, an indication of localization Nandkishore,
Gopalakrishnan, and Huse (2014); Johri, Nandkishore, and Bhatt (2015). This
can occur despite other quantities, such as the level statistics, do not show
signatures of localization for couplings to the baths which decrease
exponentially with the system size. Coupling MBL systems to local dissipative
boundary driving has shown that in the ergodic phase one can observe both
diffusive and subdiffusive transport Žnidarič, Scardicchio, and Varma (2016);
Žnidarič  _et al._ (2017); Žnidarič and Ljubotina (2018). Even when an
interacting and disordered system is coupled to a bath that drives it to the
infinite temperature state, the relaxation dynamics can be significantly
slowed down, even possibly assuming a stretched exponential form Medvedyeva,
Prosen, and Žnidarič (2016); Levi  _et al._ (2016); Fischer, Maksymenko, and
Altman (2016); Everest  _et al._ (2017). It is thus clear that strongly
disordered and interacting systems can show a large variety of phenomenology
depending on the system and bath parameters.

Recently in Refs. Vakulchyk  _et al._ , 2018; Xu, Guo, and Poletti, 2018, the
authors studied a disordered XXZ spin chain coupled to a bath which, alone,
would tend to drive the system towards a highly delocalized and coherent Dicke
state. Such dissipator has been first proposed in Diehl  _et al._ , 2008, in
which it was shown how it could be realized in ultracold atoms experiments,
and later on it was implemented with a trapped-ion universal simulator
Schindler  _et al._ (2013). In Refs. Vakulchyk  _et al._ , 2018; Xu, Guo, and
Poletti, 2018 it was shown that despite the presence of such delocalizing
dissipator, for large enough disorder and interactions, it is possible to find
signatures compatible with the localization of the steady state. Notably, the
level spacing of the density matrix follows a Poisson distribution, and the
natural orbitals of the single-particle density matrix can be all
exponentially localized in steady states.

In this work we go more in-depth in this investigation. First we show that for
a small disorder the density matrix can show a very different amount of off-
diagonal terms in the spin configuration basis depending on the magnitudes of
disorder and interaction. In particular, for the isotropic case of the XXZ
chain and weak disorder, the steady state approaches a Dicke state. For large
disorders and/or large interactions, the steady state density matrix is,
instead, mostly diagonal. For strong disorder, we find two regions of the
parameter space in which the natural orbitals of the single-particle density
matrix are localized, either with weak or with strong interaction. We also
analyze the steady-state density matrix in the basis of the eigenstates of the
system Hamiltonian and we find that for moderately large disorders the steady-
state occupation in eigenenergy basis follows an approximately exponentially
decay, indicating that the steady state is mainly composed of few (localized)
many-body eigenstates of the Hamiltonian. In the regime of large interactions
instead, the steady state has a more evenly spread overlap with the
eigenstates of the Hamiltonian, indicating a clear difference from the
previous regime. We also show that increasing the dissipation strength for
this particular bath increases the coherence in the steady-state density
matrix thus significantly weakening the signatures of localization. We
quantify this with a concatenated inverse participation ratio which signals
whether the eigenstates of the Hamiltonian are localized and, at the same
time, the steady state is not a coherent (hence delocalized) mixture of these
energy eigenstates.

##  II Model

We consider a dissipative spin chain whose evolution is described by a master
equation

$\displaystyle\frac{\partial\hat{\rho}}{\partial t}$
$\displaystyle=\mathcal{L}\left[\hat{\rho}\right]=-\frac{{\rm
i}}{\hbar}\left[\hat{H},\hat{\rho}\right]+\mathcal{D}\left[\hat{\rho}\right],$
(1)

where $\mathcal{L}$ is the system’s Lindbladian, $\mathcal{D}$ the portion of
the master equation which includes the effects of dissipation and $\hat{H}$
the Hamiltonian. We consider a prototypical MBL Hamiltonian given by an XXZ
spin-$1/2$ chain with local disorders

$\displaystyle\hat{H}=\sum_{l=1}^{L-1}\left[J\left(\hat{\sigma}^{x}_{l}\hat{\sigma}^{x}_{l+1}+\hat{\sigma}^{y}_{l}\hat{\sigma}^{y}_{l+1}\right)+\Delta\hat{\sigma}^{z}_{l}\hat{\sigma}^{z}_{l+1}\right]+\sum_{l=1}^{L}h_{l}\hat{\sigma}^{z}_{l},$
(2)

where the elements of $\hat{\sigma}_{l}^{\alpha}$ are given by the Pauli
matrices for $\alpha=x,\;y,\;z$, $J$ is the tunneling amplitude, $\Delta$ is
the interaction strength and the random field $h_{l}$ is uniformly distributed
in $\left[-W,W\right]$ with $W$ characterizing the disorder strength. In the
following we use open boundary conditions. Throughout this article We consider
a system size of $L=10$ and work in the zero total magnetization sector. Note
that this Hamiltonian exactly maps to an interacting spinless fermionic chain
after Jordan-Wigner transformation Jordan and Wigner (1928); Lieb, Schultz,
and Mattis (1961). In the absence of dissipation, this model has been shown to
have many-body localized or ergodic phases depending on the strengths of
disorder and interaction Pal and Huse (2010); Bera  _et al._ (2015, 2017);
Lezama  _et al._ (2017); Lin  _et al._ (2018); Luitz, Laflorencie, and Alet
(2015); Serbyn, Papić, and Abanin (2015), although for a better understanding
of the properties of this system in the thermodynamic limit, further studies
on the role of finite size effects may be needed Serbyn, Papić, and Abanin
(2015); Šuntajs  _et al._ (2020); Abanin  _et al._ (2021); Panda  _et al._
(2019); Kiefer-Emmanouilidis  _et al._ (2021). We consider a dissipator in
Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) form Gorini, Kossakowski, and
Sudarshan (1976); Lindblad (1976) given by

$\displaystyle\mathcal{D}\left[\hat{\rho}\right]$
$\displaystyle=\gamma\sum_{l=1}^{L-1}\left(\hat{V}_{l,l+1}\hat{\rho}\hat{V}_{l,l+1}^{\dagger}-\frac{1}{2}\left\\{\hat{V}_{l,l+1}^{\dagger}\hat{V}_{l,l+1},\hat{\rho}\right\\}\right),$
(3)

where $\gamma$ indicates the coupling strength, $\\{.,.\\}$ indicates the
anticommutator, and the jump operators $\hat{V}_{l,l+1}$ are

$\displaystyle\hat{V}_{l,l+1}=\left(\hat{\sigma}^{+}_{l}+\hat{\sigma}^{+}_{l+1}\right)\left(\hat{\sigma}^{-}_{l}-\hat{\sigma}^{-}_{l+1}\right),$
(4)

where $\hat{\sigma}^{+}_{l}=\left(\hat{\sigma}^{x}_{l}+{\rm
i}\hat{\sigma}^{y}_{l}\right)/2$ and
$\hat{\sigma}^{-}=\left(\hat{\sigma}_{l}^{x}-{\rm
i}\hat{\sigma}_{l}^{y}\right)/2$. Such a dissipator has been first proposed in
bosonic systems Diehl  _et al._ (2008, 2010); Tomadin, Diehl, and Zoller
(2011), and subsequently further studied with bosons and also with fermions
and spins Yi  _et al._ (2012); Bardyn  _et al._ (2013). While in Ref. Diehl
_et al._ , 2008 the authors proposed a realization for bosonic systems by
immersion of the system in a superfluid bath, an actual experimental
realization was performed using an ion-trap universal quantum simulator
Schindler  _et al._ (2013). It is important to stress that proving that a
disordered open many-body system, with a large number of small energy gaps in
the spectrum, can be faithfully described by a Markovian master equation in
Eq. (1) with dissipator given in Eq. (3) is not a trivial task, however a
realization via Trotterization of the evolution operator Lloyd (1996) on a
quantum simulator is in any case possible as long as one uses small enough
time steps to approximate the exact evolution with the Trotterized one. It is
thus possible to study also experimentally the model from Eqs. (1-3).

The interest of analyzing the effects of dissipator given by Eq. (3) in a
system with Hamiltonian (2) is that, for large enough disorder, the system has
eigenstates with a large and sudden variability of local magnetization, but
the dissipator we have chosen $\hat{V}_{l,l+1}$ favors a balanced
magnetization on neighboring sites. Hence, the model (1-3) can be an important
test-bed to analyze possible properties of disordered many-body open quantum
systems. We stress that this system has been studied for a single particle
problem in Refs. Yusipov  _et al._ , 2017; Vershinina  _et al._ , 2017 where
it was shown that for large enough disorder the system showed signatures of
localization. Here, as in Refs. Vakulchyk  _et al._ , 2018; Xu, Guo, and
Poletti, 2018, we add the non-trivial effect of interactions. To compute the
steady-state $\hat{\rho}_{s}$ we used either exact diagonalization or time
evolutions and we concentrated on the sector with zero total magnetization so
as to probe the strong effects of the interaction.

##  III Results

Figure 1: Single realizations for the absolute value of steady-state density
matrix in the spin configurations basis with increasing interaction strength
(a, b) $\Delta=0.5J$, (c, d) $\Delta=1J$, (e, f) $\Delta=7J$. Left panels (a,
c, d) have weak disorder $W=0.1J$ whereas right panels (b, d, f) have disorder
strength of $W=12J$. The dissipation strength $\hbar\gamma=0.1J$ and the
system size $L=10$.  Figure 2: Averaged absolute value of density matrix over
100 disorder realizations in the spin configurations basis with increasing
interaction strength (a, b) $\Delta=0.5J$, (c, d) $\Delta=1J$, (e, f)
$\Delta=7J$. Left panels (a, c, d) have weak disorder $W=0.1J$ whereas right
panels (b, d, f) have disorder strength of $W=12J$. The dissipation strength
$\hbar\gamma=0.1J$ and the system size $L=10$.

To begin, we would like to give the reader an overview of the possible steady
state properties of the system for different magnitudes of disorder and
interaction. In Fig. 1 we plot the absolute values of the elements of the
density matrix in the basis formed by all possible spin configurations. In
these plots we represent the density matrix for a randomly chosen disorder
distribution, and the results are typical. In the left panels (a), (c), and
(e) we consider a small disorder strength $W=0.5J$, while in the right columns
we have significantly larger disorder $W=12J$. The different rows represent
possible interaction values: $\Delta=0.5J$ for panels (a) and (b), $\Delta=J$
for panels (c) and (d), and $\Delta=7J$ for panels (e) and (f). For weak
interaction and weak disorder, Fig. 1(a), one observes large coherence,
however, the steady state becomes even more striking for weak disorder and
$\Delta=J$ as the state approximates the Dicke state $|\psi_{D}\rangle$
(symmetric superposition of all the possible spin configurations) which is the
exact steady state in the absence of disorder Xu, Guo, and Poletti (2018). For
larger interaction or larger disorder, Figs. 1(b, d, e, f), the off-diagonal
elements are significantly reduced and the steady state is much more diagonal.
However we can observe that for a given disorder realization the density
matrix appears to have different structures in presence of large disorder and
weak interaction, Fig. 1(b), or vice versa, Fig. 1(e), or even both large
interaction and disorder, Fig. 1(f). In Figs. 2 we perform a parallel of Fig.
1 where we first average the density matrix over 100 disorder realizations
before plotting the absolute value of its elements. From Figs. 2 we can deduce
that for weak disorders the density matrices of different realizations are
very similar to each other, as if the steady state is not strongly affected by
the presence of disorder. We also observe that, on average, the density
matrices are mostly diagonal for other cases.

One way used to characterize isolated disordered many-body systems is via the
single-particle density matrix
$\rho_{\mathrm{sp}}^{j,k}=\langle\psi_{n}|\hat{\sigma}^{+}_{j}\hat{\sigma}^{-}_{k}|\psi_{n}\rangle$
for each eigenstate $|\psi_{n}\rangle$ of the Hamiltonian Bera  _et al._
(2015, 2017); Lezama  _et al._ (2017); Lin  _et al._ (2018). In particular, it
was shown that for large enough disorder with either weak or strong
interaction, the natural orbitals of the single-particle density matrix, i.e.,
$\Psi_{\alpha}$ such that
$\rho_{\mathrm{sp}}\Psi_{\alpha}=n_{\alpha}\Psi_{\alpha}$ where $n_{\alpha}$
is the normalized occupation spectrum, all become exponentially localized. As
done in Ref. Xu, Guo, and Poletti, 2018, one can characterize in the same
manner the steady state of an open and disordered many-body quantum system
from the single-particle density matrix $\rho_{\mathrm{sp}}^{j,k}={\rm
tr}(\hat{\sigma}^{+}_{j}\hat{\sigma}^{-}_{k}\hat{\rho}_{\rm s})$. We then
compute a weighted inverse participation ratio

$\displaystyle\mathcal{I}(\hat{\rho}_{\rm
s})=\sum_{\alpha,l}n_{\alpha}|\Psi_{\alpha}(l)|^{4}$ (5)

and its average over the disorder realizations $\bar{\mathcal{I}}$.

Figure 3: Color map of averaged weighted inverse participation ratio
$\bar{\mathcal{I}}$ as a function of disorder strength $W$ and interaction
strength $\Delta$. Parameters are $L=10$, $\hbar\gamma=0.1J$, and $100$
disorder realizations.

In Fig. 3 we show a color map of the averaged IPR, $\bar{\mathcal{I}}$, versus
the interaction $\Delta$ and disorder strength $W$. For weak disorders we
observe that $\bar{\mathcal{I}}$ is small indicating the presence of
delocalized natural orbitals in the steady state. For intermediate values of
interaction, around $\Delta/J\approx 4$ we observe that such delocalized
orbitals are relevant in the steady state even for larger values of disorders.
For large disorder $W$ we observe that either weak or strong interaction
$\Delta$ results in large IPR.

Figure 4: Averaged diagonal elements of density matrix
$\bar{\rho}^{{}_{E}}_{i,i}$ over 100 disorder realizations versus the energy
level $i$ for increasing value of the energy. In (a) $\Delta=0.5J$, $W=7J$,
(b) $\Delta=0.5J$, $W=12J$, and (c) $\Delta=15J$, $W=7J$. The grey shaded
region demonstrates the standard deviation from the averaged results. The
values of $\bar{\mathcal{I}}$ for the cases shown in the different panels are
(a) $\bar{\mathcal{I}}=0.60\pm 0.01$, (b) $\bar{\mathcal{I}}=0.73\pm 0.01$ and
(c) $\bar{\mathcal{I}}=0.713\pm 0.008$.

At this point we try to characterize whether the steady states of the system
have similar properties in the regions with large $\bar{\mathcal{I}}$, whether
for small or large interaction $\Delta$ in the presence of strong disorder. To
do so we first define the steady state density matrix in the energy basis
$\rho^{{}_{E}}_{i,j}$ as

$\displaystyle\hat{\rho}_{s}=\sum_{i,j}\rho^{{}_{E}}_{i,j}|\psi_{i}\rangle\langle\psi_{j}|$
(6)

where $|\psi_{i}\rangle$ is the $i$-th eigenstates of the Hamiltonian. In Fig.
4 we then show the diagonal terms of $\rho^{{}_{E}}_{i,i}$ in ascending order
of the eigenenergies. The black curve represents the average over 100 disorder
realizations $\bar{\rho}^{{}_{E}}_{i,i}$ while the gray shading is the
standard deviation, which follows closely the functional form of the average
values. We observe that for $W=7J$ and $\Delta=0.5J$, Fig. 4(a), the steady
state has an exponentially decreasing overlap with the eigenstates of the
Hamiltonian. This indicates that few (localized) states dominate the steady
state. For the Hamiltonian chosen and the particular type of non-thermal
dissipator we are studying, the most occupied states tend to be those with
higher energy. For larger disorder $W=12J$ and $\Delta=0.5J$, Fig. 4(b), the
exponential decay is less pronounced. In Fig. 4(c), for $W=7J$ and
$\Delta=15J$ the slope of the profile has an overall small tilt due to
interactions. This indicates that the steady state is closer to an even
mixture of energy eigenstates. Hence, while both for weak and large
interaction $\bar{\mathcal{I}}$ is large, the steady state can have
significantly different properties.

Figure 5: (a) Averaged concatenated inverse participation ratio
$\bar{\mathcal{I}}_{\rm C}$ as a function of the dissipation strength
$\gamma$. (b) Averaged inverse participation ratio for $\bar{\mathcal{I}}$ as
a function of the dissipation strength $\gamma$ for various disorder strength
$W=12J$ (cyan circles), $W=7J$ (orange triangles), $W=5J$ (green squares) with
a common interaction strength $\Delta=0.5J$. The red dashed line marks the
finite size value of $\mathcal{I}$ for a Dicke state. Other parameters are
$L=10$ and $100$ disorder realizations.

When a disordered system is coupled to a bath, the localization properties
could be due to the fact the Hamiltonian is localized or that the dissipator
dominates the dynamics and drives the system towards a steady state which is
localized, or a combination of the two effects. To be able to characterize
whether the localization properties of a steady state are more prominently
inherited from the Hamiltonian, we introduce the concatenated inverse
participation ratio $\mathcal{I}_{\rm C}$ which measures whether a steady
state is, at the same time, diagonal in the energy basis of the Hamiltonian
and the energy eigenstates are localized. More precisely, we define
$\mathcal{I}_{\rm C}$ as

$\displaystyle\mathcal{I}_{\rm
C}(\hat{\rho}_{s})=\sum_{i,j}N_{i}|\langle\psi_{j}|\varphi_{i}\rangle|^{4}\mathcal{I}(\Ket{\psi_{j}}\Bra{\psi_{j}})$
(7)

where $\mathcal{I}(\Ket{\psi_{j}}\Bra{\psi_{j}})$ is the inverse participation
ratio of the natural orbitals of the $j$-th eigenstate of the Hamiltonian
$|\psi_{j}\rangle$ given in Eq. (5) Bera  _et al._ (2015). In Eq. (7),
$|\varphi_{i}\rangle$ is the $i$-th eigenstates of the steady-state density
operator,
$\hat{\rho}_{s}=\sum_{i}N_{i}|\varphi_{i}\rangle\langle\varphi_{i}|$. In Eq.
(7) $|\langle\psi_{j}|\varphi_{i}\rangle|^{4}$ can give a large contribution
only when $|\varphi_{i}\rangle$ has significant overlap with only a single
$|\psi_{j}\rangle$ 111We assume no degeneracies, and
$\mathcal{I}(\Ket{\psi_{j}}\Bra{\psi_{j}})$ is large if that eigenstate is
localized. Hence, $\mathcal{I}_{\rm C}$ can indeed inform us on whether the
steady state is made of an incoherent mixture of localized energy eigenstates.
In Fig. 5 we show both the averaged $\mathcal{I}_{\rm C}$ from Eq. (7),
$\bar{\mathcal{I}}_{\rm C}$ in panel (a), and the averaged $\mathcal{I}$ from
Eq. (5), $\bar{\mathcal{I}}$ in panel (b), both averaged over $100$ disorder
realizations. In both panels, we showed the inverse participation ratios for
different disorder strengths versus the magnitude of the dissipation rate
$\gamma$. We observe that for small $\gamma$, and large enough disorder, the
$\mathcal{I}_{\rm C}$ is large and it decreases as $\gamma$ increases. From
comparing $\mathcal{I}_{\rm C}$ with $\bar{\mathcal{I}}$, we observe that they
have a similar trend, indicating that the localization properties that we
observe in the steady state are, in a sense, inherited from the Hamiltonian.
In the limit of large $\gamma$ both inverse participation ratios obtain a
minimum. We note that $\bar{\mathcal{I}}$ tends towards a finite value which
may seem counter-intuitive given the delocalizing nature of the dissipator.
This can be explained by the fact that even for the Dicke state (which is the
state that the dissipator alone would produce),
$\mathcal{I}(|\psi_{D}\rangle\langle\psi_{D}|)$ reaches a small yet finite
value which we indicate, in Fig. 5(b) by the red dashed line. This asymptotic
value would decrease as the system size increases.

Figure 6: Single realizations of density matrix in spin configurations basis
(a, b) and in eigenbasis of the system (c, d) an intermediate dissipation
strength $\hbar\gamma=1J$ in left panels (a, c) and large dissipation strength
$\hbar\gamma=10J$ in right panels (c, d). Common parameters are $\Delta=0.5J$,
$W=12J$ and $L=10$.

The fact that in this model the dissipator tends to delocalize the system can
also be observed in Fig. 6 in which we show the modulus of the elements of the
density matrix both in the spin configuration, panels (a, b), and in the
energy basis, panels (c, d). In panels (a, c) we considered an intermediate
value of the dissipation rate $\hbar\gamma=J$, and in panels (b, d) we chose a
much larger value $\hbar\gamma=10J$. The latter case shows a significantly
increased coherence in the steady-state density matrix.

##  IV Conclusions

We have studied a disordered interacting spin chain coupled to an external
bath. While disorder and interaction can lead the system to be localized, the
bath tends to induce homogeneous states for the system. The rich variety of
attainable steady states in this system is exemplified in Figs. (1, 2), in
which we can observe that for small disorder there is large coherence even far
from the diagonal, and for $\Delta=J$ the steady state well approximates the
Dicke state. Larger values of disorder or interaction make the density matrix
close to diagonal both in the spin configuration basis and in the energy
eigenbasis of the system, especially after averaging out over many disorder
realizations. It is thus an interesting set-up to study disordered many-body
open quantum systems, as already shown in Refs. Vakulchyk  _et al._ , 2018;
Xu, Guo, and Poletti, 2018.

While Refs. Vakulchyk  _et al._ , 2018; Xu, Guo, and Poletti, 2018 focus on
possible signatures of many-body localization in open quantum systems, this
article focuses on the localization-delocalization properties within the open
quantum systems and gives an attempt to measure the localization for steady
state directly with respect to the system Hamiltonian. This is done by first
extending the studies of IPR by probing the phase diagram of the system in
steady states. We show that the averaged inverse participation ratio of the
natural orbitals of the single-particle density matrix can be large when the
disorder is large enough, both for weak and strong interactions.

We then revisit the meaning of the localization in many-body open quantum
systems by analyzing the steady-state density matrix in the Hamiltonian
eigenbasis and we find that for intermediate-large disorders the steady state
is composed of an exponentially decaying distribution of energy eigenstates,
each of which is localized. However, as disorder or interaction increases, the
distribution of eigenstates flattens out and more eigenstates have a
significant overlap with the steady state. This informs us that while the
inverse participation ratio can be very large, the steady state can be very
different, and only in regimes with intermediate-large disorders, and weak
interaction, we obtain a steady state in which a few energy eigenstates are
more prominent.

Last, we proposed the concatenated inverse participation ratio which takes
into account both the “diagonality” of the density matrix in energy eigenbasis
of the system Hamiltonian as well as the localization properties of the basis
states. This allows us to probe whether the signatures of localization are
intrinsic of the system Hamiltonian or more from an interplay of Hamiltonian
and dissipative terms. For the model studied here we observe that when the
inverse participation ratio is large, the steady state is also close to
diagonal in the energy basis, indicating that its localization properties are
inherited from the Hamiltonian. Increasing the dissipation rates increases the
coherence and decreases the localization properties of the steady state. In
the future, larger systems sizes should be studied, and other types of non-
thermal baths, so as to learn more about possible phases of matter emerging
from the interplay of interaction, disorder, and tailored dissipators. It
would also be interesting to study the relaxation dynamics of two-time
correlations and out-of-time-order correlations for many-body localized
systems Fan  _et al._ (2017) with dissipation.

###### Acknowledgements.

D.P. acknowledges support from the Ministry of Education of Singapore AcRF MOE
Tier-II (Project No. MOE2018-T2-2-142) and fruitful discussions with D.
Abanin, B. Grémaud, F. Heidrich-Meisner, and R. Chen. Part of the
computational work was performed on resources of the National Supercomputing
Centre, Singapore NSC .

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Anderson (1958) P. W. Anderson, “Absence of Diffusion in Certain Random Lattices,” Phys. Rev. 109, 1492–1505 (1958).
  * Basko, Aleiner, and Altshuler (2006) D. Basko, I. Aleiner, and B. Altshuler, “Metal–insulator transition in a weakly interacting many-electron system with localized single-particle states,” Ann. Phys. (N.Y.) 321, 1126–1205 (2006).
  * Oganesyan and Huse (2007) V. Oganesyan and D. A. Huse, “Localization of interacting fermions at high temperature,” Phys. Rev. B 75, 155111 (2007).
  * Pal and Huse (2010) A. Pal and D. A. Huse, “Many-body localization phase transition,” Phys. Rev. B 82, 174411 (2010).
  * Giamarchi and Schulz (1987) T. Giamarchi and H. J. Schulz, “Localization and Interaction in One-Dimensional Quantum Fluids,” Europhys. Lett. 3, 1287–1293 (1987).
  * Giamarchi and Schulz (1988) T. Giamarchi and H. J. Schulz, “Anderson localization and interactions in one-dimensional metals,” Phys. Rev. B 37, 325–340 (1988).
  * Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, “Many-Body Localization and Thermalization in Quantum Statistical Mechanics,” Annu. Rev. Condens. Matter Phys. 6, 15–38 (2015).
  * Abanin and Papić (2017) D. A. Abanin and Z. Papić, “Recent progress in many-body localization: Recent progress in many-body localization,” Ann. Phys. (Berlin) 529, 1700169 (2017).
  * Abanin  _et al._ (2019) D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, “Colloquium: Many-body localization, thermalization, and entanglement,” Rev. Mod. Phys. 91, 021001 (2019).
  * Schreiber  _et al._ (2015a) M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, “Observation of many-body localization of interacting fermions in a quasirandom optical lattice,” Science 349, 842–845 (2015a).
  * Schreiber  _et al._ (2015b) M. Schreiber, S. S. Hodgman, P. Bordia, H. P. Lüschen, M. H. Fischer, R. Vosk, E. Altman, U. Schneider, and I. Bloch, “Observation of many-body localization of interacting fermions in a quasirandom optical lattice,” Science 349, 842–845 (2015b).
  * Bordia  _et al._ (2016) P. Bordia, H. P. Lüschen, S. S. Hodgman, M. Schreiber, I. Bloch, and U. Schneider, “Coupling Identical one-dimensional Many-Body Localized Systems,” Phys. Rev. Lett. 116, 140401 (2016).
  * Choi  _et al._ (2016) J.-y. Choi, S. Hild, J. Zeiher, P. Schauß, A. Rubio-Abadal, T. Yefsah, V. Khemani, D. A. Huse, I. Bloch, and C. Gross, “Exploring the many-body localization transition in two dimensions,” Science 352, 1547–1552 (2016).
  * Bordia  _et al._ (2017a) P. Bordia, H. Lüschen, S. Scherg, S. Gopalakrishnan, M. Knap, U. Schneider, and I. Bloch, “Probing Slow Relaxation and Many-Body Localization in Two-Dimensional Quasiperiodic Systems,” Phys. Rev. X 7, 041047 (2017a).
  * Bordia  _et al._ (2017b) P. Bordia, H. Lüschen, U. Schneider, M. Knap, and I. Bloch, “Periodically driving a many-body localized quantum system,” Nat. Phys. 13, 460–464 (2017b).
  * Lüschen  _et al._ (2017) H. P. Lüschen, P. Bordia, S. S. Hodgman, M. Schreiber, S. Sarkar, A. J. Daley, M. H. Fischer, E. Altman, I. Bloch, and U. Schneider, “Signatures of Many-Body Localization in a Controlled Open Quantum System,” Phys. Rev. X 7, 011034 (2017).
  * Smith  _et al._ (2016) J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, “Many-body localization in a quantum simulator with programmable random disorder,” Nat. Phys. 12, 907–911 (2016).
  * Roushan  _et al._ (2017) P. Roushan, C. Neill, J. Tangpanitanon, V. M. Bastidas, A. Megrant, R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. Fowler, B. Foxen, M. Giustina, E. Jeffrey, J. Kelly, E. Lucero, J. Mutus, M. Neeley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. White, H. Neven, D. G. Angelakis, and J. Martinis, “Spectroscopic signatures of localization with interacting photons in superconducting qubits,” Science 358, 1175–1179 (2017).
  * Guo  _et al._ (2020) Q. Guo, C. Cheng, Z.-H. Sun, Z. Song, H. Li, Z. Wang, W. Ren, H. Dong, D. Zheng, Y.-R. Zhang, R. Mondaini, H. Fan, and H. Wang, “Observation of energy-resolved many-body localization,” Nat. Phys. , 1–6 (2020).
  * Nandkishore, Gopalakrishnan, and Huse (2014) R. Nandkishore, S. Gopalakrishnan, and D. A. Huse, “Spectral features of a many-body-localized system weakly coupled to a bath,” Phys. Rev. B 90, 064203 (2014).
  * Johri, Nandkishore, and Bhatt (2015) S. Johri, R. Nandkishore, and R. N. Bhatt, “Many-Body Localization in Imperfectly Isolated Quantum Systems,” Phys. Rev. Lett. 114, 117401 (2015).
  * Nandkishore (2015) R. Nandkishore, “Many-body localization proximity effect,” Phys. Rev. B 92, 245141 (2015).
  * Žnidarič, Scardicchio, and Varma (2016) M. Žnidarič, A. Scardicchio, and V. K. Varma, “Diffusive and subdiffusive spin transport in the ergodic phase of a many-body localizable system,” Phys. Rev. Lett. 117, 040601 (2016).
  * Medvedyeva, Prosen, and Žnidarič (2016) M. V. Medvedyeva, T. Prosen, and M. Žnidarič, “Influence of dephasing on many-body localization,” Phys. Rev. B 93, 094205 (2016).
  * Levi  _et al._ (2016) E. Levi, M. Heyl, I. Lesanovsky, and J. P. Garrahan, “Robustness of Many-Body Localization in the Presence of Dissipation,” Phys. Rev. Lett. 116, 237203 (2016).
  * Fischer, Maksymenko, and Altman (2016) M. H. Fischer, M. Maksymenko, and E. Altman, “Dynamics of a Many-Body-Localized System Coupled to a Bath,” Phys. Rev. Lett. 116, 160401 (2016).
  * Hyatt  _et al._ (2017) K. Hyatt, J. R. Garrison, A. C. Potter, and B. Bauer, “Many-body localization in the presence of a small bath,” Phys. Rev. B 95, 035132 (2017).
  * Everest  _et al._ (2017) B. Everest, I. Lesanovsky, J. P. Garrahan, and E. Levi, “Role of interactions in a dissipative many-body localized system,” Phys. Rev. B 95, 024310 (2017).
  * Luitz, Huveneers, and De Roeck (2017) D. J. Luitz, F. Huveneers, and W. De Roeck, “How a Small Quantum Bath Can Thermalize Long Localized Chains,” Phys. Rev. Lett. 119, 150602 (2017).
  * Nandkishore and Gopalakrishnan (2017) R. Nandkishore and S. Gopalakrishnan, “Many body localized systems weakly coupled to baths: Many body localized systems weakly coupled to baths,” Ann. Phys. (Berlin) 529, 1600181 (2017).
  * Žnidarič  _et al._ (2017) M. Žnidarič, J. J. Mendoza-Arenas, S. R. Clark, and J. Goold, “Dephasing enhanced spin transport in the ergodic phase of a many-body localizable system,” Ann. Phys. (Berlin) 529, 1600298 (2017).
  * Žnidarič and Ljubotina (2018) M. Žnidarič and M. Ljubotina, “Interaction instability of localization in quasiperiodic systems,” PNAS 115, 4595–4600 (2018).
  * Karlsson, Hopjan, and Verdozzi (2018) D. Karlsson, M. Hopjan, and C. Verdozzi, “Disorder and interactions in systems out of equilibrium: The exact independent-particle picture from density functional theory,” Phys. Rev. B 97, 125151 (2018).
  * Vakulchyk  _et al._ (2018) I. Vakulchyk, I. Yusipov, M. Ivanchenko, S. Flach, and S. Denisov, “Signatures of many-body localization in steady states of open quantum systems,” Phys. Rev. B 98, 020202 (2018).
  * Xu, Guo, and Poletti (2018) X. Xu, C. Guo, and D. Poletti, “Interplay of interaction and disorder in the steady state of an open quantum system,” Phys. Rev. B 97, 140201 (2018).
  * Wu  _et al._ (2019) L.-N. Wu, A. Schnell, G. D. Tomasi, M. Heyl, and A. Eckardt, “Describing many-body localized systems in thermal environments,” New J. Phys. 21, 063026 (2019).
  * Rubio-Abadal  _et al._ (2019) A. Rubio-Abadal, J.-y. Choi, J. Zeiher, S. Hollerith, J. Rui, I. Bloch, and C. Gross, “Many-Body Delocalization in the Presence of a Quantum Bath,” Phys. Rev. X 9, 041014 (2019).
  * Wybo, Knap, and Pollmann (2020) E. Wybo, M. Knap, and F. Pollmann, “Entanglement dynamics of a many-body localized system coupled to a bath,” Phys. Rev. B 102, 064304 (2020).
  * Diehl  _et al._ (2008) S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Büchler, and P. Zoller, “Quantum states and phases in driven open quantum systems with cold atoms,” Nat. Phys. 4, 878–883 (2008).
  * Schindler  _et al._ (2013) P. Schindler, M. Müller, D. Nigg, J. T. Barreiro, E. A. Martinez, M. Hennrich, T. Monz, S. Diehl, P. Zoller, and R. Blatt, “Quantum simulation of dynamical maps with trapped ions,” Nat. Phys. 9, 361–367 (2013).
  * Jordan and Wigner (1928) P. Jordan and E. Wigner, “Über das Paulische Äquivalenzverbot,” Z. Physik 47, 631–651 (1928).
  * Lieb, Schultz, and Mattis (1961) E. Lieb, T. Schultz, and D. Mattis, “Two soluble models of an antiferromagnetic chain,” Ann. Phys. (N.Y.) 16, 407–466 (1961).
  * Bera  _et al._ (2015) S. Bera, H. Schomerus, F. Heidrich-Meisner, and J. H. Bardarson, “Many-Body Localization Characterized from a One-Particle Perspective,” Phys. Rev. Lett. 115, 046603 (2015).
  * Bera  _et al._ (2017) S. Bera, T. Martynec, H. Schomerus, F. Heidrich-Meisner, and J. H. Bardarson, “One-particle density matrix characterization of many-body localization,” Ann. Phys. (Berlin) 529, 1600356 (2017).
  * Lezama  _et al._ (2017) T. L. M. Lezama, S. Bera, H. Schomerus, F. Heidrich-Meisner, and J. H. Bardarson, “One-particle density matrix occupation spectrum of many-body localized states after a global quench,” Phys. Rev. B 96, 060202 (2017).
  * Lin  _et al._ (2018) S.-H. Lin, B. Sbierski, F. Dorfner, C. Karrasch, and F. Heidrich-Meisner, “Many-body localization of spinless fermions with attractive interactions in one dimension,” SciPost Physics 4, 002 (2018).
  * Luitz, Laflorencie, and Alet (2015) D. J. Luitz, N. Laflorencie, and F. Alet, “Many-body localization edge in the random-field Heisenberg chain,” Phys. Rev. B 91, 081103 (2015).
  * Serbyn, Papić, and Abanin (2015) M. Serbyn, Z. Papić, and D. A. Abanin, “Criterion for Many-Body Localization-Delocalization Phase Transition,” Phys. Rev. X 5, 041047 (2015).
  * Šuntajs  _et al._ (2020) J. Šuntajs, J. Bonča, T. Prosen, and L. Vidmar, “Quantum chaos challenges many-body localization,” Phys. Rev. E 102, 062144 (2020).
  * Abanin  _et al._ (2021) D. A. Abanin, J. H. Bardarson, G. De Tomasi, S. Gopalakrishnan, V. Khemani, S. A. Parameswaran, F. Pollmann, A. C. Potter, M. Serbyn, and R. Vasseur, “Distinguishing localization from chaos: Challenges in finite-size systems,” Annals of Physics , 168415 (2021).
  * Panda  _et al._ (2019) R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor, and M. Žnidarič, “Can we study the many-body localisation transition?” Europhys. Lett. 128, 67003 (2019).
  * Kiefer-Emmanouilidis  _et al._ (2021) M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, “Slow delocalization of particles in many-body localized phases,” Phys. Rev. B 103, 024203 (2021).
  * Gorini, Kossakowski, and Sudarshan (1976) V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, “Completely positive dynamical semigroups of N-level systems,” J. Math. Phys. 17, 821–825 (1976).
  * Lindblad (1976) G. Lindblad, “On the generators of quantum dynamical semigroups,” Commun. Math. Phys. 48, 119–130 (1976).
  * Diehl  _et al._ (2010) S. Diehl, A. Tomadin, A. Micheli, R. Fazio, and P. Zoller, “Dynamical Phase Transitions and Instabilities in Open Atomic Many-Body Systems,” Phys. Rev. Lett. 105, 015702 (2010).
  * Tomadin, Diehl, and Zoller (2011) A. Tomadin, S. Diehl, and P. Zoller, “Nonequilibrium phase diagram of a driven and dissipative many-body system,” Phys. Rev. A 83, 013611 (2011).
  * Yi  _et al._ (2012) W. Yi, S. Diehl, A. J. Daley, and P. Zoller, “Driven-dissipative many-body pairing states for cold fermionic atoms in an optical lattice,” New J. Phys. 14, 055002 (2012).
  * Bardyn  _et al._ (2013) C.-E. Bardyn, M. A. Baranov, C. V. Kraus, E. Rico, A. \\.Imamoğlu, P. Zoller, and S. Diehl, “Topology by dissipation,” New J. Phys. 15, 085001 (2013).
  * Lloyd (1996) S. Lloyd, “Universal Quantum Simulators,” Science 273, 1073–1078 (1996).
  * Yusipov  _et al._ (2017) I. Yusipov, T. Laptyeva, S. Denisov, and M. Ivanchenko, “Localization in Open Quantum Systems,” Phys. Rev. Lett. 118, 070402 (2017).
  * Vershinina  _et al._ (2017) O. S. Vershinina, I. I. Yusipov, S. Denisov, M. V. Ivanchenko, and T. V. Laptyeva, “Control of a single-particle localization in open quantum systems,” Europhys. Lett. 119, 56001 (2017).
  * Note (1) We assume no degeneracies.
  * Fan  _et al._ (2017) R. Fan, P. Zhang, H. Shen, and H. Zhai, “Out-of-time-order correlation for many-body localization,” Science Bulletin 62, 707–711 (2017).
  * (64) https://www.nscc.sg/.