Hardy perturbations of subordinated Bessel heat kernels
Krzysztof Bogdan Tomasz Jakubowski Konstantin Merz


Heat kernel bounds — ]Hardy perturbations of 
 subordinated Bessel heat kernels


K. Bogdan]Krzysztof Bogdan
[Krzysztof Bogdan]Department of Pure and Applied Mathematics, Wrocław University of Science and Technology, Hoene-Wrońskiego 13C, 50-376 Wrocław, Poland
krzysztof.bogdan@pwr.edu.pl

T. Jakubowski]Tomasz Jakubowski
[Tomasz Jakubowski]Department of Pure and Applied Mathematics, Wrocław University of Science and Technology, Hoene-Wrońskiego 13C, 50-376 Wrocław, Poland
tomasz.jakubowski@pwr.edu.pl

K. Merz]Konstantin Merz
[Konstantin Merz]Institut für Analysis und Algebra, Technische Universität Braunschweig, Universitätsplatz 2, 38106 Braun­schweig, Germany 
k.merz@tu-bs.de

 

[2020]Primary 47D08, 60J35






§ ABSTRACT 
  Motivated by the spectral theory of relativistic atoms,  we prove matching upper and lower bounds for the transition density of Hardy perturbations of subordinated Bessel heat kernels. The analysis is based on suitable supermedian functions, in particular invariant functions.


This research was supported in part by grant Opus 2023/51/B/ST1/02209 of National Science Center, Poland (K.B.), and the PRIME programme of the German Academic Exchange Service DAAD with funds from the German Federal Ministry of Education and Research BMBF (K.M.).


[
    [
    
=====








§ INTRODUCTION AND MAIN RESULT



Bessel operators are fundamental for mathematical analysis, theoretical physics, and applied sciences. We are interested in their fractional powers with added critically singular Hardy potential,

    (-d^2dr^2-2ζr ddr)^α/2-κ r^-α

in L^2(_+,r^2ζdr) with _+=(0,∞).
Here and in what follows, 

    ζ∈(-1/2,∞),
        α∈(0,2]∩(0,2ζ+1),

and κ∈ is called the coupling constant of the Hardy potential.
Below, we call (the potential) κ r^-α attractive (mass-creating) if κ>0 and repulsive (mass-killing) if κ<0.
Our interest in the operators (<ref>) comes from the spectral theory of atoms. 
We pursue a program initiated in <cit.> to describe the different rates of spectral dissipation of mass in various angular momentum channels in , as parameterized by different values of ζ, α, and κ. From the point of view of theoretical physics, the most important is the three-dimensional pseudo-relativistic case  α=1, with ζ=1+ℓ,  ℓ=0,1,…, and Coulomb potential κ r^-1. Here,  the operator (<ref>) arises as one of the direct summands of the direct sum decomposition of (-Δ)^1/2-κ/|x| in L^2(^3) into different angular momentum channels, indexed by ℓ∈_0. Here,
κ defines the strength of coupling  between the nucleus and electrons in an atom. The passage from the state space ^3, or, more generally,  with d∈, to the half-line _+ in (<ref>) is inspired by the intertwining in <cit.>. In fact, in the present paper, we abstract from the setting of <cit.> and study, in maximal generality, a semigroup on _+ associated with (<ref>).
The spectral analysis of <cit.> in  will be picked up again in the forthcoming paper <cit.>, where we also clarify the definition of (<ref>) as a self-adjoint operator.
The main result of the present paper, Theorem <ref> below, provides a crucial technical ingredient to be used in <cit.>, namely non-explosion and sharp heat kernel bounds for a semigroup associated with (<ref>).




To state the result, we briefly introduce necessary notation, but defer precise definitions to the ensuing subsections. In particular, here and below, we use a specific parameterization of the coupling constant κ in (<ref>). Thus, for ζ∈(-1/2,∞), we let 

    Ψ_ζ(η)κ=Ψ_ζ(η) := 
        2^αΓ(2ζ+1-η/2)Γ(α+η/2)/Γ(η/2)Γ(2ζ+1-η-α/2)    if α∈(0,2∧ (2ζ+1)) and η∈(-α,2ζ+1),
          
    
     
          (2ζ-1-η)η    if α=2<2ζ+1 and η∈.

Here, as usual, A∧ B:=min{A,B} and A∨ B:=max{A,B}∧∨.
We introduce further notation as we proceed; it is summarized in the [index] at the end of the paper.
Clearly, Ψ_ζ(η) is symmetric about η=(2ζ+1-α)/2.
In view of <cit.>,  its maximum is 

    κ_ c(ζ,α):=Ψ_ζ((2ζ+1-α)/2).

For instance, in the three-dimensional, pseudorelativistic situation with ζ=1+ℓ, 

    κ_ c(1+ℓ,1) = 2Γ(1+ℓ/2)^2/Γ((ℓ+1)/2)^2.

In particular, for ℓ=0, we get κ_ c(1,1)=2/π while for ℓ=1, we get κ_ c(2,1)=π/2.
In the general situation, Ψ_ζ(η) is strictly positive for η∈(0,2ζ+1-α), zero for η∈{0,2ζ+1-α}, and tends to -∞ as η→ -M or η→2ζ+1-α+M. Here M:=α if α<2 and M:=∞ for α=2. Thus, the parameterization κ=Ψ_ζ(η) only produces
κ∈ (-∞,κ_ c(ζ,α)].
For this range,
we may and often do restrict η to the interval (-M,(2ζ+1-α)/2].
Moreover, for ζ>ζ'>-1/2, we have Ψ_ζ(η)>Ψ_ζ'(η) for all 
η in the intersection of the domains of Ψ_ζ and Ψ_ζ'.
By the Hardy inequality, in fact, by the ground state representation <cit.>, the quadratic form corresponding to (<ref>) on L^2(_+,r^2ζ dr) is bounded from below if κ=Ψ_ζ(η) for some η
and it is unbounded from below if κ>κ_ c(ζ,α). 
In the former case, the form is bounded from below by zero. This observation motivates to call κ_ c(ζ,α) the critical coupling constant.

Summarizing, in most of the discussion below, (<ref>) is substituted by

    (-d^2dr^2-2ζr ddr)^α/2-Ψ_ζ(η) r^-α,

for η admissible, that is, in the domain of Ψ_ζ, see (<ref>).
We call 2ζ+1 the effective dimension.
This is because the reference measure for the Hilbert space in which (<ref>) is defined is r^2ζdr=r^2ζ+1dr/r, where dr/r is the Haar measure on the multiplicative group (0,∞).
As we shall see, η conveniently parametrizes Hardy potentials, and, of course, α/2 is the order of subordination. All the operators (<ref>) are homogeneous of order α (with respect to dilations), but
the interplay of the three parameters critically influences their potential-theoretic properties.
We let

    p_ζ^(α)(t,r,s),    t, r,s>0,

be the transition density associated to (<ref>) with κ=0; see (<ref>) and (<ref>) for actual definitions in the cases α=2 and α<2, respectively. Here, t>0 is time, but r,s∈ (0,∞) are positions in space.
Moreover, we let

    p_ζ,η^(α)(t,r,s),    t, r,s>0,

be the Schrödinger perturbation of p_ζ^(α) by the Hardy potential Ψ_ζ(η)/r^α, 
as in Definition <ref>.
In <cit.>, we will show that p_ζ,η^(α) equals the heat kernel of (<ref>) in the sense of quadratic forms.
The main contribution of the present paper are the following matching upper and lower, or sharp, bounds for p_ζ,η^(α), which are proved in Section <ref>.

  
  If ζ∈(-1/2,∞), α∈(0,2)∩(0,2ζ+1), and η∈(-α,2ζ+1-α/2], then
  
    p_ζ,η^(α)(t,r,s)
          ∼_ζ,α,η(1∧r/t^1/α)^-η (1∧s/t^1/α)^-η  p_ζ^(α)(t,r,s),

  and p_ζ,η^(α)(t,r,s) 
is a jointly continuous function of t,r,s>0.

Here, we write A ∼ B∼, or B ≲ A ≲ B, if A, B ≥ 0 and c^-1B ≤ A ≤ cB for a constant c ∈ (0,∞). We write A ∼_τ B if c = c_τ depends on τ, etc. The dependence on fixed parameters, like ζ,α,η is usually omitted.

In view of the sharp Hardy inequality in <cit.>, one may duly expect a blowup for supercritical coupling constants κ, i.e., κ>κ_ c(ζ,α). 
Indeed, the interpretation and proof of the following corollary to Theorem <ref> are given in Subsection <ref> below.

  
  Let ζ∈(-1/2,∞) and α∈(0,2). Then  the heat kernel associated with (<ref>) blows up if κ >κ_ c(ζ,α).


 
  
  For α=2, we have the identity
  
    p_ζ,η^(2)(t,r,s) = (rs)^-ηp_ζ-η^(2)(t,r,s).

  It is proved, e.g., by Metafune, Negro, and Spina <cit.> by a change variables; see also the beginning of Section <ref> for a sketch of proof. For α=2, (<ref>) and  (<ref>) yield that p_ζ,η^(2)(t,r,s) is jointly continuous on (0,∞)^3 and obeys (<ref>), provided ∼ is replaced by the notation ≍ explained in the paragraph below (<ref>). Moreover, the blow up result in Corollary <ref> extends to α=2.




Below in this paper we develop  potential-theoretic techniques, which we refer to as the integral method, to study kernels appearing in <cit.> and <cit.>. This may be considered as a complementary study largely independent of <cit.> and <cit.>. Namely, we abstract from the setting of  and the angular momentum channels parameterized by ℓ in <cit.> by projecting, transferring, or intertwining <cit.> onto the half-line _+. This allows us to use sub-Markovian semigroups to capture delicate oscillations that influence the different rates of dissipation of mass for different spherical harmonics. This method is not only relevant for this context but is also of independent interest. The addition of the Hardy potential has an analogous effect on the heat kernel, changing the dissipation rate but preserving scaling.

In fact, additional motivation for our work 
comes from the fact that models with scaling are of fundamental interest for PDEs, functional analysis, probability, and physics; for further connections, see, e.g., 
<cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>.

The main methodological inspiration for this work comes from Bogdan, Grzywny, Jakubowski, and Pilarczyk <cit.> and Jakubowski and Wang <cit.>. However, the present development differs in several important ways and is more complex. Namely, the projection onto ℝ_+ unfortunately eliminates both the origin and rotational invariance from the setting. Furthermore, we address η > 0 and η < 0 simultaneously. Moreover, the integral method for negative perturbations involves rather delicate compensation. 
Finally, our resolution is expected to have significant consequences for the multi-channeled spectral resolution of the Laplacian and related operators.

  
Let us note that Cho, Kim, Song, and Vondraček <cit.> (see also Song, Wu, and Wu <cit.>) and Jakubowski and Maciocha <cit.> recently proved heat kernel bounds for the fractional Laplacian with Hardy potential on _+ and Dirichlet boundary condition on (_+)^c. As in our analysis, <cit.> use the integral method introduced in  <cit.>, but, in contrast to our approach, they cannot rely on subordination since the Dirichlet fractional Laplacian is 
different from the fractional Dirichlet Laplacian, see, e.g., <cit.>. 
Consequently, the results of <cit.> and <cit.> are different than ours, as can be readily seen by comparing the respective parameterizations of the coupling constants, cf. <cit.> and (<ref>).


In the following subsections, we make our setting more precise. To this end, we first define the subordinated Bessel heat kernels p_ζ^(α) and recall some of their properties. 
Then  we give a precise definition of p_ζ,η^(α).
At the end of this section, we outline implications of Theorem <ref> and the structure of the rest of the paper.
Recall that in <cit.>, we will prove that p_ζ,η^(α) equals the heat kernel of (<ref>) in the sense of quadratic forms for κ=Ψ_ζ(η), but we mention this fact only for motivation and of course do not use it in the present paper.



 §.§ Subordinated Bessel heat kernels


We define the Bessel heat kernel

    p_ζ^(2)(t,r,s)
        p_ζ^(2)(t,r,s)     : = (rs)^1/2-ζ/2texp(-r^2+s^2/4t)I_ζ-1/2(rs/2t),    r,s,t>0.

Here, for z∈∖(-∞,0], I_ν(z) denotes the modified Bessel functions of the first kind of order ν∈ <cit.>.I_ν
Note that

    p_ζ^(2)(t,r,0)
      = 2^-2ζ/Γ(2ζ+1/2) t^-2ζ+1/2exp(-r^2/4t)

by the series expansion <cit.> of the modified Bessel function.

The kernel p_ζ^(2)(t,r,s) with the reference (speed) measure r^2ζdr on _+ is the transition density of the Bessel process with index ζ-1/2 reflected at the origin. We remark that the Bessel process with index ζ-1/2 killed at the origin has the transition density (<ref>), but with I_ζ-1/2(·) replaced with I_|ζ-1/2|(·). When ζ≥1/2, the Bessel process never hits the origin, i.e., no condition (reflecting or killing) takes effect or needs to be imposed at the origin. See, e.g., <cit.>. We also remark that p_ζ^(2)(t,r,s) is the heat kernel of the nonnegative Bessel operator of index ζ-1/2 in the Liouville form,

    _ζ = -d^2/dr^2 - 2ζ/rd/dr  in L^2(_+,r^2ζdr).

More precisely, we understand _ζ as the Friedrichs extension of the corresponding quadratic form on C_c^∞(_+) if ζ∈[1/2,∞), but we understand _ζ as the Krein extension of the corresponding quadratic form on C_c^∞(_+) if ζ∈(-1/2,1/2]. In both cases, _ζ is nonnegative and in particular self-adjoint. We refer, e.g., to <cit.>, <cit.>, or <cit.> for further details and references regarding the spectral theory of _ζ, in particular for operator and quadratic form domains and (form) cores. In that setting,

    p_ζ^(2)(t,r,s) = -t_ζ(r,s),    t,r,s>0.

This information is important to prove Theorem <ref> when α=2.

Next, we define, for α∈(0,2), the α2-subordinated Bessel heat kernels. To that end, recall that for α∈(0,2) and t>0, by Bernstein's theorem, the completely monotone function [0,∞)∋λ↦-tλ^α/2  is the Laplace transform of a probability density function _+∋τ↦σ_t^(α/2)(τ). That is,

    -tλ^α/2 = ∫_0^∞ dτ σ_t^(α/2)(τ) -τλ,    t>0, λ≥0,

see, e.g., <cit.>.
We give some useful properties of and sharp estimates for σ_t^(α/2)(τ), and references in <cit.>. Using (<ref>), we define the α2-subordinated Bessel heat kernel with reference measure r^2ζdr on _+ as

    p_ζ^(α)(t,r,s)
        p_ζ^(α)(t,r,s)     : = ∫_0^∞ dτ σ_t^(α/2)(τ)  p_ζ^(2)(τ,r,s),    r,s,t>0.


Crucially, p_ζ^(α)(t,r,s) is a probability transition density. In fact, it also defines a strongly continuous contraction semigroup on L^2(_+,r^2ζdr). We record this in the following proposition, which is proved, e.g., in <cit.>.


  
  Let ζ∈(-1/2,∞), α∈(0,2], and t,t',r,s>0. Then  p_ζ^(α)(t,r,s)>0,
  
    ∫_0^∞ ds  s^2ζ p_ζ^(α)(t,r,s) = 1, 
       ∫_0^∞ dz  z^2ζ  p_ζ^(α)(t,r,z) p_ζ^(α)(t',z,s) = p_ζ^(α)(t+t',r,s), 
        p_ζ^(α)(t,r,s) = t^-2ζ+1/α p_ζ^(α)(1,r/t^1/α,s/t^1/α),

  and {p_ζ(t,·,·)}_t>0 is a strongly continuous contraction semigroup on L^2(_+,r^2ζdr).





 §.§ Hardy-type Schrödinger perturbations



Let ζ∈(-1/2,∞), α∈(0,2]∩(0,2ζ+1), and recall that M=α when α∈(0,2) and M=∞ when α=2. For η∈(-M,2ζ+1), we let

    q(z)
      q(z) := Ψ_ζ(η)/z^α,    z>0.

Because of symmetry of Ψ_ζ(η) about η=(2ζ+1-α)/2, as a rule we assume η∈(-M,(2ζ+1-α)/2] when considering q.
In  the Appendix <ref>, for the convenience of the reader, we discuss rather general Schrödinger perturbations of transition densities. The focus of this paper is, however, on the following cases.


  
  For ζ∈(-1/2,∞), α∈(0,2)∩(0,2ζ+1), and η∈ (-α,2ζ+1-α/2], let
  p_ζ,η^(α) be
  the Schrödinger perturbation of p_ζ^(α) by q in (<ref>), defined in the Appendix <ref>.
  For α=2 with α<2ζ+1, and η∈(-∞,2ζ-1/2], let p_ζ,η^(2)(t,r,s) be the heat kernel of the self-adjoint operator _ζ-Ψ_ζ(η)/r^2, with the domain described in <cit.>.

Thus, in the case α=2, we typically refer to the literature. For α∈ (0,2), we rely on the discussion
 in the Appendix <ref>: a Schrödinger perturbation adds mass to the semigroup if q>0 and decreases mass if q<0. Thus, for η=0, we have q=0 and p_ζ,0^(α)=p_ζ^(α). For η∈(0,2ζ+1-α/2], we have q> 0 and

    p_ζ,η^(α)
        p_ζ,η^(α)(t,r,s)
            := ∑_n≥0 p_t^(n,D)(r,s),   where  
          p_t^(0,D)(r,s) := p_ζ^(α)(t,r,s),   and
    
        p_t^(n,D)(r,s)     := ∫_0^t dτ∫_0^∞ dz  z^2ζ p_ζ^(α)(t,r,z) q(z) p_t-τ^(n-1,D)(z,s),    n∈.

Each term p_t^(n,D)(r,s) may be understood as an iteration of Duhamel's formula below,
hence the superscript D.
Of course, we have the following bound

    p_ζ,η^(α)(t,r,s) ≥ p_ζ^(α)(t,r,s)   for η∈(0,2ζ+1-α/2].

For η∈(-α,0), we have q<0 and p_ζ,η^(α)≤ p_ζ^(α); see the Appendix <ref>, in particular Theorem <ref>, for the  construction and properties of negative Schrödinger perturbations of transition densities using the more comfortable and general time-inhomogeneous setting. See the remarks before Theorem <ref>, too.
We also note that to deal with η<0 in Subsection <ref>, 
we use the results on η>0 from Subsection <ref>. The structure of Section <ref> reflects this connection.

We also note that p_ζ,η^(α) may also be defined by the Feynman–Kac formula (for bridges), but the singularity of z^-α in (<ref>) makes this problematic for large negative Ψ_ζ(η), i.e., η close to -α. 
Actually, a viable (probabilistic) approach to negative η, based on the Feynman-Kac formula and approximation by the killed semigroup, is suggested by <cit.>, but we prefer to keep the paper more analytic, and propose the potential-theoretic approach of  the Appendix <ref>. 

In the following, we write h_β(r):=r^-βh_β(r) for r>0 and β∈, and abbreviate 
    h:=h_η.
h(r)
We record some fundamental properties of p_ζ,η^(α).


  
  Let ζ∈(-1/2,∞), α∈(0,2]∩(0,2ζ+1), η∈(-M,2ζ+1-α/2], and r,s,t,t'>0. Then the following statements hold. 

  
  We have the Chapman–Kolmogorov equation,
  
    ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(t,r,z)   p_ζ,η^(α)(t',z,s)
            = p_ζ,η^(α)(t+t',r,s).

  
  We have the Duhamel formula
  
    p_ζ,η^(α)(t,r,s)
              = p_ζ^(α)(t,r,s) + ∫_0^t dτ∫_0^∞ dz  z^2ζ p_ζ^(α)(τ,r,z)   q(z)   p_ζ,η^(α)(t-τ,z,s) 
        = p_ζ^(α)(t,r,s) + ∫_0^t dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(τ,r,z)   q(z)   p_ζ^(α)(t-τ,z,s).

  
  We have the scaling relations p_t^(n,D)(r,s)=t^-2ζ+1/αp_1^(n,D)(r/t^1/α,s/t^1/α) and
  
    p_ζ,η^(α)(t,r,s)
        = t^-2ζ+1/α p_ζ,η^(α)(1,r/t^1/α,s/t^1/α).

  
  Let η∈[0,2ζ+1-α/2]. Then the function h(r)=r^-η is supermedian with respect to p_ζ,η^(α), i.e.,
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t,r,s) h(s)
             ≤ h(r).

  In particular, p_ζ,η^(α)(t,r,s)<∞ for all t,r>0 and almost all s>0.



  For α=2, these claims follow by using the explicit expression in (<ref>).
  Thus, suppose α∈(0,2) from now on.
  The scaling (<ref>) follows from (<ref>), (<ref>) when η≥0 and Theorem <ref> when η<0, and induction.
  For η≥0, the other three properties follow, e.g., from <cit.> and <cit.> together with the computation (<ref>).
  For η<0, the Chapman–Kolmogorov equation and the Duhamel formula follow by the construction in Theorem <ref>.





 §.§ Implications of Theorem <ref>



In the pioneering paper <cit.>, Bogdan, Grzywny, Jakubowski, and Pilarczyk proved estimates for the heat kernel of the homogeneous Schrödinger operator (-Δ)^α/2-κ/|x|^α on ^d, sometimes called Hardy operator, with α∈(0,2∧ d) and 0≤κ≤κ^*,
where

    κ^*=2^αΓ((d+α)/4)^2 /Γ((d-α)/4)^2 
       = Ψ_d-1/2(d-α/2).

Among others, their bounds were crucial for the systematic study of the L^p-Sobolev spaces generated by powers of (-Δ)^α/2-κ/|x|^α in <cit.>—see the founding paper <cit.> dealing with the case α=2—and to prove the strong Scott conjecture concerning the electron distribution of large relativistic atoms <cit.>.

In <cit.>, we made a first step into a more detailed analysis of the Hardy operator by taking its spherical symmetry into account. More precisely, for d∈, ℓ∈ L_d⊆_0, and m∈ M_d,ℓ⊆^d-2 (for d≥3)[Explicit descriptions of L_d and M_d,ℓ are not important here, but can be found in <cit.>.], let Y_ℓ,m denote the orthonormal basis of (hyper)spherical harmonics in L^2(^d-1). In <cit.>, we considered the function space 

    V_ℓ,m:={u(|x|)|x|^ℓ Y_ℓ,m(x/|x|):  u∈ L^2(_+,r^d+2ℓ-1dr)}

of functions with given angular momentum ℓ.
For [u]_ℓ,m(x):=u(|x|)|x|^ℓ Y_ℓ,m(x/|x|) with u∈ L^2(_+,r^d-1+2ℓdr), we showed

    ⟨[u]_ℓ,m,((-Δ)^α/2-κ/|x|^α)[u]_ℓ,m⟩_L^2(^d)
          = ⟨ u,((-d^2/dr^2-2ζ/rd/dr)^α/2-κ/r^α)u⟩_L^2(_+,r^d+2ℓ-1 dr),

whenever α<d+2ℓ, κ=Ψ_d-1+2ℓ/2(η) and all η∈(-α,d+2ℓ). In particular, the two sides of (<ref>) are nonnegative if and only if κ≤Ψ_d-1+2ℓ/2((d+2ℓ-α)/2).

Moreover, we constructed the generalized ground states and proved a ground state representation or Hardy identity for the quadratic forms (<ref>), see <cit.>. This refined the ground state representation of (-Δ)^α/2-κ/|x|^α in L^2(^d), first proved in <cit.> by Fourier transform[For α=2, the ground state representation of -Δ-κ/|x|^2 was known earlier; see, e.g., <cit.>.]; see also <cit.> for a more abstract approach.
These generalized ground states are just h(r)=r^-η and are crucial for our estimates of 
p_ζ,η^(α)(t,r,s).
The main result of our forthcoming paper <cit.> states that p_ζ,η^(α)(t,r,s) is indeed the heat kernel of the restriction of (-Δ)^α/2-Ψ_d+2ℓ-1/2(η)/|x|^α to V_ℓ,m.

Let us now comment on open problems and connections.
The above discussion indicates that the estimates in Theorem <ref>
will be important to advance the study of L^p-Sobolev spaces generated by powers of Hardy operators.
Our bounds should also lead to limits at the origin and self-similar solutions for the considered kernels, see <cit.>.
Extensions to other homogeneous operators, including those with gradient perturbations, are of interest, too, see Metafune, Negro, and Spina <cit.>
and Kinzebulatov, Semënov, and Szczypkowski <cit.>.

Moreover, Theorem <ref> and the results in our forthcoming paper <cit.> are paramount for a more precise description of the location of electrons in large, relativistically described atoms. Let us explain the connection and future plans in more detail. Frank, Merz, Siedentop, and Simon <cit.> consider large atoms by accounting for relativistic effects close to the nucleus using the so-called Chandrasekhar operator. The main result of <cit.> is that, in the limit where the number of electrons and the nuclear charge goes to infinity, the ground state density, i.e., the probability density of finding one electron, converges to the so-called hydrogenic density, i.e., the sum of absolute squares of the eigenfunctions of √(-Δ+1)-1-κ/|x|, where κ is the (rescaled) interaction strength between the nucleus and the electrons. In fact, <cit.> shows a finer result stating that the ground state density of electrons with prescribed angular momentum ℓ∈_0 converges to the sum of absolute squares of the eigenfunctions of the operator

    (-d^2dr^2-2ζr ddr+1)^α/2-1-κ r^-α  in L^2(_+,r^2ζdr),

with α=1 and ζ=ℓ+1/2. The operator in (<ref>) is a bounded perturbation of (<ref>), but has, unlike (<ref>), infinitely many negative eigenvalues[This is because of the summand +1 in (-d^2dr^2-2ζr ddr+1)^α/2, which breaks the homogeneity and necessitates a distinction between large and small frequencies, the latter being responsible for the emergence of eigenvalues.].
Moreover, <cit.> provides pointwise estimates for the hydrogenic density using the heat kernel estimates for the homogeneous operator (-Δ)^α/2-κ/|x|^α established  in <cit.>.
 
In this connection, we note that the estimates in the present paper and in <cit.>, together with the arguments developed in <cit.>, will lead to bounds for the sum of absolute squares of eigenfunctions (<ref>). In particular, we can show that for all t>0, there is c_t>0 such that

    ∑_n≥1|ϕ_n(r)|^2
      ≤ c_t r^-2η,    0<r≤ t,

where η is in the one-to-one correspondence to κ given by (<ref>) and 
{ϕ_n}_n∈ are the eigenfunctions corresponding to the negative eigenvalues of the operator in (<ref>). We are optimistic that, for ζ and α fixed, the upper bound in (<ref>) is optimal. 
In fact, the short-distance upper bound for the hydrogenic density of √(1-Δ)-1-κ/|x|, obtained in <cit.>, is optimal due to the lower bound on the ground state of √(1-Δ)-1-κ/|x|, which was proved by Jakubowski, Kaleta, and Szczypkowski in <cit.> by using lower bounds for the corresponding heat kernel in <cit.>. While upper heat kernel bounds for (<ref>) follow from those in the present paper, <cit.>, and subordination, proving corresponding lower bounds will require different techniques.






 §.§ Structure of the paper



In Section <ref>, we state further pointwise bounds and a 3G inequality for the subordinated Bessel heat kernel p_ζ^(α).
In Section <ref>, we prove that the function h(r)=r^-η is invariant under p_ζ,η^(α). We also state further integral equalities and inequalities involving p_ζ,η^(α) and h(r).
In Section <ref>, we prove Theorem <ref> and Corollary <ref>.
Finally, in  the Appendix <ref>, we provide an abstract discussion of Schrödinger perturbations of transition densities, which may be of independent interest.



  §.§.§ Acknowledgments

K.M. thanks Haruya Mizutani for his hospitality at Osaka University, where parts of this research were carried out.
We also thank Volker Bach, Rupert Frank, Wolfhard Hansen, Jacek Małecki, Haruya Mizutani, and Grzegorz Serafin for helpful discussions.









§ BOUNDS FOR SUBORDINATED BESSEL HEAT KERNELS





 §.§ Pointwise bounds and explicit expressions


For s,r ∈ℝ, t>0, and α∈(0,2], let p^(α)(t,r,s) be the stable density on the real line, i.e., 

    p^(α)(t,r,s) 
      := 1/2π∫_-∞^∞-i(s-r)z-t|z|^α dz.

Note that

    p^(2)(t,r,s) = 1/√(4π t)exp(-(r-s)^2/4t)

is the Gauß–Weierstraß  kernel and, by <cit.>, 

    p^(α)(t,r,s) ∼t/(t^1/α + |r-s|)^1+α  for α∈(0,2).



From <cit.> and <cit.>, we record the following sharp upper and lower bounds for p_ζ^(α)(t,r,s). 


  
  Let ζ∈(-1/2,∞). Then, there are c,c'>0 such that for all r,s,t>0,
  
    
    
    p_ζ^(2)(t,r,s)
             ≍_ζ t^-1/2exp(-(r-s)^2/c t)/(rs+t)^ζ
       ≍_ζ(1∧r/t^1/2)^ζ (1∧s/t^1/2)^ζ (1/rs)^ζ· t^-1/2·exp(-(r-s)^2/c' t).

  ≍Moreover, for all α∈(0,2) and all r,s,t>0,
  
    p_ζ^(α)(t,r,s)
           ∼_ζ,αt/|r-s|^1+α(r+s)^2ζ + t^1+α/α(t^1/α+r+s)^2ζ
       ∼_ζ,α p^(α)(t,r,s) (t^1/α + r+s)^-2ζ.





Here and below the notation ≍_ζ combines an upper bound and a lower bound similarly as ∼_ζ, but the constants in exponential factors (e.g., the constants c and c' in (<ref>)–(<ref>)) may be different in the upper and the lower bounds. Note that we allow the constants in the exponential factors to depend on ζ. Then, for instance, (<ref>) is equivalent to the statement that there are positive c_j,ζ, j∈{1,2,3,4}, with

    c_1,ζ t^-1/2exp(-(r-s)^2/c_2,ζ t)/(rs+t)^ζ≤ p_ζ^(2)(t,r,s)
        ≤ c_3,ζ t^-1/2exp(-(r-s)^2/c_4,ζ t)/(rs+t)^ζ.


For an explicit expressions for p_ζ^(α)(t,r,s) in the physically important case α=1, we refer, e.g., to Betancor, Harboure, Nowak, and Viviani <cit.>.


Below, we will also use the following convenient estimate.


  
  Let ζ∈(-1/2,∞) and α∈(0,2). Then, for every σ,λ∈[0,2ζ+1] such that σ+λ∈[0,2ζ+1], we have
  
    p_ζ^(α)(τ,z,s) ≲τ^1/α+z+s/τ^1/α+|z-s|τ^-λ/α s^-σ z^σ+λ-2ζ-1.

  In particular if z<s/2 or s  ≲τ^1/α z, then
  
    p_ζ^(α)(τ,z,s) ≲τ^-λ/α s^-σ z^σ+λ-2ζ-1.




  Since 0<σ≤ 2ζ+1, by (<ref>) and (<ref>), we have
  
    p_ζ^(α)(τ,z,s) 
           ∼τ/(τ^1/α+|z-s|)^1+α (τ^1/α+z+s)^-2ζ
        = τ(τ^1/α+z+s)/(τ^1/α+|z-s|)^1+α  (τ^1/α+z+s)^-λ(τ^1/α+z+s)^-σ(τ^1/α+z+s)^σ+λ-2ζ-1
       ≤τ^1/α+z+s/τ^1/α+|z-s|τ^-λ/α s^-σ z^σ+λ-2ζ-1.

  For z<s/2 or s/2 ≤τ^1/α z, we have τ^1/α+z+s/τ^1/α+|z-s|≲ 1 and consequently p_ζ^(α)(τ,z,s) ≲τ^-λ/α s^-σ z^σ+λ-2ζ-1.




 §.§ Comparability results for p_ζ^(α)


The following lemma is proved in <cit.> and allows to compare two kernels p_ζ^(α) at different times and positions with each other.


  
  Let ζ∈(-1/2,∞) and α∈(0,2].
  
  
  * Let z,s>0, 0<C≤ 1, and τ∈[C,C^-1]. Then,  there exist c_j=c_j(ζ,C)∈(0,∞), j∈{1,2}, such that
    
    p_ζ^(2)(1,c_1 z,c_1 s)    ≲_C,ζ p_ζ^(2)(τ,z,s) ≲_C,ζ p_ζ^(2)(1,c_2 z,c_2 s), 
    
            p_ζ^(α)(τ,z,s)    ∼_C,ζ,α p_ζ^(α)(1,z,s),   α<2.

  

  
  * Let C,τ>0 and 0<z≤ s/2<∞. Then,  there is c=c(ζ,C)∈(0,∞) with
    
      
      
    p_ζ^(α)(τ,z,s)
                 ≲_C,ζ,α p_ζ^(α)(τ,c,c s)_{τ>C}
          + (τ^-1/2-cs^2/τ/(τ+s^2)^ζ_α=2 + τ/s^2ζ+1+α+τ^2ζ+1+α/α_α∈(0,2))_{τ<C},
            
    
            p_ζ^(α)(τ,z,s)
               ≲_ζ,α s^-(2ζ+1).

    

  
  * Let 0<τ≤1, 0<z≤ s/2, and s≥ C>0. Then, there is c=c(ζ,C)∈(0,∞) with
    
    p_ζ^(α)(τ,z,s) ≲_C,ζ,α p_ζ^(α)(1,c,c s).


  
  * There is c=c(ζ)∈(0,∞) such that for all r,s>0,
    
    min{p_ζ^(α)(1,1,r),p_ζ^(α)(1,1,s)}≲_ζ,α p_ζ^(α)(1,c r,c s).


  
  * Let r,s,z,t>0 with |z-s|>|r-s|/2. Then,  there is c=c(ζ)∈(0,∞) with
    
    p_ζ^(α)(t/2,z,s) ≲_ζ,α p_ζ^(α)(t,cr,cs).

  


The constants in exponential factors in Lemma <ref> may change from place to place when α=2.



The following lemma describes the concentration of mass of p_ζ^(α) at the origin. It is similar to <cit.>.

  
  Let ζ∈(-1/2,∞) and α∈(0,2]. Then, there is c=c_ζ,α∈(0,1) such that for any R>1,
  
    c < ∫_0^2R dw  w^2ζ  p_ζ^(α)(τ,z,w) ≤ 1,   τ∈(0,1], z<R.




  The upper bound follows from the normalization of p_ζ^(α)(τ,z,w).
  We now prove the lower bound in (<ref>).
  First, consider α∈(0,2). Let x=z ∨τ^1/α. Then, for w ∈ (x,x+τ^1/α) ⊂ (0,2R), we have τ^1/α + |z-w| ∼τ^1/α  and τ^1/α+z+w ∼ w ∼ x. Hence, by (<ref>)  and (<ref>), we get
   
    ∫_0^2R dw  w^2ζ  p_ζ^(α)(τ,z,w) ≳∫_x^x+τ^1/α dw τ/(τ^1/α + |z-w|)^1+α (τ^1/α+z+w)^-2ζ w^2ζ
         ≳∫_x^x+τ^1/α dw τ^-1/α = 1.

  Applying (<ref>), by the same proof, we get the result for α=2.
  



 §.§ A 3G inequality for p_ζ^(α) with α∈(0,2)


For the proof of the continuity of the heat kernel for all α∈(0,2) and η∈(-α,(2ζ+1-α)/2] in Theorem <ref>, we will use the following 3G inequality for p_ζ^(α)(t,r,s). 

  
  Let ζ∈(-1/2,∞) and α∈(0,2). For t,τ,t,s,r,z>0, 

    p_ζ^(α)(t,r,z) p_ζ^(α)(τ,z,s)/p^(α)(t+τ,r,s)≲p_ζ^(α)(t,r,z)/(τ^1/α+z+s)^2ζ + p_ζ^(α)(τ,z,s)/(t^1/α+z+r)^2ζ.


It suffices to use (<ref>) and the 3G inequality in <cit.>

    p^(α)(t,r,z) p^(α)(τ,z,s) ≲ p^(α)(t+τ,r,s) [p^(α)(t,r,z) + p^(α)(τ,z,s) ].




There is a similar 3G inequality in <cit.>. However, the present 3G inequality in Lemma <ref> and that in <cit.> do not imply each other.





§ INTEGRAL ANALYSIS OF P_Ζ,Η^(Α)



The integral analysis in this section is a potential-theoretic study based on suitable supermedian functions. In particular, we show that the function s^-η is invariant for the transition density p^(α)_ζ,η,
see (<ref>).



The main results of this section are Theorems <ref> and <ref>, which are key elements in the proof of Theorem <ref>. Let us now sketch the main ideas behind the integral analysis presented in this section. The method, developed in <cit.> and further improved in <cit.>, builds on the results of <cit.>. 

For flexibility in what follows, we introduce an extra parameter β acting as a proxy for η. To
prove Theorem <ref>, we essentially show that Ψ_ζ(β)(_ζ)^-α/2 s^-α-β = s^-β, see (<ref>). 
As a consequence, we obtain (<ref>) for η = 0; see Lemma <ref>. Next, by applying Duhamel's formula and Lemma <ref>, we derive (<ref>), which is almost (<ref>). To fully establish (<ref>), we prove the convergence of the integrals in (<ref>); see Corollary <ref> and Lemma <ref>.

We point out Proposition <ref> as a crucial step in our integral analysis, in which we estimate the mass of p_ζ,η^(α). It is the key to proving the estimates stated in Corollary <ref>; the case η<0 being simpler since the convergence of the integrals follows from the inequality p_ζ,η^(α)≤ p_ζ^(α).

The result stated in Theorems <ref> and <ref> may be viewed as follows. Denote _ζ,η^(α) := _ζ^α/2 - Ψ_ζ(η)s^-α. If we let t→∞ in (<ref>), then we formally get 

    (Ψ_ζ(β) -Ψ_ζ(η)) (_ζ,η^(α))^-1 s^-β-α=  s^-β,
 
which is another expression of 
_ζ,η^(α) s^-β= (Ψ_ζ(β) -Ψ_ζ(η))  s^-β-α. The latter right-hand side changes sign at β=η, so we may informally conclude that for the operator _ζ,η^(α), the function s^-β is subharmonic, harmonic and superharmonic if β <η, β=η and β>η, respectively. In fact, in Theorems <ref> and <ref>, we prove that s^-β is supermedian for p_ζ,η^(α) if β≥η and h(s) = s^-η is even invariant for p_ζ,η^(α), which is crucial to prove the lower bounds. 
As we shall see, the estimates of Bessel and subordinated Bessel heat kernels also play an important role in our development.

When proving ground state representations for (<ref>) in <cit.>, we arrived at the invariant function h(r) by integrating p_ζ^(α)(t,r,s) against suitable functions in space and time, as suggested by the road map in <cit.>.
In the following, we call two parameters β,γ∈ admissible, whenever they satisfy γ∈(-1,∞) and β∈(1,(2ζ-γ)/α).
For admissible β,γ, we recall from <cit.> that the function[The function h_β,γ is also well-defined for β∈(0,(2ζ-γ)/α].]




    h_β,γ(r)
            := ∫_0^∞dt/t  t^β∫_0^∞ ds  s^γ p_ζ^(α)(t,r,s)
          = C^(α)(β,γ,ζ) r^αβ+γ-2ζ,    r>0,

is finite where

    C^(α)(β,γ,ζ)
          := Γ(β)/Γ(αβ/2) C(αβ/2,γ,ζ),   with
    
      C(β,γ,ζ)
          := 2^-2βΓ (β) Γ(1/2 (γ+1)) Γ(1/2 (2ζ-2β-γ))/Γ(1/2(2ζ-γ)) Γ(1/2(2β+γ+1)).

Moreover, the Fitzsimmons ratio (see <cit.>, <cit.>) is

    q_β,γ(r) := (β-1) h_β-1,γ(r)/h_β,γ(r),    r>0.

In particular, we recover

    q(z) = (β-1)h_β-1,γ(z)/h_β,γ(z),    z>0

from (<ref>), whenever β,γ are admissible, to wit, γ∈(-1,∞) and β∈(1,(2ζ-γ)/α), and are chosen so that αβ+γ-2ζ=-η. This is, e.g., the case when β=(2ζ-γ-η)/α with ζ∈(-1/2,∞), α∈(0,2]∩(0,2ζ+1), η∈[0,2ζ+1-α/2], and γ∈(-1,ζ-(α+1)/2). 
We will use the function h_β,γ and its variant h_β,γ^(+) in Section <ref> to obtain invariant functions for p_ζ,η^(α). When doing so, the additional flexibility provided by β and γ will be helpful.  




 §.§ Continuity


For α∈(0,2), we will use the following lemma to prove that p_ζ,η^(α) is a strongly continuous contraction semigroup and that p_ζ,η^(α)(t,r,s) is continuous in t,r,s>0.



 
 Let ζ∈(-1/2,∞), α∈(0,2] and δ∈ (0,2ζ+1). Then,
  
    ∫_0^∞ ds  s^2ζ p_ζ^(α)(τ,r,s) s^-δ∼ r^-δ_{τ<r^α} + τ^-δ/α_{τ>r^α} =r^-δτ^-δ/α

  and
  
    ∫_0^t dτ ∫_0^∞ ds  s^2ζ p_ζ^(α)(τ,r,s) s^-δ ∼ 
              t (r^α∨ t)^-δ/α     for δ<α
    log(1+t/r^α)      for δ=α
    
              r^-δ  (t∧ r^α)       for δ>α.




  We only prove this lemma for α∈(0,2); the proof for α=2 is analogous and, in fact, easier using the explicit expression for p_ζ^(2)(t,r,s).
  Formula (<ref>) follows from integrating (<ref>). Hence, we only need to prove (<ref>).  By scaling we may and do assume r=1. We first show that for τ>1,
  
    ∫_0^∞ ds  s^2ζ  p_ζ^(α)(τ,1,s)s^-δ∼τ^-δ/α.

  Indeed, since τ>1, we have |s-1|< s+1 <3τ^1/α for s<2τ^1/α. Hence, by (<ref>),
  
    ∫_0^2τ^1/α ds  s^2ζ  p_ζ^(α)(τ,1,s) s^-δ∼τ^-(2ζ+1)/α∫_0^2τ^1/α ds    s^2ζ-δ∼τ^-δ/α.

  For s>2τ^1/α, we have s+1 ∼ |s-1| ∼ s. Thus,
  
    ∫_2τ^1/α^∞ ds  s^2ζ p_ζ^(α)(τ,1,s) s^-δ∼τ∫_2τ^1/α^∞ ds   s^-δ-1-α∼τ^-δ/α,

  which yields (<ref>).
  Now, suppose τ≤ 1. Note that s+1+τ^1/α∼ s ∨ 1. Additionally, for s<1/2, p^(α)(τ,1,s) ∼τ. Hence, by (<ref>),
  
    ∫_0^∞ ds  s^2ζ  p_ζ^(α)(τ,1,s) s^-δ≲∫_1/2^∞ ds  p^(α)(τ,1,s) +τ∫_0^1/2 ds  s^2ζ-δ≲ 1,

  where we used that p^(α)(τ,1,s) is the density of a probability distribution.
  Finally, for 1<s<1+τ^1/α, p^(α)(τ,1,s) ∼τ^-1/α. Thus,
  
    ∫_0^∞ ds  s^2ζ  p_ζ^(α)(τ,1,s) s^-δ≳∫_1^1+τ^1/α ds τ^-1/α = 1,

  which, together with (<ref>), yields (<ref>).



We now prove that t↦ p_ζ,η^(α)(t,·,·) is a strongly continuous contraction semigroup. This property will also be important to prove the pointwise continuity of p_ζ,η^(α)(t,r,s) with respect to all r,s,t>0.


  
  Let ζ∈(-1/2,∞), α∈(0,2], and η∈(-α,2ζ+1-α/2].
  Then {p_ζ,η^(α)(t,·,·)}_t≥0 is a strongly continuous contraction semigroup on L^2(_+,r^2ζdr).



  For α=2, Proposition <ref> follows, e.g., from the explicit expression for p_ζ,η^(2) in (<ref>).
  Thus, we consider α∈(0,2) in the following.
  Since we verified the semigroup property in Proposition <ref>, it suffices to verify the strong continuity and the contraction property.
  By the symmetry of p_ζ,η^(α), the supermedian property (<ref>), and a Schur test, p_ζ,η^(α) is a contraction on L^2(_+,r^2ζdr) for every t>0.
  It remains to prove the strong continuity of p_ζ,η^(α)(t,·,·). Since C_c^∞(_+) is dense in L^2(_+,r^2ζdr) and p_ζ^(α)(t,·,·) is strongly continuous (see Lemma <ref>), it follows from Duhamel's formula, that p_ζ,η^(α)(t,·,·) is strongly continuous, if lim_t↘0T_tϕ_L^2(_+,r^2ζdr)=0 holds for every nonnegative function ϕ∈ C_c^∞(_+), where T_t is the integral operator acting as
  
    (T_tϕ)(r)
        = ∫_0^t dτ∫_0^∞ dz  z^2ζ p_ζ^(α)(τ,r,z)q(z) ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t-τ,z,s)ϕ(s).

  By the definition of p_ζ,η^(α), the integral kernel of T_t is maximal for η=(2ζ+1-α)/2 (this is clear for η≥ 0 by (<ref>); for η<0, see Lemma <ref>). Thus, it suffices to consider this case. Since ϕ∈ C_c^∞(_+), we have ϕ≲_ϕ h^*, with h^*(r)=r^-(2ζ+1-α)/2. By the supermedian property
  (<ref>) and Lemma <ref> with δ =  ζ+(1+α)/2 ∈ (0,2ζ+1),
  
    T_tϕ_L^2(_+,r^2ζdr)^2
          ≲_ϕT_t h^*_L^2(_+,r^2ζdr)^2 
         ≲∫_0^∞ dr  r^2ζ(∫_0^t dτ∫_0^∞ dz  z^2ζ p_ζ^(α)(τ,r,z)q(z)h^*(z))^2 
         ≲∫_0^∞ dr  r^-(α+1)(t∧ r^α)^2.

  Here we used αβ+γ-2ζ=-ζ-(1+α)/2<-α, which allows us to use the third line in (<ref>).
  The conclusion follows from an application of the dominated convergence theorem.




 §.§ The case η∈(0,(2ζ+1-α)/2]


One of the main ingredients in the proof of heat kernel estimates is the invariance of the ground state expressed in terms of the heat kernel. 
This is the content of the following important theorem.


  
  Let ζ∈(-1/2,∞), α∈(0,2]∩(0,2ζ+1), η∈(0,(2ζ+1-α)/2], and r,t>0. Then,
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t,r,s)s^-η
        = r^-η.

  Moreover, for β∈[0,2ζ+1-α-η), we have
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t,r,s)s^-β ds 
          = r^-β + (1-Ψ_ζ(β)/Ψ_ζ(η))∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(τ,r,s)  q(s)   s^-β.



For α=2, Theorem <ref> follows, e.g., from the explicit expression for p_ζ,η^(2) in (<ref>). Thus, we may freely assume α∈(0,2) in the following, although all lemmas, propositions, and theorems in this subsection also hold true for α=2.

To prove Theorem <ref>, we need several auxiliary statements, which follow almost exactly as in <cit.>.
In Proposition <ref>, we estimate the mass of p_ζ,η^(α) by 1+(r/t^1/α)^-η.
Auxiliary bounds to that end are contained in Lemmas <ref> and <ref>. Lemmas <ref> and <ref> contain the heart of the proof of Theorem <ref>. Here, we integrate the zeroth and n-th summands of the perturbation series (<ref>) against monomials. We note the role of the auxiliary parameter β, in particular (<ref>), in proving (<ref>). 
We also note that (<ref>) plays an important role to prove the upper bound in (<ref>) in Theorem <ref> for small distances to the origin; see Lemmas <ref> and <ref>.
We apply Lemma <ref> in Corollary <ref> to obtain a weaker version of Theorem <ref>, which will be useful to prove the upper bounds for p_ζ,η^(α); see, in particular, Lemma <ref>.
These results together with the finiteness of certain space and space-time integrals over p_ζ,η^(α) in Lemma <ref> eventually enable us to conclude the proof of Theorem <ref>.
In Corollary <ref>, we collect some convenient consequences of Theorem <ref> that will be used to prove the bounds for p_ζ,η^(α) when η>0.
We finish this section with Lemma <ref>, where we show the finiteness of p_ζ,η^(α)(t,r,s)>0 for all r,s,t>0.

Let us now give the details and start with the following auxiliary bound.

  
  For ζ∈(-1/2,∞), α∈(0,2], η∈(0,(2ζ+1-α)/2], and r,t,R>0,
  
    ∫_0^1 dτ∫_0^R dz  z^2ζ p_ζ,η^(α)(τ,r,z)q(z)
        ≤(2R)^η/c r^-η,

  where c=c_R is the constant appearing in Lemma <ref>.



  By Lemma <ref>, there is c=c_R∈(0,1) such that
  
    c ≤∫_0^2R p_ζ^(α)(1-τ,z,s)s^2ζ ds ≤ 1,
           0<z<R, τ∈(0,1].

  Thus, by Duhamel (<ref>) and the supermedian property (<ref>), we have
  
    ∫_0^1 dτ∫_0^R dz  z^2ζ p_ζ,η^(α)(τ,r,z)q(z) 
         ≤1/c∫_0^2R ds  s^2ζ∫_0^1 dτ∫_0^R dz  z^2ζ p_ζ,η^(α)(τ,r,z)q(z)p_ζ^(α)(1-τ,z,s) 
         ≤1/c∫_0^2Rds  s^2ζ p_ζ,η^(α)(1,r,s)
          ≤(2R)^η/c∫_0^2R ds  s^2ζ p_ζ,η^(α)(1,r,s) s^-η≤(2R)^η/c r^-η.

  This concludes the proof.


We record the following Chapman–Kolmogorov-type inequality for partial sums of the perturbation series (<ref>).


  
  Let ζ∈(-1/2,∞), α∈(0,2], η∈(0,(2ζ+1-α)/2], and r,s,t>0.
  For n∈_0={0,1,2,...}, let
  
    P_t^(n,D)(r,s) := ∑_k=0^n p_t^(k,D)(r,s),

  with p_t^(n,D)(r,s) defined in (<ref>). Then, for all 0<τ<t and n∈_0,
  
    ∫_0^∞ dz  z^2ζ P_τ^(n,D)(r,z)p_ζ^(α)(t-τ,z,s)
        ≤ P_t^(n,D)(r,s).




  By scaling, it suffices to consider 0<τ<1=t.
  For n=0, the statement is just the Chapman–Kolmogorov equality, so let n∈.
  Using the definition of P_t^(n,D) and p_t^(n,D), we get
  
    ∫_0^∞ dz  z^2ζ P_τ^(n,D)(r,z) p_ζ^(α)(1-τ,z,s) 
          = ∫_0^∞ dz  z^2ζ p_τ^(0,D)(r,z) p_ζ^(α)(1-τ,z,s) 
          + ∑_k=1^n ∫_0^∞ dz  z^2ζ∫_0^τ dρ∫_0^∞ dw  w^2ζ p_ρ^(k-1,D)(r,w) q(w)p_ζ^(α)(τ-ρ,w,z) p_ζ^(α)(1-τ,z,s) 
          = p_ζ^(α)(1,r,s)
          + ∑_k=1^n∫_0^τ dρ∫_0^∞ dw  w^2ζ p_ρ^(k-1,D)(r,w) q(w)p_ζ^(α)(1-ρ,w,s) 
         ≤ p_ζ^(α)(1,r,s)
            + ∑_k=1^n∫_0^1 dρ∫_0^∞ dw  w^2ζ p_ρ^(k-1,D)(r,w) q(w)p_ζ^(α)(1-ρ,w,s) 
          = P_1^(n,D)(r,s).

  This concludes the proof.



  
  Let ζ∈(-1/2,∞), α∈(0,2], η∈(0,(2ζ+1-α)/2], and r,s,t>0.
  Then  p_ζ,η^(α)(t,r,s)<∞ for all r>0 and almost all s>0, and there is a constant M≥1 such that
  
    ∫_0^∞ p_ζ,η^(α)(t,r,s)s^2ζ ds
        ≤ M[1+(t^-1/αr)^-η].




  By the scaling (<ref>), it suffices to consider t=1.
  For n∈, we consider the n-th partial sum P_t^(n,D)(r,s) = ∑_k=0^n p_t^(n,D)(r,s), with p_t^(n,D)(r,s) defined in (<ref>). Then, by Duhamel's formula (<ref>) and P_1^(n-1,D)(r,s)≤ P_1^(n,D)(r,s)≤ p_ζ,η^(α)(1,r,s),
  
    P_1^(n,D)(r,s)
             ≤ p_ζ^(α)(1,r,s)
            + ∫_0^1 dτ∫_0^R dz  z^2ζ p_ζ,η^(α)(τ,r,z)(r,z)q(z)p_ζ^(α)(1-τ,z,s) 
          + Ψ_ζ(η)/R^α∫_0^1 dτ∫_R^∞ dz  z^2ζ P_τ^(n,D)(r,z)p_ζ^(α)(1-τ,z,s).

  By Lemma <ref>,
  
    ∫_0^1 dτ∫_R^∞ dz  z^2ζ P_τ^(n,D)(τ,r,z)p_ζ^(α)(1-τ,z,s) ≤ P_1^(n,D)(r,s).

  By induction and the heat kernel bounds in Proposition <ref>, for each n∈_0, we have p_1^(n,D)(r,s)<∞. Thus, for R≥(2Ψ_ζ(η))^1/α∨2 and all r,s>0, we have, by (<ref>),
  
    P_1^(n,D)(r,s) ≤ 2p_ζ^(α)(1,r,s) + 2∫_0^1 dτ∫_0^R dz  z^2ζ p_ζ,η^(α)(τ,r,z)q(z)p_ζ^(α)(1-τ,z,s).

  Integrating both sides of (<ref>) against s^2ζ ds and using ∫_0^∞ p_ζ^(α)(1,r,s)s^2ζ ds=1 for all r>0 by (<ref>), we get, upon letting n→∞ on the left-hand side of (<ref>),
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(1,r,s)
        
        ≲ 1 + ∫_0^1 dτ∫_0^R dz  z^2ζ p_ζ,η^(α)(τ,r,z)q(z)
        ≲ 1 + R^η r^-η.

  Here, we used Lemma <ref> to obtain the final inequality in (<ref>).




  
  
  Let ζ∈(-1/2,∞), α∈(0,2], η∈[0,(2ζ+1-α)/2],
  β∈(0,2ζ+1-α), n∈_0, and r>0.
  Then,
  
    ∫_0^∞ dτ∫_0^∞ ds  s^2ζ p_τ^(n,D)(r,s) s^-β-α
        = Ψ_ζ(η)^n/Ψ_ζ(β)^n+1 r^-β.





  We prove the statement using induction. By (<ref>), we get
  
    ∫_0^∞ dt  ∫_0^∞ ds  s^2ζ  p_ζ^(α)(t,r,s)s^-α-β    = C^(α)(1,2ζ-α-β,ζ) r^-β
        = 2^-αΓ(2ζ+1-α-β/2) Γ(β/2)/Γ(α+β/2) Γ(2ζ+1-β/2) r^-β  = 1/Ψ_ζ(β) r^-β,

  which gives (<ref>) for n=0. Now, suppose (<ref>) holds for some n∈_0. Then, by the definition of p_τ^(n,D)(r,s) in the formula (<ref>), a time-translation τ↦τ+t, (<ref>), and the induction hypothesis,
  
    ∫_0^∞ dτ∫_0^∞ ds  s^2ζ p_τ^(n+1,D)(r,s) s^-β-α
          = ∫_0^∞ dτ∫_0^∞ ds  s^2ζ∫_0^τ dt ∫_0^∞ dz  z^2ζ p_t^(n,D)(r,z)q(z)p_ζ^(α)(τ-t,z,s) s^-β-α
          = ∫_0^∞ dt ∫_0^∞ dz  z^2ζ p_t^(n,D)(r,z)q(z) ∫_0^∞ dτ ∫_0^∞ ds  s^2ζ p_ζ^(α)(τ,z,s)s^-β-α
          = ∫_0^∞ dτ' ∫_0^∞ dz  z^2ζ p_τ'^(n,D)(r,z)Ψ_ζ(η)z^-α·z^-β/Ψ_ζ(β) =Ψ_ζ(η)^n+1/Ψ_ζ(β)^n+2r^-β.

  This concludes the proof.





  
  Let ζ∈(-1/2,∞), α∈(0,2]∩(0,2ζ+1), η∈(0,(2ζ+1-α)/2], and β∈(η,2ζ+1-α-η). Then, for all r>0,
  
    ∫_0^∞ dτ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(τ,r,s)s^-β-α 
        = r^-β/Ψ_ζ(β)-Ψ_ζ(η).




  Since Ψ_ζ(β)>Ψ_ζ(η) (by the symmetry of Ψ_ζ(η) about η=(2ζ+1-α)/2), the claim follows from p_ζ,η^(α)(t,·,·)=∑_n≥0p_t^(n,D), Lemma <ref>, and geometric series.


For η=0, i.e., Ψ_ζ(η)=0, Theorem <ref> is verified by the following lemma.


  
  Let ζ∈(-1/2,∞), α∈(0,2], and
  β∈[0,2ζ+1-α).
  Then, for all r,t>0,
  
    ∫_0^∞ ds  s^2ζ p_ζ^(α)(t,r,s) s^-β
        = r^-β
        - ∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ^(α)(τ,r,s)  Ψ_ζ(β)/s^α  s^-β.




  The assertion for β=0 follows from (<ref>) and Ψ_ζ(0)=0.
  For β>0, by (<ref>) and Chapman-Kolmogorov,
  
    r^-β    =  ∫_0^t dτ ∫_0^∞ ds  s^2ζ  p_ζ^(α)(τ,r,s) Ψ_ζ(β)/s^αs^-β
          + ∫_0^∞ dτ ∫_0^∞ ds  s^2ζ  p_ζ^(α)(t+τ,r,s) Ψ_ζ(β)/s^αs^-β
        = ∫_0^t dτ ∫_0^∞ ds  s^2ζ  p_ζ^(α)(τ,r,s) Ψ_ζ(β)/s^αs^-β
        + ∫_0^∞ dτ ∫_0^∞ ds  s^2ζ∫_0^∞ dz   z^2ζ p_ζ^(α)(t,r,z) p_ζ^(α)(τ,z,s)  Ψ_ζ(β)/s^αs^-β
        = ∫_0^t dτ ∫_0^∞ ds  s^2ζ  p_ζ^(α)(τ,r,s) Ψ_ζ(β)/s^αs^-β + ∫_0^∞ dz  z^2ζ p_ζ^(α)(t,r,z) z^-β,

  where in the last line we used again (<ref>).




  
  Let ζ∈(-1/2,∞), α∈(0,2], η∈(0,(2ζ+1-α)/2), and β∈(0,2ζ+1-α-η). Then, for all r,t>0,
  
    ∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(τ,r,s)s^-β-α < ∞,  
       ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(τ,r,s)s^-β < ∞.




  By (<ref>), for β∈ (η,2ζ+1-α-η), we have 
  
    ∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(τ,r,s)s^-β-α≤1/Ψ_ζ(β)- Ψ_ζ(η) r^-β.

  Now let β∈ (0,η]. Note that β < (2ζ+1-α)/2 and (2ζ+1-α)/2∈ (η,2ζ+1-α-η). By (<ref>) and (<ref>),
  
    ∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(τ,r,s)s^-β-α
       ≤∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(τ,r,s)(1+s^-(2ζ+1-α)/2-α) 
       ≤ M∫_0^t [1+(τ^-1/αr)^-η] dτ + 1/Ψ_ζ((2ζ+1-α)/2)- Ψ_ζ(η) r^-(2ζ+1-α)/2 <∞,

  Hence, (<ref>) follows by (<ref>) and (<ref>). By Lemma <ref>, we have 
  
    ∫_0^∞ ds  s^2ζ p_ζ^(α)(τ,r,s)s^-β≤ r^-β.

  Therefore, by Duhamel's formula, (<ref>), and (<ref>),
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(τ,r,s)s^-β
       ≤∫_0^∞ ds  s^2ζ p_ζ^(α)(τ,r,s)s^-β + ∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(τ,r,s)q(s)s^-β < ∞,

  which ends the proof of (<ref>).
 

We are now ready to give the

  Let t>0 and r>0. By Duhamel's formula,
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t,r,s) s^-β  =  ∫_0^∞ ds s^2ζ p_ζ^(α)(t,r,s) s^-β
                + ∫_0^t dτ∫_0^∞ dz z^2ζp_ζ,η^(α)(τ,r,z) q(z) (∫_0^∞ ds s^2ζp_ζ^(α)(t-τ, z, s) s^- β) .

  Next, by Lemma <ref>,
  
    Ψ_ζ(β) ∫_0^t du∫_0^∞ ds  s^2ζ p_ζ,η^(α)(u,r,s) s^-α - β
        = Ψ_ζ(β) ∫_0^t du ∫_0^∞ ds  s^2ζ p_ζ^(α)(u,r,s) s^-α - β
        + Ψ_ζ(β)  ∫_0^t du ∫_0^∞ ds  s^2ζ∫_0^u dτ∫_0^∞ dz z^2ζ p_ζ,η^(α)(τ,r,z) q(z) p_ζ^(α)(u-τ, z, s) s^-α - β
        = r^-β - ∫_0^∞ s^2ζ p_ζ^(α)(t,r,s) s^-β
        + ∫_0^t dτ∫_0^∞ dz  z^2ζ   p_ζ,η^(α)(τ,r,z) q(z) (Ψ_ζ(β)∫_0^t-τ du ∫_0^∞ ds  s^2ζ p_ζ^(α)(u, z, s) s^-α - β).

  We add (<ref>) and (<ref>) and apply Lemma <ref> to the terms in parentheses to get
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t,r,s) s^-β + Ψ_ζ(β) ∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(τ,r,s) s^-α - β
        = r^-β + Ψ_ζ(η) ∫_0^t dτ∫_0^∞ s^2ζ p_ζ,η^(α)(τ,r,s) s^-α - β.


  Let η∈ (0,(2ζ+1-α)/2). Then  (<ref>) follows by (<ref>) and Lemma <ref>. Furthermore, let β↗δ. By (<ref>), Lemma <ref>, and the Lebesgue convergence theorem, we get (<ref>).

  Now, consider η = (2ζ+1-α)/2 and let β∈ (0,η). Note that p_ζ,η^(α) is the perturbation of p_ζ,β^(α) by (Ψ_ζ(η) - Ψ_ζ(β))s^-α (see, e.g., <cit.>). Hence, by Duhamel's formula, Fubini's theorem, and (<ref>) applied to p_ζ,β^(α), we get
  
    (Ψ_ζ(η)-Ψ_ζ(β)) ∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(τ,r,s) s^-α - β
        =  ∫_0^∞  dz  z^2ζ[(Ψ_ζ(η)-Ψ_ζ(β)) ∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(τ,r,s) s^-α p_ζ,β^(α)(t-τ,s,z)] z^-β
        = ∫_0^∞  dz  z^2ζ[ p_ζ,η^(α)(t,r,z) - p_ζ,β^(α)(t,r,z)] z^-β
        = ∫_0^∞  dz  z^2ζ p_ζ,η^(α)(t,r,z)  z^-β - r^-β,

  which gives (<ref>). Next, let β < η. Since s^-β < 1+ s^-η, by (<ref>) and (<ref>), 
  
    ∫_0^∞ ds  s^2ζ  p_ζ,η^(α)(t,r,s) s^-β≤∫_0^∞ ds  s^2ζ  p_ζ,η^(α)(t,r,s) (1+ s^-η) <∞.

  By (<ref>),
  
    ∫_0^∞ ds  s^2ζ  p_ζ,η^(α)(t,r,s) s^-β dy ≥ r^-β.

  Hence, letting β↗η, by Lebesgue's convergence theorem we get 
  
    ∫_0^∞ ds  s^2ζ  p_ζ,η^(α)(t,r,s) s^-η dy ≥ r^-η,

  which together with (<ref>) gives (<ref>). This concludes the proof of Theorem <ref>.




Using Theorem <ref> and Proposition <ref>, we now prove an approximate supermedian property for

    H(t,r) := 1 + (t^-1/αr)^-η,    r,t>0

and

    H(r) := H(1,r) = 1+r^-η

with respect to p_ζ,η^(α)(t,·,·).

  
  
  Let ζ∈(-1/2,∞), α∈(0,2], and η∈(0,(2ζ+1-α)/2]. Then,  there is a constant M>0 such that for all r,t>0,
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t,r,s)H(t,s)
        ≤ (M+1) H(t,r).

  Furthermore, for all β∈[0,η], 
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t,r,s) (s/t^1/α)^-β≤ M H(t,r).

  Finally,
  
    ∫_0^1 dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z)
        ≤ M H(r).




  Formula (<ref>) follows from (<ref>) and (<ref>).
  Formula (<ref>) follows from s^-β≤ 1+s^-η=H(s), Formula (<ref>) and scaling. To prove Formula (<ref>), we use (<ref>) and Duhamel's formula to obtain
  
    MH(r)
             ≥∫_0^∞ ds  s^2ζ p_ζ,η^(α)(1,r,s) 
        = ∫_0^∞ ds  s^2ζ(p_ζ^(α)(1,r,s) + ∫_0^1 dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) p_ζ^(α)(τ,z,s) ) 
        = 1 + ∫_0^1 dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z).

  This concludes the proof.


Finally, we show the finiteness of p_ζ,η^(α)(t,r,s) for all r,s,t>0.


  
  Let ζ∈(-1/2,∞), α∈(0,2], η∈(0,(2ζ+1-α)/2]. Then p_ζ,η^(α)(t,r,s)<∞ for all r,s,t>0.  



  By the scaling (<ref>), it suffices to consider t=1.
  By (<ref>), we know p_ζ,η^(α)(1,r,s)<∞ for all r>0 and almost all s>0. We now show finiteness for a given r>0 and all s>0. By (<ref>), i.e.,
  
    P_1^(n,D)(1,r,s) ≤ 2p_ζ^(α)(1,r,s) + 2∫_0^1 dτ∫_0^R dz  z^2ζ p_ζ,η^(α)(τ,r,z)q(z)p_ζ^(α)(1-τ,z,s),

  with R=2∧(2Ψ_ζ(η))^1/α and P_t^(n,D)(r,s) = ∑_k=0^n p_t^(n,D)(r,s), it suffices to show
  
    ∫_0^1 dτ∫_0^R dz  z^2ζ p_ζ,η^(α)(τ,r,z)q(z)p_ζ^(α)(1-τ,z,s) < ∞.

  If s≥2R, we use the uniform boundedness of p_ζ^(α)(1-τ,z,s) and (<ref>), which allows to bound the remaining integral by
  
    ∫_0^1dτ∫_0^∞ ds  z^2ζ p_ζ,η^(α)(τ,r,z)q(z)≲ H(r)<∞.

  We assume s<2R from now on.
  We first consider the integral over τ∈(0,1-1/2Ψ_ζ(η)(s/2)^α). Then  the left-hand side of (<ref>) is again finite by the uniform boundedness of p_ζ^(α)(1-τ,z,s) and by (<ref>).
  It remains to consider the integral over τ∈(1-1/2Ψ_ζ(η)(s/2)^α,1). We now distinguish between |z-s|≶ s/2. For |z-s|>s/2, we bound
  
    p_ζ^(α)(1-τ,z,s) ≲1/s^2ζ+1+α

  and use (<ref>). If |z-s|≤ s/2, then z≥ s/2. Thus, by (<ref>) and the Chapman–Kolmogorov equation, we can bound the integral in question by
  
    Ψ_ζ(η)(2/s)^α∫_1-1/2Ψ_ζ(η)(s/2)^α^1 dτ∫_s/2^R dz  p_ζ,η^(α)(τ,r,z)p_ζ,η^(α)(1-τ,z,s) 
         ≤Ψ_ζ(η)(2/s)^α·1/2Ψ_ζ(η)(s/2)^α· p_ζ,η^(α)(1,r,s)
          = 1/2 p_ζ,η^(α)(1,r,s),

  which can be absorbed by the left-hand side of (<ref>) upon taking n→∞. This concludes the proof.




 §.§ The case η<0



In this section, we extend Theorem <ref> to η < 0. Given the explicit expression for p_ζ,η^(2), this is straightforward for α = 2. Hence, the main goal of this section is to generalize the previous integral analysis for α∈ (0, 2) to those coupling constants Ψ_ζ(η) with η∈ (-α, 0). Our approach relies on compensating the kernels p_ζ^(α)(t,r,s) to add extra integrability for large t; see (<ref>).


 


We begin the integral analysis by constructing generalized ground states for α=2 and then use subordination to study α∈(0,2) in the spirit of <cit.>. 
We proceed similarly as in <cit.> by integrating p_ζ^(α) against suitable functions of space and time. This gave us the functions h_β,γ in (<ref>) discussed at the beginning of Section <ref> above. However, since we wish to consider a wider range of parameters, we modify our construction slightly by compensating the kernel p_ζ^(α)(t,r,s) with p_ζ^(α)(t,r,0) when we integrate against monomials growing too fast for large t. The following calculation is similar but more delicate than <cit.>. 


  
  Let ζ∈(-1/2,∞), α∈(0,2], γ∈(-1,∞), and 0<β<(2ζ-γ+2)/α with β≠ (2ζ-γ)/α.
  Then,
  
    h_β,γ^(+)(s)
              := ∫_0^∞dt/t  t^β∫_0^∞ dr  r^γ( p_ζ^(α)(t,r,0) - p_ζ^(α)(t,r,s) ) 
        = C^(α)(β,γ,ζ) s^αβ+γ-2ζ,

  with C^(α)(β,γ,ζ) as in (<ref>).
  
  If, additionally, β=(2ζ-γ-η)/α for some η∈(-α,0), then
  
    q(s)
        = Ψ_ζ(η)/s^α
        = (β-1) h_β-1,γ(s)/h_β,γ^(+)(s)

  where h_β-1,γ(s) is as in (<ref>).


Note that the case αβ=2ζ-γ, which is excluded in the above lemma, corresponds to η=0, as can be seen by considering the right-hand sides of (<ref>) and (<ref>).
The constant on the right-hand side of (<ref>) agrees with that in (<ref>). The subtraction of p_ζ^(α)(t,r,0) in (<ref>), however, allows us to extend the range of admissible β, compare with (<ref>).



  Formula (<ref>) immediately follows from (<ref>), which we prove now. Consider first α=2. For β∈(0,(2ζ-γ)/2), the calculations were already carried out in <cit.>. Thus, consider β∈((2ζ-γ)/2,(2ζ-γ+2)/2) from now on. It is for this range that we require the compensation by p_ζ^(2)(t,r,0) to make the t-integral convergent at t=∞; at t=0, no compensation is necessary in view of the exponential factor -c(r-s)^2/t in p_ζ^(2)(t,r,s) in (<ref>) and similarly for p_ζ^(2)(t,r,0) in (<ref>). We now compute  h_β,γ^(+)(s). By scaling,
  h_β,γ^(+)(s)· s^-γ-2β+2ζ = C^(+)(β,γ,ζ)
  with
  
    C^(+)(β,γ,ζ)
              = ∫_0^∞dt/t  t^β∫_0^∞ dr  r^γ(p_ζ^(2)(t,r,0)-p_ζ^(2)(t,r,1)).

  Since γ>-1, the r-integral gives, as in <cit.>,
  
    ∫_0^∞ dr  r^γ(p_ζ^(2)(t,r,0)-p_ζ^(2)(t,r,1)) 
          = -2^γ -2 ζΓ(γ +1/2) t^1/2 (γ -2 ζ )(  _1F_1(ζ -γ/2;ζ +1/2;-1/4 t)-1)/Γ(ζ +1/2).

  Here, _1F_1(a;b;z) with a,b∈, b∉{0,-1,-2}, and z∈{w∈: |w|<1} denotes Kummer's confluent hypergeometric function <cit.>._1F_1(a;b;z)
  The right-hand side of (<ref>) behaves like t^-ζ+γ/2-1 as t→∞ by <cit.>. To study the behavior at t=0, we use <cit.>, i.e.,
  
    _1F_1(a;b;z) = z  _1F_1(b-a;b;-z)

  and
  
    lim_z→∞ _1F_1(a;b;z) ·1/zz^a-b 
        = 1/Γ(a)Γ(b)

  to infer the existence of some c_γ,ζ∈ such that
  
    _1F_1(ζ -γ/2;ζ +1/2;-1/4 t) =
       c_γ,ζ t^-γ/2+ζ + 𝒪(t^-γ/2+ζ+1)      at t=0.

  Thus, the right-hand side behaves like 1+t^-ζ+γ/2 as t→0.
  These asymptotics and similar computations as in <cit.> thus give
  
    C^(+)(β,γ,ζ)
        = -2^-2 β -1Γ (β) Γ(γ +1/2) Γ(-β +ζ -γ/2)/Γ(β +γ/2+1/2) Γ(-γ/2+ζ)

  for β∈(ζ-γ/2,ζ-γ/2+1).
  This proves (<ref>) for α=2.

  Now consider α∈(0,2). By subordination and the previous computations,
  
    h_β,γ^(+)(s)
              = ∫_0^∞ dr  r^γ∫_0^∞ dτ[p_ζ^(2)(τ,r,s)-p_ζ^(2)(τ,r,0)] ∫_0^∞dt/t σ_t^(α/2)(τ) · t^β
        = Γ(β)/Γ(αβ/2)∫_0^∞ dr  r^γ∫_0^∞ dτ [p_ζ^(2)(τ,r,s)-p_ζ^(2)(τ,r,0)] ·τ^-1+αβ/2
        = Γ(β)/Γ(αβ/2) C(αβ/2,γ,ζ) s^αβ+γ-2ζ
            = C^(α)(β,γ,ζ) s^αβ+γ-2ζ,

  by the definition of C^(α)(β,γ,ζ) in (<ref>).


Here and in the following, we could also work with h_β,γ^(+) defined by (<ref>) when β∈(0,(2ζ-γ)/α). Using the new ansatz (<ref>) involving the difference of heat kernels allows us to increase the range of admissible β and, importantly, work under the less restrictive assumption β<(2ζ-γ+2)/α.


We are now ready to start the actual integral analysis. We begin with a variant of Lemma <ref>.


  
  Let ζ∈(-1/2,∞), α∈(0,2], β∈(-α,0), and r,t>0. Then,
  
    ∫_0^∞ ds  s^2ζ p_ζ^(α)(t,r,s) s^-β
          = r^-β - ∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ^(α)(τ,r,s) Ψ_ζ(β)/s^α s^-β.




  Let β,γ be such that h_β,γ^(+)(r) = C^(α)(β,γ,ζ) r^-β = h_β,γ with h_β,γ^(+)(r) as in (<ref>) and h_β,γ(r) as in (<ref>).
  We start with the Ψ_ζ(β) term. By the expression (<ref>) for the Hardy potential, and (<ref>), we have
  
    Ψ_ζ(β)/s^α· h_β,γ(s)
          = (β-1) ∫_0^∞ dϑ∫_0^∞ dz ϑ^β-2 z^γ p_ζ^(α)(ϑ,s,z).

  Integrating (<ref>) against s^2ζp_ζ^(α)(τ,r,s)_[0,t](τ) ds dτ yields
  
    ∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ^(α)(τ,r,s) Ψ_ζ(β)/s^α· h_β,γ(s) 
          = (β-1) ∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ^(α)(τ,r,s) ∫_0^∞ dϑ∫_0^∞ dz ϑ^β-2 z^γ p_ζ^(α)(ϑ,s,z) 
          = ∫_0^t dτ∫_0^∞ dϑ∫_0^∞ dz  p_ζ^(α)(τ+ϑ,r,z) ∂_ϑϑ^β-1 z^γ
          = -∫_0^t dτ∫_0^∞ dϑ∫_0^∞ dz  [∂_τ p_ζ^(α)(τ+ϑ,r,z)] ϑ^β-1 z^γ
          = ∫_0^∞ dϑ∫_0^∞ dz [p_ζ^(α)(ϑ,r,z) - p_ζ^(α)(ϑ+t,r,z)] ϑ^β-1 z^γ,

  where we used the semigroup property and (β-1)ϑ^β-2=(ϑ^β-1)' in the second step, integrated by parts in the third step[The boundary term at ϑ=∞ vanishes due to the decay of order ϑ^-(2ζ+1)/α for the heat kernel (see (<ref>) and (<ref>)), which suppresses ϑ^β-1 at infinity, since β<(2ζ-γ)/α<(2ζ+1)/α for γ>-1. The boundary term at ϑ=0 vanishes since β>1 and the heat kernels are non-singular at the temporal origin.], and integrated with respect to dτ using the fundamental theorem of calculus. Now we add and subtract p_ζ^(α)(ϑ,0,z) in the integral on the right-hand side of (<ref>), and use (<ref>) to obtain
  
    ∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ^(α)(τ,r,s) Ψ_ζ(β)/s^α· h_β,γ(s) 
          = ∫_0^∞ dϑ∫_0^∞ dz [(p_ζ^(α)(ϑ,r,z) - p_ζ^(α)(ϑ,0,z)) + (p_ζ^(α)(ϑ,0,z) - p_ζ^(α)(ϑ+t,r,z))] ϑ^β-1 z^γ
          = h_β,γ(r) - ∫_0^∞ dϑ∫_0^∞ dz ∫_0^∞ dw  w^2ζ p_ζ^(α)(t,r,w) (p_ζ^(α)(ϑ,w,z)-p_ζ^(α)(ϑ,0,z))  ϑ^β-1 z^γ
          = h_β,γ(r) - ∫_0^∞ dw  w^2ζ p_ζ^(α)(t,r,w) h_β,γ^(+)(w).

  In the last two steps, we used the definition (<ref>) of h_β,γ^(+), the normalization (<ref>), and the semigroup property for p_ζ^(α)(t+ϑ,·,·).


We finish off with a complement of Theorem <ref> for repulsive Hardy potentials
Recall that M=∞ if α=2 and M=α if 0<α<2.

  
  Let ζ∈(-1/2,∞), α∈(0,2]∩(0,2ζ+1), η∈(-M,0), β∈(-M,0), and r,t>0. Then
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t,r,s)s^-β
          = r^-β + (Ψ_ζ(η)-Ψ_ζ(β)) ∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(τ,r,s)s^-β-α.

  In particular, for any r,t>0,
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t,r,s)s^-η = r^-η.




  For α=2, the statements follow from explicit calculations using (<ref>). Thus, suppose α∈(0,2).
  By Duhamel's formula (<ref>) and (<ref>),
  
    ∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(τ,r,s) Ψ_ζ(β)/s^α s^-β
          = ∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ^(α)(τ,r,s) Ψ_ζ(β)/s^α  s^-β
            + ∫_0^t dτ∫_0^∞ ds  s^2ζ-β∫_0^τ dϑ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(ϑ,r,z)q(z) p_ζ^(α)(τ-ϑ,z,s) Ψ_ζ(β)/s^α
          = r^-β - ∫_0^∞ ds  s^2ζ p_ζ^(α)(t,r,s)s^-β
            + ∫_0^t dϑ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(ϑ,r,z)q(z)[z^-β - ∫_0^∞ ds  s^2ζ p_ζ^(α)(t-ϑ,z,s)s^-β].

  Here we applied a time-shift τ↦τ+ϑ and (<ref>) to both lines after the first equality. Bringing the z^-β term to the left-hand side and using
  
    ∫_0^∞ ds  s^2ζ∫_0^t dϑ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(ϑ,r,z)q(z) p_ζ^(α)(t-ϑ,z,s) s^-β
          = ∫_0^∞ ds  s^2ζ (p_ζ,η^(α)(t,r,s) - p_ζ^(α)(t,r,s))s^-β,

  which follows from Duhamel's formula, we obtain
  
    (Ψ_ζ(β)-Ψ_ζ(η))∫_0^t dτ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(τ,r,s)s^-β-α
          = r^-β - ∫_0^∞ ds  s^2ζ p_ζ^(α)(t,r,s)s^-β
            - ∫_0^∞ ds  s^2ζ[p_ζ,η^(α)(t,r,s) - p_ζ^(α)(t,r,s)] s^-β
          = r^-β - ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t,r,s) s^-β.

  This proves (<ref>). Finally, (<ref>) follows from (<ref>) by taking β=η.





§ PROOF OF THEOREM <REF>



In this section, we prove the upper and lower bounds for, and continuity of p_ζ,η^(α)(t,r,s). Moreover, we prove that p_ζ,η^(α) blows up for supercritical coupling constants. 
Before we start the proof of Theorem <ref>, we sketch the proof of (<ref>), which
essentially follows from a change of variables. It was carried out, e.g., by Metafune, Negro, and Spina <cit.>. 


  Recall that p_ζ^(2) is the heat kernel of the nonnegative Bessel operator _ζ in (<ref>); see (<ref>). Using the unitary operator U_ζ:L^2(_+,r^2ζdr)→ L^2(_+,dr) defined by
  
    L^2(_+,r^2ζdr) ∋ u ↦ (U_ζ u)(r)=r^ζu(r) ∈ L^2(_+,dr),

  we have
  
    _ζ - Ψ_ζ(η)/r^2    = U_ζ^*(-d^2/dr^2 + (2ζ-1)^2-1-4Ψ_ζ(η)/4r^2)U_ζ
        = U_η^* U_ζ-η^*(-d^2/dr^2 + (2(ζ-η)-1)^2-1/4r^2)U_ζ-ηU_η
        = U_η^* _ζ-η U_η  in L^2(_+,r^2ζdr).

  Here we used (2ϑ-1)^2-1=(2ζ-1)^2-1-4Ψ_ζ(η) for ϑ∈{ζ-η,1-ζ+η} since Ψ_ζ(η)=η(2ζ-1-η) for α=2. Thus, the change of variables in (<ref>) yields
  
    exp(-t(_ζ-Ψ_ζ(η)/r^2))
        = exp(-tU_η^*_ζ-ηU_η)
        = U_η^* exp(-t_ζ-η) U_η.

  Hence, since p_ζ-η^(2) is the heat kernel of _ζ-η, we see that
  
    p_ζ,η^(2)(t,r,s) = (rs)^-η p_ζ-η^(2)(t,r,s)

  is the heat kernel of _ζ-Ψ_ζ(η)/r^2 in L^2(_+,r^2ζdr), where _ζ is defined in (<ref>).


In the remainder of this section, we focus on α∈(0,2). 
To prove Theorem <ref>, we build on the ideas of <cit.> and <cit.> for the cases η∈(-α,0) and η∈(0,(2ζ+1-α)/2], respectively. Let us first discuss positive η. The main tools we use are Duhamel's formula, the Chapman–Kolmogorov equation, and the scaling given in Proposition <ref>. The first step is to get estimates for p_ζ,η^(α)(1,r,s) with arbitrary r>0 and s≳1. In Lemma <ref> below, we show p_ζ,η^(α)(1,r,s) ≲ H(r)p_ζ^(α)(1,r,s) by properly splitting the integrals in Duhamel's formula and applying the estimates (<ref>) and (<ref>).
For 0<r∨ s≲1, we show p_ζ,η^(α)(1,r,s)≲ H(r) H(s) using p_ζ,η^(α)(2,r,s) = ∫__+ dz  z^2ζ p_ζ,η^(α)(1,r,z) p_ζ,η^(α)(1,z,s) (_z<1 + _z>1) (Chapman–Kolmogorov). This requires bounds for p_ζ,η^(α)(1,r,z) for r≲1 and z≶1.
For r≲1≲ z, we use Lemma <ref> below.
For r∨ z≲1, we use Lemmas <ref>–<ref> below. More precisely, in Lemma <ref>, we prove the preliminary estimate p_ζ,η^(α)(1,r,s) ≲ r^-η s^μ-2ζ-1, which we systematically improve in Lemma <ref> and in the crucial Lemma <ref> to get the upper bounds stated in Lemma (<ref>).
To get the lower bounds in Theorem <ref>, we use the upper bounds stated in Lemma <ref>, the Chapman–Kolmogorov equation, and the invariance of the function h(s) = s^-η. 

For η<0, the Duhamel formula is almost useless. Our proofs are generally based on the Chapman–Kolmogorov equation and the method called "self-improving estimates". This method is applied in the key estimates of the mass M(t,r) of p_ζ,η^(α), defined in (<ref>), which are contained in Proposition <ref>. To that end, we proceed as follows. In the first step, we show M(1,r) < δ M(1,3^1/αr) + A r^-η for sufficiently small δ and some A>0. Next, we iterate this inequality to obtain M(1,r) < δ^n M(1,3^n/αr) + A_n r^-η for some bounded sequence A_n, thereby giving the desired estimates for M(1,r). Proposition <ref> yields the upper estimates of p_ζ,η^(α)(1,r,s) for r,s ≲ 1. The case r,s ≳ 1 is quite obvious because of the bound p_ζ,η^(α)≤ p_ζ^(α). In contrary to positive perturbations, the most challenging case is r ≲ 1 ≲ s. The key estimate is given in Lemma <ref>, which together with <ref> allows to apply the method of "self-improving estimates" in the proof of the inequality  (<ref>). The proof of the lower bounds is based on Lemma <ref> and the estimates of the first summand in the perturbation series.




 §.§ The case η∈(0,(2ζ+1-α)/2]




  §.§.§ Upper bound


The proof of the upper bound in (<ref>) will be concluded in Lemma <ref> after a series of preparatory lemmas involving H(r) = H(1,r) = 1+r^-η defined in (<ref>) above.


In the proofs below, we use the scaling symmetry of p_ζ,η^(α)(t,r,s) to reduce the analysis to the case t=1. Our first step to study p_ζ,η^(α)(1,r,s) is to use Duhamel's formula. In order to discuss "small" and "large" distances to the origin, we will, motivated by the uncertainty principle, introduce a cut-off in the time parameter, defined as

    g(s)
      g(s) := s^α/2^α+1Ψ_ζ(η),    s>0.

The prefactor [2^α+1Ψ_ζ(η)]^-1 is chosen such that certain parts of the integrals appearing in Duhamel's formula are absorbed by p_ζ,η^(α). 
In turn, the time cut-off (<ref>) suggests to distinguish distances to the origin that are smaller or larger than 2(2Ψ_ζ(η))^1/α.
By the scaling symmetry of p_ζ,η^(α)(t,r,s), we could have also considered values of t different from one, which would then require to consider different time and spatial cut-offs.


  All theorems, propositions, and lemmas of this subsection continue to hold for α=2. The only necessary change is that the spatial arguments r,s>0 of the heat kernel p_ζ^(2)(t,r,s) for α=2 should be, in bounds, replaced with cr and cs for some c=c(ζ,α,η)∈(0,∞).


In the following lemma, we examine the region r<2(2Ψ_ζ(η))^1/α, s>0.


  
  Let ζ∈(-1/2,∞),
  α∈(0,2),
  η∈(0,(2ζ+1-α)/2], s>0, and 0<r≤ 2(2Ψ_ζ(η))^1/α. Then,  there is M>0 such that
  
    p_ζ,η^(α)(1,r,s)
             ≤ 2 (p_ζ^(α)(1,r,s) + M p_ζ^(α)(1,1,s)) H(r) 
          + 2∫_g(s)∧ 1^1 dτ∫_s/2^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z) q(z) p_ζ^(α)(τ,z,s) 
          + ∫_0^1/2dτ∫_0^s/2 dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) p_ζ^(α)(τ,z,s).




  By Duhamel's formula (<ref>),
  
    p_ζ,η^(α)(1,r,s) = p_ζ^(α)(1,r,s) + I_1 + I_2 + I_3,

  where
  
    I_1     := ∫_0^1∧ g(s) dτ∫_s/2^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) p_ζ^(α)(τ,z,s), 
    
        I_2     := ∫_1∧ g(s)^1 dτ∫_s/2^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) p_ζ^(α)(τ,z,s), 
    
        I_3     := ∫_0^1 dτ∫_0^s/2 dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) p_ζ^(α)(τ,z,s).

  We leave I_2 untouched. To bound I_1, we use q(z)≤ q(s/2), p_ζ^(α)(τ,z,s)≤ p_ζ,η^(α)(τ,z,s), and the semigroup property to obtain
  
    I_1
             ≤ q(s/2) ∫_0^1∧ g(s)dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)  p_ζ,η^(α)(τ,z,s) 
        =  p_ζ,η^(α)(1,r,s) q(s/2) ∫_0^1∧ g(s)dτ≤1/2 p_ζ,η^(α)(1,r,s).

  Since p_ζ,η^(α)(1,r,s)<∞ by Lemma <ref>, this term can be absorbed by the left-hand side of (<ref>). To treat I_3, we split
  
    I_3
              = ∫_0^1/2dτ∫_0^s/2 dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) p_ζ^(α)(τ,z,s) 
          + ∫_1/2^1 dτ∫_0^s/2 dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) p_ζ^(α)(τ,z,s).

  We leave the first summand untouched. To estimate the second summand, we use p_ζ^(α)(τ,z,s)≲ p_ζ^(α)(1,1,s) for 1/2<τ<1 and 0<z<s/2 (by (<ref>) and (<ref>) in Lemma <ref>), (<ref>), Duhamel's formula (<ref>), and Proposition <ref>, and obtain
  
    ∫_1/2^1 dτ∫_0^s/2 dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) p_ζ^(α)(τ,z,s) 
         ≲ p_ζ^(α)(1,1,s) ∫_0^1 dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) 
          = p_ζ^(α)(1,1,s)∫_0^1dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z)∫_0^∞ dy  y^2ζ p_ζ^(α)(τ,z,y) 
          = p_ζ^(α)(1,1,s)∫_0^∞ dy  y^2ζ∫_0^1dτ∫_0^∞ dz  z^2ζ  p_ζ,η^(α)(1-τ,r,z)q(z)p_ζ^(α)(τ,z,y) 
          = p_ζ^(α)(1,1,s) ∫_0^∞ dy  y^2ζ ( p_ζ,η^(α)(1,r,y) - p_ζ^(α)(1,r,y))
          ≤ M p_ζ^(α)(1,1,s) H(r).

  This concludes the proof of (<ref>).


Lemma <ref> will be particularly important to study the region r∨ s≤ 2(2Ψ_ζ(η))^1/α. Although it could also be used to study r≤2(2Ψ_ζ(η))^1/α≤ s, we will use different methods here.


  
  Let ζ∈(-1/2,∞), α∈(0,2), η∈(0,(2ζ+1-α)/2]. Then,  there is C>0 such that for all s≥2(2Ψ_ζ(η))^1/α and r>0,
  
    p_ζ,η^(α)(1,r,s) ≤  C H(r)  p_ζ^(α)(1,r,s).




  We use (<ref>) for R=(2Ψ_ζ(η))^1/α, to get 
  
    p_ζ,η^(α)(1,r,s)
        ≤ 2 p_ζ^(α)(1,r,s)
        + 2 ∫_0^1dτ∫_0^R dz  z^2ζ p_ζ,η^(α)(τ,r,z) q(z) p_ζ^(α)(1-τ,z,s).

  It suffices to estimate the second summand by p_ζ^(α)(1,r,s)H(r).
  For 0<τ<1, 0<z≤ s/2, s≥ 2R>0, we estimate p_ζ^(α)(1-τ,z,s)≲_η p_ζ^(α)(1,1,s) ((<ref>) in Lemma <ref>) and use (<ref>) to bound
  
    ∫_0^1dτ∫_0^R dz  z^2ζ p_ζ,η^(α)(τ,r,z) q(z) p_ζ^(α)(1-τ,z,s)
        ≲ p_ζ^(α)(1,1,s)  H(r).

  Since for r≤2R, we have p_ζ^(α)(1,1,s) ∼ p_ζ^(α)(1,r,s), we obtain (<ref>) for 0≤ r≤ 2R≤ s. 
  On the other hand, by the symmetry of p_ζ,η^(α)(1,r,s) in r and s, (<ref>) and (<ref>) give
  
    p_ζ,η^(α)(1,r,s) ≲ p_ζ^(α)(1,r,s) + min{p_ζ^(α)(1,1,s)  H(r),p_ζ^(α)(1,1,r)  H(s)},    r,s≥ 2R.

  Since H(s)∼1∼ H(r) for all r,s≥ 2R, and p_ζ^(α)(1,r,1)∧ p_ζ^(α)(1,1,s)≲_η p_ζ^(α)(1,r,s) for all r,s>0 ((<ref>) in Lemma <ref>), we get (<ref>) for r,s≥ 2R.
  This ends the proof.


Lemma <ref>, together with (<ref>) in Theorem <ref> and (<ref>) in Corollary <ref>, allows us to understand the remaining region r,s≤ 2(2Ψ_ζ(η))^1/α.


  
  Let ζ∈(-1/2,∞), α∈(0,2), η∈(0,(2ζ+1-α)/2]. Then, for every μ∈(0,η), there is a constant C_μ>0 such that
  
    p_ζ,η^(α)(1,r,s)
        ≤ C_μ( H(r)s^μ-(2ζ+1)∧ H(s)r^μ-(2ζ+1)),
           0<r,s≤2(2Ψ_ζ(η))^1/α.




  By symmetry of p_ζ,η^(α)(1,r,s), we may assume r<s without loss of generality. We use (<ref>) in Lemma <ref>. 
  First note that p_ζ^(α)(1,r,s)∨ p_ζ^(α)(1,1,s)≲1 for r,s≲ 1.
  Since 0<μ≤ 2ζ+1, by (<ref>) with λ=0 and σ = 2ζ+1-μ, we get p_ζ^(α)(τ,z,s) ≲ s^μ-2ζ-1z^-μ. Therefore, by Lemma <ref>,  (<ref>) in Theorem <ref> and (<ref>) in Corollary <ref>, we get  
  
    p_ζ,η^(α)(1,r,s)       ≲ H(r) + s^μ-(2ζ+1)∫_0^1dτ∫_0^∞dz  z^2ζ p_ζ,η^(α)(1-τ,r,z) z^-μ-α
         ≤ H(r) + s^μ-(2ζ+1)1/Ψ_ζ(η)-Ψ_ζ(μ)∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1,r,z)z^-μ
         ≲ H(r) + s^μ-(2ζ+1)·H(r)/Ψ_ζ(η)-Ψ_ζ(μ)≲ H(r)s^μ-(2ζ+1),
 
  where in the last line we used that s≲ 1 and μ<2ζ+1.




Lemma <ref> is not yet precise enough. The following two lemmas, together with Proposition <ref> and (<ref>) in Theorem <ref> will help us sharpen it. The improved estimates for r,s≤ 2(2Ψ_ζ(η))^1/α are contained in Lemma <ref>.


  
  Let ζ∈(-1/2,∞),
  α∈(0,2),
  η∈(0,(2ζ+1-α)/2], and r>0.
  Then, for all ν∈(η,2ζ+1-η) and δ∈(0,η),
  
    ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(t,r,z)z^-ν≲_ν,δ r^-ν[1+(r/t^1/α)^-δ].

  In particular,
  
    ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(t,r,z)z^-ν≲_ν,δ r^-ν[1+r^-δ],    t∈(0,1].



This lemma can be seen as an extension of (<ref>) in Corollary <ref>, which treated all ν∈(0,η].


  By scaling, it suffices to consider t=1.
  We write R=2(2Ψ_ζ(η))^1/α and distinguish between r≤ R and r>R. We first consider r≤ R and divide the integral in question into two parts.
  Let μ=η-δ so that μ∈(0,η) and μ-ν<0.
  By Lemma <ref>, for r,z≤ R, we have
  
    p_ζ,η^(α)(1,r,z)
          ≲ r^-ηz^μ-(2ζ+1)_{r≤ z} + r^μ-(2ζ+1)z^-η_{z≤ r}.

  
  Then, by Lemma <ref> and (<ref>),
  
    ∫_0^R dz  z^2ζ p_ζ,η^(α)(1,r,z) z^-ν
         ≲_R ∫_0^r dz  z^2ζ-ν· r^μ-(2ζ+1)· z^-η + ∫_r^R dz  z^2ζ-ν· r^-η· z^μ-(2ζ+1)
         ≲ r^-ν-δ,    r≤ R.

  For the complementary integral, we use that Proposition <ref> and ν+δ>η imply
  
    ∫_R^∞ dz  z^2ζ p_ζ,η^(α)(1,r,z)z^-ν≤ R^-ν∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1,r,z)
        ≲_R H(r) ≲_R r^-ν-δ,    r≤ R.

  This concludes the proof for r≤ R.

  For r≥ R we split the z-integration at z=r/2. For z>r/2, we get
  
    ∫_r/2^∞ dz  z^2ζ p_ζ,η^(α)(1,r,z)z^-ν≤ r^-ν∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1,r,z)
          ≲_R r^-ν,   r≥ R.

  For z<r/2, we get, using Lemma <ref> and p_ζ^(α)(1,r,z)≲ r^-(2ζ+1+α) (Proposition <ref>),
  
    ∫_0^r/2 dz  z^2ζ p_ζ,η^(α)(1,r,z)z^-ν≲∫_0^r/2  dz  z^2ζ-ν p_ζ^(α)(1,r,z) H(z)
            ≲ r^-ν-α.

  The proof is concluded.


Let us introduce

    H_ν(r)
      H_ν(r) = 1 + r^-ν,    r, ν>0.

Thus, H_η(r)=H(r), with H(r) as in (<ref>).
By Proposition <ref>, (<ref>) in Corollary <ref>, and Lemma <ref>, we conclude that for all ν∈[0,2ζ+1-η) and δ∈(0,η),

    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t,r,s) H_ν(s)
      ≲_δ H_(ν+δ)∨η(r) ,    t∈(0,1] ,    r>0 .


The following result improves and extends this estimate to all t>0.


  
  Let ζ∈(-1/2,∞),
  α∈(0,2),
  η∈(0,(2ζ+1-α)/2]. Let ν∈[0,2ζ+1-η). Then,
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t,r,s) H_ν(s/t^1/α)
        ≲_ν H(r/t^1/α),    t>0,    r>0.

  In particular,
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t,r,s) H_ν(s)
        ≲_ν t^-ν/α H(r/t^1/α),    t∈(0,1],    r>0.




  By scaling (<ref>) of p_ζ,η^(α), it suffices to prove (<ref>) for t=1.
  If ν≤η, then the claim follows from (<ref>) in Corollary <ref>. 

  Thus, suppose ν=η+ξα with ξ∈(0,(2ζ+1-2η)/α). We distinguish between ξ∈(0,1) and ξ≥1.

  If ξ<1, then (<ref>) in Theorem <ref> and (<ref>) in Corollary <ref> imply for η<μ≤η+α,
  
    ∫_0^1 dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(τ,r,z) H_μ(z)
             ≲ H(r) + ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(1,r,s)s^-μ+α
       ≲ H(r),    r>0.

  Let 0<ϵ<(1-ξ)α∧η so that ν+ϵ=η+ξα+ϵ<η+α. By the semigroup property, ∫_0^1 dτ=1, and Formulae (<ref>) and (<ref>), if ν+ϵ∈(η,η+α) then there is C>0 such that
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(1,r,s) H_ν(s) 
          = ∫_0^1 dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(τ,r,z)∫_0^∞ ds  s^2ζ p_ζ,η^(α)(1-τ,z,s) H_ν(s)
         ≲∫_0^1 dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(τ,r,z) H_ν+ϵ(z)
            ≤ CH(r),    r>0,

  as desired. This concludes the case ξ∈(0,1).

  Now let ξ∈[1,(2ζ+1-2η)/α). By Corollary <ref>, if η+α<μ<2ζ+1-η then
  
    ∫_0^1 dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(τ,r,z) H_μ(z)
        ≲ H_μ-α(r),    r>0.

  We fix ϵ∈(0,1) such that δ:=(1-ϵ)α<2ζ+1-ν-η, and δ<η. By the semigroup property, ∫_0^1dτ=1, (<ref>), and (<ref>), for η+α≤μ≤ν, we have, if η+α≤μ≤ν,
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(1,r,s) H_μ(s) 
          = ∫_0^1 dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(τ,r,z) ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(1-τ,z,s) H_μ(s) 
         ≲∫_0^1 dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(τ,r,z) H_μ+δ(z)
          ≲ H_μ-ϵα(r),    r>0.

  We choose n∈ such that ξα∈[ϵ(nα-1),ϵ(nα-1)+α]∩[α,2ζ+1-2η).
  By the semigroup property, (<ref>), and (<ref>), if η+(ξ-nϵ)α+ϵ∈ (η,η+α), then
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(n+1,r,s) H_ν(s) 
          = ∫_0^∞ dz_1  z_1^2ζ ... ∫_0^∞ dz_n  z_n^2ζ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(1,r,z_1)  ...  p_ζ,η^(α)(1,z_n,s) H_ν(s) 
         ≤ C^n∫_0^∞ dz_1  z_1^2ζ p_ζ,η^(α)(1,r,z_1) H_ν-nϵα(z_1) 
          = C^n ∫_0^∞ dz_1 z_1^2ζ p_ζ,η^(α)(1,r,z_1) H_η+(ξ-nϵ)α(z_1) 
         ≤ C^n+1 H(r),
             r>0,

  recalling ν=η+ξα. This ends the proof in view of p_ζ,η^(α)(n+1,r,s) ∼_n  p_ζ,η^(α)(1,r,s).


Lemmas <ref> and <ref> allow us to improve the estimate for r,s≤ 2(2Ψ_ζ(η))^1/α from Lemma <ref>.


  
  Let ζ∈(-1/2,∞), α∈(0,2), η∈(0,(2ζ+1-α)/2]. Let 0<r,s≤2(2Ψ_ζ(η))^1/α. Then, for each δ∈(0,η),
  
    p_ζ,η^(α)(1,r,s)
        ≲_δ(H(r)s^-η-δ∧ H(s) r^-η-δ).



 
  We will use Lemma <ref> and the symmetry of p_ζ,η^(α). The first line on the right-hand side of (<ref>) is clearly bounded by H(r)s^-η-δ. Now, we estimate the second line on the right-hand side of  (<ref>). Note that for τ∈(1/2,1), we have p_ζ^(α)(τ,z,s)≲ 1 (see (<ref>)). Hence, by (<ref>), we have
 
    ∫_1/2^1 dτ∫_s/2^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) p_ζ^(α)(τ,z,s) 
         ≲∫_0^1 dτ∫_s/2^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) 
            ≲ H(r)
            ≲ H(r)s^-η-δ.

Next, let ε∈ (0,αδ). By (<ref>) with σ = η+δ and λ =α-ε and Lemma <ref>, we get 
  
    ∫_0^1/2dτ∫_0^s/2 dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) p_ζ^(α)(τ,z,s) 
            + ∫_g(s)^1/2 dτ∫_s/2^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) p_ζ^(α)(τ,z,s) 
         ≲ s^-η-δ∫_0^1/2dτ τ^(ε-α)/α∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z) z^-(2ζ+1 -η-δ+ε)
         ≲ s^-η-δ∫_0^1/2 dτ τ^(ε-α)/αH(r/(1-τ)^1/α) · (1-τ)^-(2ζ+1 -η-δ+ε)/α
         ≲ H(r) s^-η-δ.

  This concludes the proof. 



Combining the previous results yields the desired upper bounds in Theorem <ref>.

  
  Let ζ∈(-1/2,∞), α∈(0,2), η∈(0,(2ζ+1-α)/2]. Then,
  
    p_ζ,η^(α)(t,r,s) ≲ p_ζ^(α)(t,r,s) H(t,r)H(t,s),    r,s,t>0.




  By scaling (Formula (<ref>)), it suffices to consider t=1. Let R=2(2Ψ_ζ(η))^1/α. Fix δ∈(0,(2ζ+1)/2-η). By the semigroup property for p_ζ,η^(α) and p_ζ^(α), Lemmas <ref> and <ref>, we have for 0<r,s≤ R,
  
    p_ζ,η^(α)(2,r,s)
              = ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1,r,z)  p_ζ,η^(α)(1,z,s) 
       ≲ H(r)H(s) ∫_0^R dz  z^2ζ· z^-2(η+δ)
            + H(r)H(s) ∫_R^∞ dz  z^2ζ p_ζ^(α)(1,r,z) p_ζ^(α)(1,z,s) 
       ≲_R H(r)H(s) (1+ p_ζ^(α)(2,r,s)),    r,s≤ R.

  Since p_ζ^(α)(2,r,s)∼ p_ζ^(α)(1,r,s)∼_R 1 by (<ref>), this concludes the proof of r,s≤ R.
  Now let r>0 and s≥ R. By Lemma <ref>,
  
    p_ζ,η^(α)(1,r,s)
        ≲ p_ζ^(α)(1,r,s)H(r)
        ∼ p_ζ^(α)(1,r,s) H(r) H(s),    s>R.

  Combining (<ref>) and (<ref>), and using the symmetry of p_ζ,η^(α) yields the claim.




  §.§.§ Lower bound


We prove the lower bound in (<ref>). By scaling (Formula (<ref>)), it suffices to consider t=1. Let q(t,r,s) := p_ζ,η^(α)(t,r,s)/(H(r) H(s)) and μ(ds) = H(s)^2s^2ζ ds. By Doob conditioning, q(t,r,s) is an integral kernel of a semigroup in L^2(_+,μ). By (<ref>) in Theorem <ref>, (<ref>), and Corollary <ref>,

    1 ≤∫_0^∞ q(1,r,s) μ(ds) ≤ M+1,    r>0.

By the upper bound for p_ζ,η^(α) in Lemma <ref>,

    q(1,r,s) ≲ p_ζ^(α)(1,r,s),    r,s>0.

This estimate and the bounds (<ref>) show that there is R>2 such that

    ∫_R^∞ q(1,r,s)μ(ds) ≤1/4,    0<r≤1.

Furthermore, there is 0<ρ<1 such that

    ∫_0^ρ q(1,r,s)μ(ds) ≤1/4,    r>0.

Combining the last two estimates shows

    ∫_ρ^R q(1,r,s)μ(ds) ≥1/2,    0<r<1.


We are now ready to prove the lower bound in (<ref>). We consider the regions 0<r<1<s, and r,s≶1 separately and start with the former. We first note that

    p_ζ^(α)(1,z,s)
           ∼1/|z-s|^1+α(z+s)^2ζ + (1+z+s)^2ζ
       ≥1/R^2ζ+1+α + s^2ζ+1+α + (1+s+R)^2ζ
       ∼ p_ζ^(α)(1,1,s),
           z∈(ρ,R),  s>ρ,

which follows from (<ref>). By the semigroup property, (<ref>), and (<ref>), we have for 0<r≤1 and s≥ρ,

    q(2,r,s)
           ≥∫_ρ^R q(1,r,z)q(1,z,s) μ(dz)
        ≥∫_ρ^R q(1,r,z)p_ζ^(α)(1,z,s)/H(ρ)^2μ(dz) 
       ≳ p_ζ^(α)(1,1,s) ∫_ρ^R q(1,r,z)μ(dz) 
       ≳ p_ζ^(α)(1,1,s),
           0<r≤1,  s≥ρ.

Similarly, we get

    q(3,r,s)
      ≳ p_ζ^(α)(1,1,s)
      ≳ p_ζ^(α)(3,r,s),    0<r<1<s,

where the second estimate follows from (<ref>) if s≥2 and from (<ref>) and (<ref>) if s≤2.

Now let r,s≤1. By the semigroup property, we obtain

    q(3,r,s)
           ≥∫_ρ^R q(1,r,z)q(2,z,s)μ(dz)
        ≳∫_ρ^R q(1,r,z) p_ζ^(α)(1,z,1) μ(dz) 
       ≳_ρ,R p_ζ^(α)(1,1,s)∫_ρ^R q(1,r,z) μ(dz)
        ∼ p_ζ^(α)(1,1,s)
        ∼ p_ζ^(α)(1,r,s),    r,s≤1,

where we used (<ref>) in the second estimate for q(2,z,s). The third and the final estimates in (<ref>) follow from (<ref>) and (<ref>) and (<ref>), since r,s,z≲1.

Finally, if r,s≥1, then q(3,r,s)≥ H(1)^2 p_ζ,η^(α)(3,r,s) ≳ p_ζ^(α)(3,r,s). Combining all the cases shows

    q(3,r,s) ≳ p_ζ^(α)(3,r,s),    r,s>0.

By the definition of q(t,r,s) = p_ζ,η^(α)(t,r,s)/(H(r) H(s)), the claimed lower bound in Theorem <ref> is proved.




  §.§.§ Continuity


To show that p_ζ,η^(α)(t,r,s) is jointly continuous, we will use Duhamel's formula and the known joint continuity of p_ζ^(α)(t,r,s).
To that end, we first prepare the following estimate.


  
  Let ζ∈(-1/2,∞), α∈(0,2), η∈[0,(2ζ+1-α)/2], T>0, and 
  
    G_η^(T)(t,r,s)
       G_η^(T)(t,r,τ,s) := ∫_0^T dz  z^2ζ· z^-α-η  p_ζ^(α)(t,r,z)(τ^1/α+s+z )^-2ζ.

Then, 
  
    G_η^(T)(t,r,τ,s)    ≲ r^-η-α s^-2ζ,         ζ≥0,  r,s,t,τ>0, 
    
          G_η^(T)(t,r,τ,s)    ≲ r^-η-α (τ^1/α+T)^-2ζ,        ζ∈(-1/2,0),  r,t,τ>0, s∈(0,T).




  Let ζ≥ 0. By Lemma <ref>,
  
    G_η^(T)(t,r,τ,s) ≤ s^-2ζ∫_0^∞ dz  z^2ζ· z^-α-η  p_ζ^(α)(t,r,z) ≲ s^-2ζ r^-η-α.

  If ζ∈ (-1/2,0), then
  
    G_η^(T)(t,r,τ,s) ≤∫_0^∞ dz  z^2ζ· z^-α-η  p_ζ^(α)(t,r,z) (τ^1/α +2T)^-2ζ≲ (τ^1/α +2T)^-2ζ r^-η-α.

  This concludes the proof.



  
  Let ζ∈ (-1/2,∞), α∈(0,2), and η∈[0,(2ζ+1-α)/2]. Then, for every r,t>0, the function (0,∞)∋ s↦ p_ζ,η^(α)(t,r,s) is continuous.



  The continuity of p_ζ^(α)(t,r,s) is well known, hence let η>0 from now on. By scaling, it suffices to consider t=1. Fix r,s>0 and let w>0. Then,
  
    p_ζ,η^(α)(1,r,s) - p_ζ,η^(α)(1,r,w) 
          = p_ζ^(α)(1,r,s) - p_ζ^(α)(1,r,w) 
            + ∫_0^1 dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z)[p_ζ^(α)(τ,z,s)- p_ζ^(α)(τ,z,w)].

  We claim that the integral vanishes as w converges to s. Indeed, let ϵ∈(0,1/2).
  Then, for τ∈(ϵ,1) and |w-s| sufficiently small, we have p_ζ^(α)(τ,z,w)∼ p_ζ^(α)(τ,z,s) by (<ref>). Thus, by dominated convergence,
  
    lim_w→ s∫_ϵ^1 dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) [p_ζ^(α)(τ,z,s) - p_ζ^(α)(τ,z,w)] = 0.

  Furthermore, by the upper heat kernel bounds in Theorem <ref> and the 3G inequality in Lemma <ref>, together with Lemma <ref>, we have, for any ξ∈(s/2,2s) and ϵ∈(0,1/2),
  
    ∫_0^ϵ dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) p_ζ^(α)(τ,z,ξ) 
         ≲∫_0^ϵ dτ∫_0^∞ dz  z^2ζ p_ζ^(α)(1-τ,r,z) H(1-τ,r)H(1-τ,z)q(z) p_ζ^(α)(τ,z,ξ) 
         ≲ p^(α)(1,r,s) H(r) ∫_0^ϵ dτ[G_η^2(r∨ξ∨1)(1-τ,r,τ,ξ) + G_0^2(r∨ξ∨1)(1-τ,r,τ,ξ) . 
                                 . + G_η^2(r∨ξ∨1)(τ,ξ,1-τ,r) + G_0^2(r∨ξ∨1)(τ,ξ,1-τ,r)] 
            + H(r) ∫_0^ϵ dτ∫_2(r∨ξ∨ 1)^∞ dz  z^-2ζ -3α-2· (1+z^-η) 
         ≲ϵ H(r)p^(α)(1,r,ξ)[(r^-η-α+r^-α)ξ^-2ζ + (ξ^-η-α+ξ^-α)r^-2ζ]_ζ≥0
            + ϵ H(r) p^(α)(1,r,ξ) [r^-η-α+r^-α/(1+r+ξ)^2ζ + ξ^-η-α+ξ^-α/((1-τ)^1/α+1+r+ξ)^2ζ]_ζ∈(-1/2,0)
            + ϵ H(r) (r+ξ+1)^-2ζ-3α-1
         ≲ϵ H(r) c(r,s),

  where c(r,s)=c_ζ,η,α(r,s) depends only on r,s>0 and the fixed parameters ζ,α, and η. Using (<ref>), we see that
  
    lim_ϵ→0∫_0^ϵ dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z)[p_ζ^(α)(τ,z,s)- p_ζ^(α)(τ,z,w)] = 0,

  uniformly in w∈(s/2,2s).
  This yields the result.



  
  Let ζ∈ (-1/2,∞), α∈(0,2), and η∈(0,(2ζ+1-α)/2]. Then, p_ζ,η^(α)(t,r,s) is jointly continuous in t,r,s>0.



  By scaling, it suffices to consider t=1. Fix r,s>0 and let r̃ and s̃ be close to r and s respectively. As in Lemma <ref>, it suffices to show
  
    ∫_0^1 dτ∫_0^∞ dz  z^2ζ| p_ζ,η^(α)(1-τ,r̃,z) p_ζ^(α)(τ,z,s̃) - p_ζ,η^(α)(1-τ,r,z) p_ζ^(α)(τ,z,s)| q(z) → 0

  as r̃→ r and s̃→ s. In addition to (<ref>) we get (with the same arguments)
  
    ∫_1-ϵ^1 dτ∫_0^∞ dz  z^2ζ p_ζ,η^(α)(1-τ,r,z)q(z) p_ζ^(α)(τ,z,s) 
         ≲∫_0^ϵ dτ∫_0^∞ dz  z^2ζ p_ζ^(α)(τ,r,z) H(τ,r)H(τ,z) q(z) p_ζ^(α)(1-τ,z,s) ≤ϵ c(r,s)

  for a constant c(r,s)>0 independent of ε. 
  By (<ref>)–(<ref>), we get
  
    ∫_0^1 dτ∫_0^∞ dz  z^2ζ| p_ζ,η^(α)(1-τ,r̃,z) p_ζ^(α)(τ,z,s̃) -  p_ζ,η^(α)(1-τ,r,z)p_ζ^(α)(τ,z,s)| q(z) 
         ≤∫_0^ϵ dτ∫_0^∞ dz  z^2ζ[ p_ζ,η^(α)(1-τ,r̃,z) p_ζ^(α)(τ,z,s̃) + p_ζ,η^(α)(1-τ,r,z) p_ζ^(α)(τ,z,s)] q(z) 
            + ∫_1-ϵ^1 dτ∫_0^∞ dz  z^2ζ[ p_ζ,η^(α)(1-τ,r̃,z) p_ζ^(α)(τ,z,s̃) + p_ζ,η^(α)(1-τ,r,z) p_ζ^(α)(τ,z,s)] q(z) 
            + ∫_ϵ^1-ϵ dτ∫_0^∞ dz  z^2ζ| p_ζ,η^(α)(1-τ,r̃,z) p_ζ^(α)(τ,z,s̃) - p_ζ,η^(α)(1-τ,r,z) p_ζ^(α)(τ,z,s)| q(z) 
         ≤ϵ c(r,s) + ∫_ϵ^1-ϵ∫_0^∞ dz  z^2ζ| p_ζ,η^(α)(1-τ,r̃,z) p_ζ^(α)(τ,z,s̃) - p_ζ,η^(α)(1-τ,r,z)p_ζ^(α)(τ,z,s) |.

  By the upper and lower heat kernel estimates in Theorem <ref>, we have, for ϵ<τ<1-ϵ, the estimates
  p_ζ,η^(α)(1-τ,r̃,z) ∼ p_ζ,η^(α)(1-τ,r,z) and
  p_ζ^(α)(τ,z,s̃) ∼ p_ζ^(α)(τ,z,s).
  By Lemma <ref> and dominated convergence, the last integral is arbitrarily small. This concludes the proof of Lemma <ref> and in particular that of the continuity statement in Theorem <ref>.




  §.§.§ Blowup



We now prove Corollary <ref> for all α∈(0,2] using Theorem <ref>. Note that the pointwise bounds stated in Theorem <ref> extend to α=2 in view of Remark <ref>.


  Let η_*:=(2ζ+1-α)/2 be the parameter corresponding to the critical coupling constant κ_ c(ζ,α) and recall our assumption κ>κ_ c(ζ,α). 
According to  the Appendix <ref>,
for t,r,s>0, we define the Schrödinger perturbation p̃_κ(t,r,s) 
of p_ζ^(α)(t,r,s) by κ/r^α:
  
    p̃_κ(t,r,s)
              := ∑_n≥0 p_κ^(n)(t,r,s),   where  
            p_κ^(0)(t,r,s) := p_ζ^(α)(t,r,s),   and
    
          p_κ^(n)(t,r,s)     := ∫_0^t dτ∫_0^∞ dz  z^2ζ p_ζ^(α)(t,r,z) κ/z^α p_κ^(n-1)(t-τ,z,s),    n∈.

  Our goal is to prove that p̃_κ(t,r,s)=∞ for s<t.
  By the definitions of p̃_κ and p_ζ,η^(α), we have p̃_κ≥ p_ζ,η_*^(α).
  In fact, according to <cit.>, p̃_κ may be considered as a Schrödinger perturbation of p_ζ,η_*^(α) by the positive potential (κ-κ_ c(ζ,α))r^-α. Using the second term in the resulting perturbation series 
and the lower heat kernel bound in Theorem <ref>, we get
  
    p̃_κ(t,r,s)
           ≥(κ-κ_ c(ζ,α)) ∫_0^t dτ∫_0^∞ dz  z^2ζ p_ζ,η_*^(α)(t-τ,r,z) z^-α p_ζ,η_*^(α)(τ,z,s)
        = ∞,

for s<t.  This concludes the proof.





 §.§ The case η∈(-α,0)


Throughout this section, we assume α∈(0,2).
Similarly as in (<ref>)–(<ref>), and inspired by Proposition <ref> and 
Theorem <ref>, we consider the total mass of the kernel p_ζ,η^(α),

    M(t,r) := ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t,r,s).

By the scaling (<ref>), we have the scaling relation

    M(t,r)
      = ∫_0^∞ ds  s^2ζ t^-(2ζ+1)/α p_ζ,η^(α)(1,t^-1/αr,t^-1/αs)
      = M(1,t^-1/αr).

Note that, by Theorem <ref>, p_ζ,η^(α) satisfies the Chapman–Kolmogorov equations (<ref>), the Duhamel formulae (<ref>), and

    0 < p_ζ,η^(α)(t,r,s) ≤ p_ζ^(α)(t,r,s)   for all r,s,t>0.

In particular, by the normalization (<ref>), we have

    0 ≤ M(t,r) ≤∫_0^∞ ds  s^2ζ p_ζ^(α)(t,r,s) = 1.





  §.§.§ Upper bound



  
  Let ζ∈(-1/2,∞), α∈(0,2), and η∈(-α,0). Then, there is C>0 such that for all r>0,
  
    M(1,r) ≤ C(1∧ r^-η).




  By Chapman–Kolmogorov and p_ζ,η^(α)(1/3,z,w)≤ c for all z,w>0 (by (<ref>)), we have
  
    p_ζ,η^(α)(1,r,s)
              = ∫_0^∞ dz  z^2ζ∫_0^∞ dw  w^2ζ p_ζ,η^(α)(1/3,r,z) p_ζ,η^(α)(1/3,z,w) p_ζ,η^(α)(1/3,w,s) 
       ≤ c ∫_0^∞ dz  z^2ζ∫_0^∞ dw  w^2ζ p_ζ,η^(α)(1/3,r,z) p_ζ,η^(α)(1/3,w,s) 
        = c M(1/3,r) M(1/3,s) = c M(1,3^1/αr) M(1,3^1/αs).

  Next, let ρ>0 be so small that δ:=c· (2ζ+1)^-1ρ^2ζ+1 < 3^η/α. Recall η<0. Then, by (<ref>), (<ref>), and (<ref>), we have for any r>0,
  
    M(1,r)
             ≤∫_0^ρ ds  s^2ζ p_ζ,η^(α)(1,r,s) + ρ^η∫_ρ^∞ ds  s^2ζ p_ζ,η^(α)(1,r,s) s^-η
       ≤∫_0^ρ ds  s^2ζ p_ζ,η^(α)(1,r,s) + A r^-η
       ≤ c∫_0^ρ ds  s^2ζ M(1,3^1/αr) M(1,3^1/αs) + A r^-η
       ≤δ M(1,3^1/αr) + A r^-η,

  where A:=ρ^η. Iterating (<ref>) yields for all r>0,
  
    M(1,r)
             ≤δ M(1,3^1/αr) + A r^-η
       ≤δ[δ M(1,3^2/αr) + A (3^1/αr)^-η] + A r^-η
       ≤δ^2 [δ M(1,3^3/αr) + A (3^2/αr)^-η] + A (1+δ 3^-η/α) r^-η
       ≤⋯
       ≤δ^n M(1,3^n/αr) + A  [1 + δ 3^-η/α + ⋯ + (δ 3^-η/α)^n-1] r^-η.

  Since δ<1 and δ 3^-η/α<1, the above geometric series converges. Hence, by (<ref>), we obtain
  
    M(1,r) ≤A/1-δ 3^-η/α r^-η,

  which concludes the proof.


Applying Proposition <ref> to (<ref>) immediately yields


  
  Let ζ∈(-1/2,∞), α∈(0,2), η∈(-α,0). Then, there is a constant C>0 such that
  
    p_ζ,η^(α)(1,r,s)
        ≤ C (1∧ r^-η) (1∧ s^-η),
           r,s>0.



We now refine these upper bounds for p_ζ,η^(α)(t,r,s). 


  
  Let ζ∈(-1/2,∞), α∈(0,2), and η∈ (-α,0). Then, for any r,s,t>0, we have
  
    ∫_|z-s|<|r-s|/2 dz  z^2ζ p_ζ^(α)(t,r,z) p_ζ^(α)(t,z,s)
        ≤p_ζ^(α)(2t,r,s)/2.




  Fix r,s,t>0. Then, by p_ζ^(α)(t,r,s) = p_ζ^(α)(t,s,r) for all r,s,t>0,
  
    ∫_|z-s|<|r-s|/2 dz  z^2ζ p_ζ^(α)(t,r,z) p_ζ^(α)(t,z,s)
        = ∫_|z-r|<|r-s|/2 dz  z^2ζ p_ζ^(α)(t,r,z) p_ζ^(α)(t,z,s).

  By Chapman–Kolmogorov,
  
    2∫_|z-s|<|r-s|/2 dz  z^2ζ p_ζ^(α)(t,r,z)p_ζ^(α)(t,z,s) 
          = ∫_|z-s|<|r-s|/2 dz  z^2ζ p_ζ^(α)(t,r,z)p_ζ^(α)(t,z,s)
          + ∫_|z-r|<|r-s|/2 dz  z^2ζ p_ζ^(α)(t,r,z)p_ζ^(α)(t,z,s) 
         ≤∫_0^∞ dz  z^2ζ p_ζ^(α)(t,r,z) p_ζ^(α)(t,z,s)
          = p_ζ^(α)(2t,r,s),

  which concludes the proof.




  
  Let ζ∈(-1/2,∞), α∈(0,2), and η∈(-α,0).
  Let also h(t,r):=(r t^-1/α)^-η.
  Then, there is a constant A>0 such that for all r,s,t>0, we have
  
    p_ζ,η^(α)(t,r,s)
        ≤∫_|z-s|<|r-s|/2 dz  z^2ζ p_ζ,η^(α)(t/2,r,z) p_ζ,η^(α)(t/2,z,s) + A h(t,r) p_ζ^(α)(t,r,s).




  By Chapman–Kolmogorov for p_ζ,η^(α), we have for any r,s,t>0,
  
    p_ζ,η^(α)(t,r,s)
              = ∫_|z-s|≤|r-s|/2 dz  z^2ζ p_ζ,η^(α)(t/2,r,z) p_ζ,η^(α)(t/2,z,s) 
          + ∫_|z-s|≥|r-s|/2 dz  z^2ζ p_ζ,η^(α)(t/2,r,z) p_ζ,η^(α)(t/2,z,s).

  By (<ref>) and Proposition <ref>, we have
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(t,r,s)
        = M(t,r) = M(1,t^-1/αr)
        ≤ c  h(t,r),    r,t>0.

  For t>0 and r,s,z>0 with |z-s|>|r-s|/2, by (<ref>) and (<ref>), we have
  
    p_ζ,η^(α)(t/2,z,s)
        ≤ p_ζ^(α)(t/2,z,s)
        ≲ p_ζ^(α)(t,r,s).

  This implies for r,s,t>0,
  
    ∫_|z-s|≥|r-s|/2dz  z^2ζ p_ζ,η^(α)(t/2,r,z) p_ζ,η^(α)(t/2,z,s) 
         ≲ p_ζ^(α)(t,r,s) ∫_|z-s|≥|r-s|/2dz  z^2ζ p_ζ,η^(α)(t/2,r,z) 
         ≲ p_ζ^(α)(t,r,s) · h(t/2,r)  
          = 2^-η/α h(t,r)p_ζ^(α)(t,r,s),

  which concludes the proof.


We are now in position to derive the upper bounds in Theorem <ref> in the following key lemma.


  
  Let ζ∈(-1/2,∞), α∈(0,2), and η∈(-α,0). Then, for all r,s,t>0,
  
    p_ζ,η^(α)(t,r,s)
        ≲_ζ,α,η (1∧ (t^-1/αr)^-η) · (1∧ (t^-1/αs)^-η) · p_ζ^(α)(t,r,s).




  By the scaling relation (<ref>), it suffices to consider t=1. Let δ=1/2 and ν=2^(-η-α)/α<1. As in Lemma <ref>, let h(t,r)=(r/t^1/α)^-η. Note that
  
    δ h(t/2,r)
        = 1/2 ((t/2)^-1/αr)^-η
        = 2^(-η-α)/α (t^-1/αr)^-η
        = ν h(t,r),    r,t>0.
 
  Let A be the constant from Lemma <ref>. We claim that for any n∈_0, we have
  
    p_ζ,η^(α)(t,r,s)
        ≤[δ^n+1 + (1+ν+⋯+ν^n) A h(t,r)] p_ζ^(α)(t,r,s),
           r,s,t>0.

  We prove (<ref>) by induction. For n=0, t>0, and r,s>0, Lemmas <ref> (together with
  p_ζ,η^(α)(t,r,s) ≤ p_ζ^(α)(t,r,s)) and <ref> imply
  
    p_ζ,η^(α)(t,r,s)
        ≤[ δ + A h(t,r) ] p_ζ^(α)(t,r,s).

  We now make the induction step. Thus, suppose (<ref>) was true for fixed n∈, i.e.,
  
    p_ζ,η^(α)(t,r,s)
        ≤[δ^n + (1+ν+⋯+ν^n-1) A h(t,r) ] p_ζ^(α)(t,r,s),
           r,s,t>0.

  Then, for any r,s,t>0, Lemmas <ref> and <ref>, together with (<ref>) yield
  
    p_ζ,η^(α)(t,r,s)
           ≤∫_|z-s|≤|r-s|/2dz  z^2ζ p_ζ,η^(α)( t2,r,z) p_ζ,η^(α)( t2,z,s) + A h(t,r) p_ζ^(α)(t,r,s) 
        ≤∫_|z-s|≤|r-s|/2dz  z^2ζ[δ^n + (1+ν+⋯+ν^n-1) A h( t2,r)] p_ζ^(α)( t2,r,z) p_ζ^(α)( t2,z,s) 
          + A h(t,r) p_ζ^(α)(t,r,s) 
        ≤[ δ^n + (1+ν+⋯+ν^n-1) A h(t/2,r) ] ·δ p_ζ^(α)(t,r,s) + A h(t,r) p_ζ^(α)(t,r,s) 
        ≤[δ^n+1 + (ν+⋯+ν^n) A h(t,r) ] · p_ζ^(α)(t,r,s) + A h(t,r)p_ζ^(α)(t,r,s) 
        = [ δ^n+1 + (1+ν+⋯+ν^n+1) A h(t,r)] · p_ζ^(α)(t,r,s),

  which proves (<ref>). Taking t=1, observing h(1,r)=r^-η and letting n→∞ in (<ref>) thus yields
  
    p_ζ,η^(α)(1,r,s)
        ≤A/1-ν r^-η p_ζ^(α)(1,r,s),
           r,s>0.

  This estimate allows us to conclude the claimed upper bound (<ref>) (with t=1). By the symmetry of p_ζ,η^(α) it suffices to assume 0<r<s.
  For 1<r<s, the claim (<ref>) follows from p_ζ,η^(α)≤ p_ζ^(α)(1,r,s).
  For 0<r<s<1 we use Corollary <ref> as well as p_ζ^(α)(1,r,s)≳1.
  Finally, for 0<r<1<s, the claim (<ref>) follows from (<ref>).
  This concludes the proof of Lemma <ref> and thereby the upper bound in Theorem <ref>.





  §.§.§ Lower bound


We begin with a lower bound on the function M(t,r), defined in (<ref>).


  
  Let ζ∈(-1/2,∞), α∈(0,2), and η∈(-α,0). Then, there is C>0 such that
  
    M(1,r) ≥ C(1∧ r^-η),    r>0.




  Let R>0 and 0<r<R/2. Then, by (<ref>), there are R-independent numbers c_1,c_2>0 such that
  
    ∫_R^∞ ds  s^2ζ p_ζ^(α)(1,r,s) s^-η≤ c_1 ∫_R^∞ ds  s^2ζ·s^-η/s^2ζ+1+α
        = c_2 R^-η-α,

  and the right-hand side vanishes as R→∞. Let C≡ C_ζ,α,η be the constant in the upper heat kernel bounds (<ref>) and choose R≥1 so large that c_2CR^-η-α≤1/2. Then,  by (<ref>) in Theorem <ref>, for ρ≥ R and r∈(0,ρ/2), we have
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(1,r,s)
            ≥ρ^η∫_0^ρ ds  s^2ζ p_ζ,η^(α)(1,r,s) s^-η
          = ρ^η( r^-η - ∫_ρ^∞ ds  s^2ζ p_ζ,η^(α)(1,r,s) s^-η) 
         ≥ρ^η(r^-η - C∫_ρ^∞ ds  s^2ζ· r^-η p_ζ^(α)(1,r,s) s^-η)
          ≥r^-η/2ρ^-η.

  Thus, we proved
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(1,r,s)
          ≥r^-η/2R^-η≥1∧ r^-η/2R^-η,    0<r<R/2.

  On the other hand, if r>R/2, then, by taking ρ=2r+1 in (<ref>), we obtain
  
    ∫_0^∞ ds  s^2ζ p_ζ,η^(α)(1,r,s)
        ≥r^-η/2(2r+1)^-η≥r^-η/2(4r)^-η≥1/4^-η+1/2,    r>R/2.

  Combining (<ref>)–(<ref>) yields (<ref>) and concludes the proof.


To get the desired lower bounds in Theorem <ref>, we compare p_ζ^(α)(t,r,s) with p_ζ,η^(α)(t,r,s) and distinguish between r∨ s≲1 and r∧ s≳1. We first focus on the case r∧ s≳1. 
To that end, recall the function p_t^(1,D)(r,s) from the perturbation series in (<ref>).
By the scaling (<ref>) of p_ζ^(α), we have the scaling

    p_t^(1,D)(r,s) = t^-2ζ+1/α p_1^(1,D)(t^-1/αr,t^-1/αs),
         r,s,t>0.


The statement of the following lemma is the bound in (<ref>), proved in Appendix <ref>. It is closely related to estimates in <cit.> or <cit.>. It is helpful to obtain lower bounds for p_ζ,η^(α)(1,r,s) when r,s≳1 from upper bounds for p_1^(1,D), since Ψ_ζ(η)<0 for η<0.


  
  Let ζ∈(-1/2,∞), α∈(0,2), and η∈(-α,0). Then, for all r,s,t>0, 
  
    p_ζ,η^(α)(t,r,s)
        ≥ p_ζ^(α)(t,r,s) exp(p_t^(1,D)(r,s)/p_ζ^(α)(t,r,s)).



Our task thus consists of finding good upper bounds for p_1^(1,D)(r,s) when r,s≳1. 
The following lemma, in combination with Lemma <ref>, is important for the derivation of lower bounds for p_ζ,η^(α)(1,r,s) when r,s≳1. 


  
  Let ζ∈(-1/2,∞), α∈(0,2). Then, for all r,s>1, we have
  
    p_1^(1,D)(r,s)
        ≲  p_ζ^(α)(1,r,s).



Note that <cit.> proved both upper and lower bounds. However, the upper bounds suffice for the subsequent analysis.


  It suffices, by Chapman–Kolmogorov and (<ref>), to estimate
  
    ∫_0^1 dτ∫_0^∞ dz  z^2ζ p_ζ^(α)(1-τ,r,z) z^-α p_ζ^(α)(τ,z,s) 
         ≤ 2^α p_ζ^(α)(1,r,s) 
          +  ∫_0^1 dτ∫_0^1/2 dz  z^2ζ p_ζ^(α)(1-τ,r,z) z^-α p_ζ^(α)(τ,z,s) 
         ≲  p_ζ^(α)(1,r,s) 
          +  ∫_0^1 dτ∫_0^1/2 dz  z^2ζ-α1/(1+r)^1+α+2ζ(1+s)^1+α+2ζ
         ≲  p_ζ^(α)(1,r,s).

  In the last estimate we used again (<ref>) and the inequality [(1+r)(1+s)]^-1-α-2ζ≲  (1+|r-s|)^-1-α(1+r+s)^-2ζ valid for all ζ>-1/2.



Recall that Ψ_ζ(η)<0 for η<0. Thus, Lemmas <ref> and <ref> lead to the following useful estimate, which suffices to prove the aspired lower heat kernel bounds for p_ζ,η^(α)(1,r,s) when r∧ s≳1.


  
  
  Let ζ∈(-1/2,∞), α∈(0,2), and η∈(-α,0). Then,  
  
    p_ζ,η^(α)(t,r,s)
        ≳ p_ζ^(α)(t,r,s),     r,s≳ t^1/α>0.




Since Ψ_ζ(η)<0, the result follows by Lemmas <ref> and <ref>.




Corollary <ref> suffices for the derivation of the lower heat kernel bounds for r∧ s≳1. The following lemma will be used to derive the lower heat kernel bounds when r∨ s≲1. In combination with Corollary <ref> it will also handle the case r≲ 1 ≲ s.


  
  Let ζ∈(-1/2,∞), α∈(0,2), and η∈(-α,0). 
  Then, there are R>1 and C_R>0 such that for all r,s>0 with r∨ s≤ R,
  
    p_ζ,η^(α)(1,r,s) ≥ C_R (r· s)^-η.




  By (<ref>) there is R_0>2· 3^1/α large enough such that for all 0<r<R_0/2,
  
    ∫_0^R_0 ds  s^2ζ p_ζ,η^(α)(1,r,s)
        ≥r^-η/2R_0^-η.

  On the other hand, the upper heat kernel bounds (<ref>) imply for all ρ_0>0 and r>0,
  
    ∫_0^ρ_0 ds  s^2ζ p_ζ,η^(α)(1,r,s)
        ≤ C r^-η∫_0^ρ_0 ds  s^2ζ p_ζ^(α)(1,r,s) s^-η≤ C r^-η·ρ_0^-η

  where we used (<ref>) and -η>0 in the last step. We now take ρ_0=(4C)^1/ηR_0^-1. Then, for r<R_0/2,
  
    ∫_ρ_0^R_0 dz  z^2ζ p_ζ,η^(α)(1,r,z)
        ≥ r^-η(1/2R_0^-η - Cρ_0^-η)
        = r^-η/4R_0^-η.

  Therefore, by Chapman–Kolmogorov for p_ζ,η^(α) and estimate (<ref>),
  
    p_ζ,η^(α)(3,r,s)
             ≥∫_ρ_0^R_0 dw  w^2ζ∫_ρ_0^R_0 dz  z^2ζ p_ζ,η^(α)(1,r,z)  p_ζ,η^(α)(1,z,w)  p_ζ,η^(α)(1,w,s) 
       ≥(rs)^-η/16R_0^-2ηinf_z,w∈(ρ_0,R_0) p_ζ,η^(α)(1,z,w).

  The infimum on the right is estimated with the help of (<ref>) by 
  
    inf_z,w∈(ρ_0,R_0) p_ζ,η^(α)(1,z,w)
        ≳_ρ_0,R_0inf_z,w∈(ρ_0,R_0) p_ζ^(α)(1,z,w)
        ≳_ρ_0,R_0 1,

  where we used (<ref>) (together with R_0≳1) to estimate the final infimum. Hence,
  
    p_ζ,η^(α)(3,r,s) ≳ (rs)^-η,    r∨ s<R_0/2.

  Thus, by the scaling property of p_ζ,η^(α) we obtain
  
    p_ζ,η^(α)(1,r,s)
        ≳_ζ,α,η,ρ_0,R_0 (rs)^-η,
           r∨ s<R_0/2· 3^1/α.

  This concludes the proof of Lemma <ref>.


We are now in position to prove the lower bounds in Theorem <ref>.


  
  Let ζ∈(-1/2,∞), α∈(0,2), and η∈(-α,0). Then,
  
    p_ζ,η^(α)(1,r,s)
        ≳_ζ,α,η (1∧ r^-η)(1∧ s^-η) p_ζ^(α)(1,r,s).




  By symmetry, it suffices to consider r≤ s. For 0<r<1/4, s>1, and 1/4≤ z≤1/2, Lemma <ref> and (<ref>) imply
  
    p_ζ,η^(α)(1,r,z) ≳ r^-η

  and
  
    p_ζ,η^(α)(1,s,z)
        ≳ p_ζ^(α)(1,s,z)
        ∼ p_ζ^(α)(1,s,r)

  where we used s≥ 2z≥ 2r, s>1, and (<ref>) in the final step. Hence, for any 0<r<1/4 and s>1, we have
  
    p_ζ,η^(α)(1,r,s)
             ≥∫_1/4^1/2 dz  z^2ζ p_ζ,η^(α)(1/2,r,z) p_ζ,η^(α)(1/2,z,s)
            ≳ r^-η p_ζ^(α)(1/2,r,s) 
       ∼ r^-η p_ζ^(α)(1,r,s).

  Finally, for r∧ s≥ 1/4, the claim follows from (<ref>), whereas for r∨ s≤ 1 the claim follows from Lemma <ref>.
  This concludes the proof of Lemma <ref> and therefore the proof of the lower bound on Theorem <ref>.





  §.§.§ Continuity



We now prove the continuity statement in Theorem <ref> for η<0.


  
  Let ζ∈(-1/2,∞), α∈(0,2), and η∈(-α,0). Then, for any fixed r>0, the function (0,∞)∋ s↦ p_ζ,η^(α)(t,r,s) is continuous.



  Fix r,s>0 and let z>0 converge to s. Recall q(r)=Ψ_ζ(η)r^-α with Ψ_ζ(η)<0. By Duhamel's formula,
  
    p_ζ,η^(α)(1,r,s) - p_ζ,η^(α)(1,r,z)
            = p_ζ^(α)(1,r,s) - p_ζ^(α)(1,r,z) 
          + ∫_0^1 dτ∫_0^∞ dw  w^2ζ p_ζ,η^(α)(1-τ,r,w)q(w)(p_ζ^(α)(τ,w,s) - p_ζ^(α)(τ,w,z)).

  For all sufficiently small ϵ, estimates (<ref>), (<ref>) and (<ref>) imply
  
    -∫_0^ϵ dτ∫_0^∞ dw  w^2ζ p_ζ,η^(α)(1-τ,r,w)q(w) p_ζ^(α)(τ,w,s) 
         ≤ -∫_0^ϵ∫_0^∞ dw  w^2ζ p_ζ^(α)(1-τ,r,w) q(w) p_ζ^(α)(τ,w,s) 
         ≲ -∫_0^ϵ∫_0^∞ dw  w^2ζ p_ζ^(α)(τ,w,s)q(w)
            ≲ϵ s^-α.

  Analogously, we have
  
    -∫_0^ϵ dτ∫_0^∞ dw  w^2ζ p_ζ,η^(α)(1-τ,r,w)q(q)p_ζ^(α)(τ,w,z)
        ≲ϵ z^-α.

  Next, for any ϵ<τ≤1 and w,s,z>0 with z→ s, we have
  p_ζ^(α)(τ,w,s)∼ p_ζ^(α)(τ,w,z) (e.g., as a consequence of (<ref>)). By dominated convergence,
  
    ∫_ϵ^1 dτ∫_0^∞ dw  w^2ζ p_ζ,η^(α)(1-τ,r,w)q(w)(p_ζ^(α)(τ,w,s)-p_ζ^(α)(τ,w,z)) →0   as z→ s.

  Combining all previous estimates establishes the claim.



  
  Let ζ∈(-1/2,∞), α∈(0,2), and η∈(-α,0). Then, p_ζ,η^(α)(t,r,s) is jointly continuous as a function of r,s,t>0.



  By scaling, it suffices to prove the continuity of p_ζ,η^(α)(1,r,s) with respect to r,s>0. As indicated in the proof of Lemma <ref>, we only need to verify
  
    ∫_0^1 dτ∫_0^∞ dw  w^2ζ| p_ζ,η^(α)(1-τ,r̃,w) p_ζ^(α)(τ,w,s̃) - p_ζ,η^(α)(1-τ,r,w) p_ζ^(α)(τ,w,s) | q(w) →0

  for any r,s,r̃,s̃>0 with r̃→ r and s̃→ s. In addition to (<ref>), we have
  
    -∫_1-ϵ^1 dτ∫_0^∞ dw  w^2ζ p_ζ,η^(α)(1-τ,r,w) q(w) p_ζ^(α)(τ,w,s) 
          = - ∫_0^ϵ dτ∫_0^∞ dw  w^2ζ p_ζ,η^(α)(τ,r,w) q(w) p_ζ^(α)(1-τ,w,s) 
         ≤ -∫_0^ϵ dτ∫_0^∞ dw  w^2ζ p_ζ^(α)(τ,r,w) q(w) p_ζ^(α)(1-τ,w,s)
          ≲ϵ r^-α

  by  the same argument as in the proof of (<ref>).

  For any ϵ<τ<1-ϵ and r,s,z,r̃,s̃>0 with r→r̃ and s→s̃, we have, by Lemma <ref> both p_ζ^(α)(τ,z,s̃) ∼ p_ζ^(α)(τ,z,s) and p_ζ,η^(α)(1-τ,r̃,z)∼ p_ζ^(α)(1-τ,r,z). Then, by dominated convergence,
  
    ∫_ϵ^1-ϵ dτ∫_0^∞ dw  w^2ζ|p_ζ,η^(α)(1-τ,r̃,w) p_ζ^(α)(τ,w,s̃) - p_ζ,η^(α)(1-τ,r,w)p_ζ^(α)(τ,w,s) | q(w) → 0

  for any r,s,r̃,s̃>0 with r̃→ r and s̃→ s. This concludes the proof of Proposition <ref> and in particular the continuity in Theorem <ref>.






§ SCHRÖDINGER PERTURBATIONS



In this Appendix, we discuss Schrödinger perturbations of transition densities needed in Theorem <ref>.
However, the setting is more general and the results are of independent interest; in particular we allow transition densities inhomogeneous in time.




 §.§ General setting and positive perturbations

Let X be a locally compact space with a countable base of open sets. Consider the Borel σ-algebra ℳ on X. Let m be a σ-finite measure on (X,ℳ). 
We further consider an arbitrary nonempty interval  on the real line , with the Borel σ-field ℬ and the Lebesgue measure ds.   We call × X the space-time.
The functions considered below are assumed to be (jointly) Borel measurable.
Let p be a 
transition density on X with time in . This means that p is a (jointly measurable) function on × X ×× X with values in [0,∞] such that p(s,x,t,y) =0 for s ≥ t and the Chapman-Kolmogorov equations hold for p, i.e.,

    p(s,x,t,y)=∫_X p(s,x,u,z) p(u,z,t,y)  m(dz),      s,t, u∈,   s<u<t,  x,y∈ X.

Below, we say that p is finite and strictly positive if 0<p(s,x,t,y)<∞ for all s,t∈, s<t and x,y∈ X.


We consider a measurable function (potential) q≥ 0 on × X and let 

    p_0(s,x,t,y)   := p(s,x,t,y)   
    
    p_n(s,x,t,y)    := ∫_s^t ∫_X p_n-1(s,x,u,z) q(u,z) p(u,z,t,y)(u,z) m(dz) du.

Here and below, s,t∈, s<t, and x,y∈ X.
We define 

    p^q(s,x,t,y) := ∑_n=0^∞  p_n(s,x,t,y)

and say that p^q is a Schrödinger perturbation of p by q. Specifically, the perturbation is positive since q≥ 0 and mass-creating since p^q≥ p. 
Of course, p^q is increasing in q≥ 0. 
By <cit.>, we have the Chapman-Kolmogorov equations

    p^q(s,x,t,y) = ∫_X p^q(s,x,u,z)  p^q(u,z,t,y) m(dz),   s<u<t.

By Tonelli,
p^q(s,x,t,y)  also satisfies the following perturbation (Duhamel) formulas

    p^q(s,x,t,y)    = p(s,x,t,y) + ∫_s^t ∫_X p^q(s,x,u,z) q(u,z) p(u,z,t,y) m(dz) du
       =
    p(s,x,t,y) + ∫_s^t ∫_X p(s,x,u,z) q(u,z) p^q(u,z,t,y) m(dz) du.

As straightforward as positive Schrödinger perturbations are, p^q may be infinite. For example, this is so for the Gaussian kernel p(s,x,t,y)=(4π (t-s))^-d/2e^-|y-x|^2/(4(t-s)) in , d>2, and q(x)=c|x|^-2, x∈, if (and only if) c >(d-2)^2/4, which underlies the famous result of Baras and Goldstein <cit.>. Needless to say, the example is 
time-homogeneous, meaning that p(s,x,s+t,y):=p_t(x,y),
s∈, t>0, 
and q(s,x):=q(x),
s∈, yield 
a (time-homogeneous) perturbation p_t^q(x,y):=p^q(0,x,t,y)=p^q(s,x,s+t,y), s∈, t>0.
It is in such a time-homogeneous setting that the results of this Appendix are used in Subsection <ref>.



 §.§ Negative perturbations


We now focus on functions q (or -q) which may be negative. To facilitate discussion, we assume p is sub-Markov, i.e., we assume ∫_X p(s,x,t,y) m(dy)≤1.
Then we can extend the definitions
(<ref>)—first to the case
of bounded, but not necessarily positive q:× X→.
Indeed, in this case we can estimate the integrals by using |q| and q_∞:=sup_s∈,x∈ X |q(s,x)|.
In particular, by induction, we get

    |p_n(s,x,t,y)|≤q_∞^n (t-s)^n/n! p(s,x,t,y),   n=0,1,….
 
We see that the perturbation series converges absolutely and

    |p^q(s,x,t,y)|≤ p^|q|(s,x,t,y)≤ e^q_∞ (t-s)p(s,x,t,y).

By <cit.>, p^q satisfies Chapman-Kolmogorov equations.
By Fubini, p^q  also satisfies the Duhamel formulas (<ref>).

Next, we shall define and study rather general negative
perturbations of transition densities. To this end, we consider integral kernels P(s,x,B) on the space-time, where s∈, x∈ X, and B is a (measurable) subset of × X. Thus, by <cit.>,

    (s,x)↦ P((s,x),B)  B, 
        B↦ P((s,x),B) × X  s  x.

Then,
Pf(s,x):=∫ f(t,y)P(s,x,dt dy)≥ 0 
is measurable for every measurable function f≥ 0 on × X. We will identify the kernel P with the operator P. In fact, every additive, positively homogeneous and monotone operator on nonnegative measurable functions
is an integral kernel, see <cit.>.


For instance, our transition density p defines, for functions f≥ 0 on × X,

    Pf(s,x) := ∫_∫_X p(s,x,t,y) f(t,y)m(dy)dt,   s∈,  x∈ X.

Here is a variation on <cit.>, which points out a sub-Markov resolvent <cit.> associated with the kernel P in (<ref>). The proof applies to  in the same way as to , with only minor adjustments; we provide it for the reader's convenience.


  
  For (<ref>), a sub-Markov resolvent P^λ, λ>0, exists with sup_λ>0 P^λ = P.



  For λ≥ 0, we let
  
    p^-λ(s,x,t,y) :=  e^-λ(t-s)p(s,x,t,y),   s,t∈,  x,y∈ X,

  and, for f≥ 0, 
  
    P^λ f(s,x) := ∫_∫_X p^-λ(s,x,t,y) f(t,y)   m(dy)   dt,    s∈,   x∈ X.

  For clarity, note that p^-λ(s,x,t,y) = 0 for s≥ t, so in (<ref>), we have ∫_=∫_s^S, where S:=sup. Of course, sup_λ>0 P^λ = P as kernels. Furthermore, for s∈, x∈ X, λ>0,
  
    λ P^λ (s,X)  = λ∫_∫_X p^-λ(s,x,t,y)   m(dy)   dt ≤λ∫_s^∞ e^-λ(t-s)  dt ≤ 1.

  Finally, by the Chapman-Kolmogorov equations, for λ>μ, s∈, x∈ X, f≥ 0, 
  
    P^μ P^λ f(s,x) +(λ-μ)^-1P^λf(s,x)
        
    
        =    ∫_∫_X  ∫_∫_X p^-μ(s,x,u,z)p^-λ(u,z,t,y)f(t,y)   m(dy)   dt   m(dz)   du+(λ-μ)^-1P^λf(s,x)     
    
        =    ∫_s^S ∫_X p(s,x,t,y) f(t,y) e^μ s e^-λ t∫_s^t e^-u(μ-λ)  du   m(dy)   dt+(λ-μ)^-1P^λf(s,x)     
    
        =      (λ-μ)^-1∫_X ∫_s^S p(s,x,t,y) f(t,y) ( e^-μ(t-s) - e^-λ(t-s))  m(dy)   dt    
        + (λ-μ)^-1P^λf(s,x)
         =  (λ-μ)^-1 P^μf(s,x).

  Thus, P^λ, λ>0, is a sub-Markov resolvent.


The notation p^-λ  in (<ref>) agrees with that in (<ref>). In particular, the resolvent equation (<ref>) with μ=0 is a variant of the Duhamel formula for q≡ -λ.

The following lemma is crucial for handling negative perturbations.

  
  If f,g≥ 0, Pf≤ Pg+1 on {f>0}, then Pf≤ Pg+1 on × X. In particular, if P|h|<∞ and Ph≤ 1 on {h>0}, then Ph≤ 1 on × X. 


Here, of course, {f>0}:={(s,x)∈× X: f(s,x)>0}. We will call the first implication in Lemma <ref> the complete maximum principle, abbreviated as CMP. It is a variant of <cit.>, but here we do not assume that f is bounded and Pf finite. CMP may be also viewed as a variant of the domination principle <cit.>, but discussing this connection would take longer than the proof below.


If f≥ 0 is bounded, then Lemma <ref> and the proof of <cit.> give the implication.
In the case of general f≥ 0 with
Pf≤ Pg+1 on {f>0}, we let f_n:=f∧ n for 
n∈.
On {f_n>0}⊂{f>0}, we have Pf_n≤ P f≤ Pg+1, so by the first part of the proof,  Pf_n≤ Pg+1 on × X. Letting n→∞, we get the first statement. For the second one, we let f:=h_+, g:=h_-.





  
  Let q≥ 0 be bounded. If p is a strictly positive finite sub-Markov transition density then so is p^-q. Furthermore, p^-q≤ p.



  Let y∈ X, t∈, and _t:={s∈: s<t}. 
  For s,u∈_t, s<u, x,z∈ X, let
  
    h(s,x,u,z) := p(s,x,u,z)p(u,z,t,y)/p(s,x,t,y).

  The function is a transition density on the space-time _t× X, so by Lemma <ref>,
  
    P_t,yf(s,x)    := ∫_s^t ∫_X p(s,x,u,z)p(u,z,t,y)/p(s,x,t,y)f(u,z)  m(dz)  du,    s<t,   x ∈ X,

is an integral kernel that  satisfies CMP. Let
  
    P f(s,x) := ∫_s^t ∫_X p(s,x,u,z)p(u,z,t,y)/p(s,x,t,y)f(u,z)q(u,z) m(dz)du.

  Then P satisfies CMP, too. Indeed, let f,g≥ 0 on _t× X and P f ≤P g + 1 on {f >0}. Let f' :=  f q, g':=g q. Since {f'>0}⊂{f>0}, we have P_t,y f'≤ P_t,yg'+1 on {f'>0}. By CMP for P_t,y, Pf=P_t,y f' ≤ P_t,yg'+ 1=Pg+1 on _t× X, which gives CMP for P.
  Next, let h(s,x) := p^-q(s,x,t,y)/p(s,x,t,y). By (<ref>), P |h|<∞. By (<ref>),
  
    h = 1 - P h.

  Of course, if h(s,x) > 0, then P h(s,x) < 1. By CMP for P, we get P h ≤ 1, so h ≥ 0 and 0≤ p^-q≤ p everywhere.
  Furthermore, we note that for all s,t∈ and x,y∈ X,
  
    p^-q(s,x,t,y)
           ≥ p(s,x,t,y)-∑_n=1^∞ p_n(s,x,t,y)
       ≥(2-e^q_∞ (t-s))p(s,x,t,y)>0

  if 0<t-s< ln 2/q_∞. By (<ref>), p^-q(s,x,t,y)>0 for all s,t∈, s<t, x,y∈ X.


We say that the perturbation of p in Lemma <ref> is negative since -q≤ 0 and mass-decreasing since p^-q≤ p. The main point of Lemma <ref> is that p^-q≥ 0, in fact, that p^-q is strictly positive for (bounded) q≥ 0. This  is a similar phenomenon as e^-x> 0 for x∈ [0,∞).
In passing, we also note that for bounded functions q,ρ≥ 0 on × X, by <cit.>,

    p^-q-ρ=(p^-q)^-ρ.



In the reminder of this Appendix, we analyze p^-q for possibly unbounded functions q≥ 0 on space-time. 
As we shall see, the following integrability assumption,

    p_1(s,x,t,y)=∫_s^t ∫_X p(s,x,u,z)q(u,z)p(u,z,t,y)m(dz)du<∞

for all s,t∈, s<t, x,y∈ X, will suffice for construction of p^-q. Note that 

    p^q<∞

implies (<ref>), but not conversely; see the above discussion of the result of Baras and Goldstein. 
The advantage of (<ref>) is that for arbitrary λ>0, it holds for q≥ 0 if and only if it holds for λ q.
For the sake of
the discussion following Definition <ref>, we note that (<ref>) holds for the time-homogeneous semigroup p(s,x,t,y):=p_ζ^(α)(t-s,x,y) and  potential  q(s,x)=x^-α there,
thanks to 0<α<2ζ+1. This is because p^λ q<∞ for small λ=Ψ_ζ(η)>0, i.e., small  η>0; see Lemma <ref>. 


Here is the main result of this Appendix.

  
  Let q:× X→ [0,∞]. Let p be a sub-Markov, strictly positive, finite transition density, and assume (<ref>).
  Then there is a strictly positive transition density p^-q≤ p satisfying (<ref>) and (<ref>).




  For n∈, we let q_n:=q∧ n. Then each p^-q_n are sub-Markov, strictly positive, finite transition densities decreasing in n. Let 
  
    p^-q:=inf_n p^-q_n.

  Of course, 0≤ p^-q≤ p. By (<ref>) and the dominated convergence theorem, we get (<ref>) and (<ref>).
  The strict positivity of p^-q follows from the lower bound
  
    p^-q (s,x,t,y) ≥ p(s,x,t,y) exp[-p_1(s,x,t,y)p(s,x,t,y)],

  where q≥ 0, s,t ∈, s<t, and x,y ∈ X. 
  For the proof of (<ref>), we first let q≥ 0 be bounded. Then the function (0,∞)∋λ↦ h(λ):=p^-λ q(s,x,t,y) is completely monotone, meaning that (-1)^n h^(n)(λ)≥ 0, n=0,1,…, λ>0.
  Indeed, (-1)^n h^(n)(λ)=n! (p^-λ q)_n(s,x,t,y) by <cit.>, which is nonnegative in view of Lemma <ref> and (<ref>).
  Since completely monotone functions are logarithmically convex, we get
  
    p^-q(s,x,t,y)=h(1)
            = h(0)exp[∫_0^1 (ln h(λ))'dλ]
          ≥ h(0)exp[∫_0^1 h'(0)/h(0)dλ] 
        = p(s,x,t,y)exp[-p_1(s,x,t,y)/p(s,x,t,y)].

  For unbounded q≥ 0, we use this, (<ref>), and (<ref>) and we let n→∞.



  (1) Estimate (<ref>) strengthens <cit.> and <cit.>.
  

  (2) In the time-homogeneous setting, (<ref>) yields

    p^-q_t (x,y) ≥ p_t(x,y) exp[-∫_0^t∫_X p_s(x,z)q(z)p_t-s(z,y)m(dz) ds/p_t(x,y)].

  Here, of course, p_t is sub-Markov, q≥ 0, and the numerator is assumed finite.










We conclude this appendix with the following example,
which illustrates the importance of CMP, and therefore the Chapman–Kolmogorov equations, for
the positivity of p^-q.


  Let X be the space containing only one point and m(dz) be the Dirac delta at this point. Let p(s,t) = (t-s)_+ and q(s) ≡ 1 (to simplify notation, we skip the space variables from the notation). Then,
  
    p_n(s,t) = (t-s)_+^2n+1/(2n+1)!,

  therefore
  
    p^q(s,t) = sinh(t-s)_+          p^-λ q(s,t) = sin(√(λ)(t-s)_+)/√(λ).

  Clearly, p^-λ q(s,t) takes on negative values, too. This does not contradict our findings because the operator P f(s) := ∫_s^∞ p(s,t) f(t) dt does not satisfy CMP, as may be verified by considering the function f(s) = - 1_[0,1)(s) + 2 · 1_[1,2](s).










'
KMV+18

[ABG+23]Armstrongetal2023
Gavin Armstrong, Krzysztof Bogdan, Tomasz Grzywny, Łukasz Leżaj, and
  Longmin Wang.
Yaglom limit for unimodal Lévy processes.
Ann. Inst. Henri Poincaré Probab. Stat., 59(3):1688–1721,
  2023.

[BD23]BuiDAncona2023
The Anh Bui and Piero D'Ancona.
Generalized Hardy operators.
Nonlinearity, 36(1):171–198, 2023.

[BDG11]Bruneauetal2011
Laurent Bruneau, Jan Dereziński, and Vladimir Georgescu.
Homogeneous Schrödinger operators on half-line.
Ann. Henri Poincaré, 12(3):547–590, 2011.

[BDK16]Bogdanetal2016
Krzysztof Bogdan, Bartłomiej Dyda, and Panki Kim.
Hardy inequalities and non-explosion results for semigroups.
Potential Anal., 44(2):229–247, 2016.

[BG60]BlumenthalGetoor1960
R. M. Blumenthal and R. K. Getoor.
Some theorems on stable processes.
Trans. Amer. Math. Soc., 95:263–273, 1960.

[BG84]BarasGoldstein1984
Pierre Baras and Jerome A. Goldstein.
The heat equation with a singular potential.
Trans. Amer. Math. Soc., 284(1):121–139, 1984.

[BGJP19]Bogdanetal2019
Krzysztof Bogdan, Tomasz Grzywny, Tomasz Jakubowski, and Dominika Pilarczyk.
Fractional Laplacian with Hardy potential.
Comm. Partial Differential Equations, 44(1):20–50, 2019.

[BH86]BliedtnerHansen1986
J. Bliedtner and W. Hansen.
Potential Theory.
Universitext. Springer-Verlag, Berlin, 1986.
An Analytic and Probabilistic Approach to Balayage.

[BHJ08]Bogdanetal2008
Krzysztof Bogdan, Wolfhard Hansen, and Tomasz Jakubowski.
Time-dependent Schrödinger perturbations of transition
  densities.
Studia Math., 189(3):235–254, 2008.

[BHNV10]Betancoretal2010
Jorge J. Betancor, Eleonor Harboure, Adam Nowak, and Beatriz Viviani.
Mapping properties of fundamental operators in harmonic analysis
  related to Bessel operators.
Studia Math., 197(2):101–140, 2010.

[BJ07]BogdanJakubowski2007
Krzysztof Bogdan and Tomasz Jakubowski.
Estimates of heat kernel of fractional Laplacian perturbed by
  gradient operators.
Comm. Math. Phys., 271(1):179–198, 2007.

[BJK20]Bogdanetal2020
Krzysztof Bogdan, Sven Jarohs, and Edyta Kania.
Semilinear Dirichlet problem for the fractional Laplacian.
Nonlinear Anal., 193:111512, 20, 2020.

[BJKP23]Bogdanetal2023
Krzysztof Bogdan, Tomasz Jakubowski, Panki Kim, and Dominika Pilarczyk.
Self-similar solution for Hardy operator.
J. Funct. Anal., 285(5):Paper No. 110014, 40, 2023.

[BKW01]Bileretal2001
Piotr Biler, Grzegorz Karch, and Wojbor A. Woyczyński.
Critical nonlinearity exponent and self-similar asymptotics for
  Lévy conservation laws.
Ann. Inst. H. Poincaré C Anal. Non Linéaire,
  18(5):613–637, 2001.

[BM23]BogdanMerz2023S
Krzysztof Bogdan and Konstantin Merz.
Subordinated Bessel heat kernels.
arXiv e-prints, page arXiv:2308.15026, August 2023.

[BM24a]BogdanMerz2024
Krzysztof Bogdan and Konstantin Merz.
Ground state representation for the fractional Laplacian with
  Hardy potential in angular momentum channels.
J. Math. Pures Appl. (9), 186:176–204, 2024.

[BM24b]BogdanMerz2024H
Krzysztof Bogdan and Konstantin Merz.
Heat kernel bounds for the fractional Laplacian with Hardy
  potential in angular momentum channels.
In preparation, 2024.

[BN22]BuiNader2022
The Anh Bui and Georges Nader.
Hardy spaces associated to generalized Hardy operators and
  applications.
NoDEA Nonlinear Differential Equations Appl., 29(4):Paper No.
  40, 40, 2022.

[CKSV20]Choetal2020
Soobin Cho, Panki Kim, Renming Song, and Zoran Vondraček.
Factorization and estimates of Dirichlet heat kernels for non-local
  operators with critical killings.
J. Math. Pures Appl. (9), 143:208–256, 2020.

[CPR13]Chaumontetal2013
Loïc Chaumont, Henry Pantí, and Víctor Rivero.
The Lamperti representation of real-valued self-similar Markov
  processes.
Bernoulli, 19(5B):2494–2523, 2013.

[CS24]ChoSong2024
Soobin Cho and Renming Song.
Fractional Laplacian with supercritical killings.
arXiv e-prints, page arXiv:2403.03298, March 2024.

[DG21]DerezinskiGeorgescu2021
Jan Dereziński and Vladimir Georgescu.
On the domains of Bessel operators.
Ann. Henri Poincaré, 22(10):3291–3309, 2021.

[DLMF23]NIST:DLMF
NIST Digital Library of Mathematical Functions.
<https://dlmf.nist.gov/>, Release 1.1.9 of 2023-03-15, 2023.
F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider,
  R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and
  M. A. McClain, eds.

[FG16]FrankGeisinger2016
Rupert L. Frank and Leander Geisinger.
Refined semiclassical asymptotics for fractional powers of the
  Laplace operator.
J. Reine Angew. Math., 712:1–37, 2016.

[Fit00]Fitzsimmons2000
P. J. Fitzsimmons.
Hardy's inequality for Dirichlet forms.
J. Math. Anal. Appl., 250(2):548–560, 2000.

[FLS08]Franketal2008H
Rupert L. Frank, Elliott H. Lieb, and Robert Seiringer.
Hardy-Lieb-Thirring inequalities for fractional Schrödinger
  operators.
J. Amer. Math. Soc., 21(4):925–950, 2008.

[FMS21]Franketal2021
Rupert L. Frank, Konstantin Merz, and Heinz Siedentop.
Equivalence of Sobolev norms involving generalized Hardy
  operators.
International Mathematics Research Notices, 2021(3):2284–2303,
  February 2021.

[FMS23a]Franketal2023T
Rupert L. Frank, Konstantin Merz, and Heinz Siedentop.
The Scott conjecture for large Coulomb systems: a review.
Lett. Math. Phys., 113(1):Paper No. 11, 2023.

[FMS23b]Franketal2023
Rupert L. Frank, Konstantin Merz, and Heinz Siedentop.
Relativistic strong Scott conjecture: a short proof.
In Density functionals for many-particle systems—mathematical
  theory and physical applications of effective equations, volume 41 of Lect. Notes Ser. Inst. Math. Sci. Natl. Univ. Singap., pages 69–79. World
  Sci. Publ., Hackensack, NJ, March [2023] 2023.

[FMSS20]Franketal2020P
Rupert L. Frank, Konstantin Merz, Heinz Siedentop, and Barry Simon.
Proof of the strong Scott conjecture for Chandrasekhar atoms.
Pure Appl. Funct. Anal., 5(6):1319–1356, December 2020.

[GT21]GrzywnyTrojan2021
Tomasz Grzywny and Bartosz Trojan.
Subordinated Markov processes: sharp estimates for heat kernels and
  Green functions.
arXiv e-prints, page arXiv:2110.01201, October 2021.

[HB23]HansenBogdan2023
Wolfhard Hansen and Krzysztof Bogdan.
Positive harmonically bounded solutions for semi-linear equations.
arXiv e-prints, page arXiv:2212.13999, December 2023.

[JKS23]Jakubowskietal2023
Tomasz Jakubowski, Kamil Kaleta, and Karol Szczypkowski.
Bound states and heat kernels for fractional-type Schrödinger
  operators with singular potentials.
Comm. Math. Phys., 403(2):795–828, 2023.

[JKS24]Jakubowskietal2024
Tomasz Jakubowski, Kamil Kaleta, and Karol Szczypkowski.
Relativistic stable operators with critical potentials.
Doc. Math., 29(1):209–245, 2024.

[JM23a]JakubowskiMaciocha2023G
Tomasz Jakubowski and Paweł Maciocha.
Ground-state representation for the fractional Laplacian on the
  half-line.
Probab. Math. Statist., 43(1):83–108, 2023.

[JM23b]JakubowskiMaciocha2023
Tomasz Jakubowski and Paweł Maciocha.
Heat kernel estimates of fractional Schrödinger operators with
  Hardy potential on half-line.
arXiv e-prints, page arXiv:2311.06683, November 2023.

[JW20]JakubowskiWang2020
Tomasz Jakubowski and Jian Wang.
Heat kernel estimates of fractional Schrödinger operators with
  negative Hardy potential.
Potential Anal., 53(3):997–1024, 2020.

[KMV+18]Killipetal2018
R. Killip, C. Miao, M. Visan, J. Zhang, and J. Zheng.
Sobolev spaces adapted to the Schrödinger operator with
  inverse-square potential.
Math. Z., 288(3-4):1273–1298, 2018.

[KRS21]Kyprianouetal2021
Andreas E. Kyprianou, Victor Rivero, and Weerapat Satitkanitkul.
Stable Lévy processes in a cone.
Ann. Inst. Henri Poincaré Probab. Stat., 57(4):2066–2099,
  2021.

[KSS21]Kinzebulatovetal2021
D. Kinzebulatov, Yu. A. Semënov, and K. Szczypkowski.
Heat kernel of fractional Laplacian with Hardy drift via
  desingularizing weights.
J. Lond. Math. Soc. (2), 104(4):1861–1900, 2021.

[Mer21]Merz2021
Konstantin Merz.
On scales of Sobolev spaces associated to generalized Hardy
  operators.
Math. Z., 299(1):101–121, 2021.

[MNS18]Metafuneetal2018
G. Metafune, L. Negro, and C. Spina.
Sharp kernel estimates for elliptic operators with second-order
  discontinuous coefficients.
J. Evol. Equ., 18(2):467–514, 2018.

[MNS23]Metafuneetal2023
Giorgio Metafune, Luigi Negro, and Chiara Spina.
Elliptic and parabolic problems for a Bessel-type operator.
In Recent Advances in Mathematical Analysis—Celebrating the
  70th Anniversary of Francesco Altomare, Trends Math., pages 397–424.
  Birkhäuser/Springer, Cham, [2023] 2023.

[MS22]MerzSiedentop2022
Konstantin Merz and Heinz Siedentop.
Proof of the strong Scott conjecture for heavy atoms: the Furry
  picture.
Annales Henri Lebesgue, 5:611–642, 2022.

[MSŻ16]Maleckietal2016
Jacek Małecki, Grzegorz Serafin, and Tomasz Żórawik.
Fourier-Bessel heat kernel estimates.
J. Math. Anal. Appl., 439(1):91–102, 2016.

[PS21]PatieSavov2021
Pierre Patie and Mladen Savov.
Spectral expansions of non-self-adjoint generalized Laguerre
  semigroups.
Mem. Amer. Math. Soc., 272(1336):vii+182, 2021.

[RS75]ReedSimon1975
Michael Reed and Barry Simon.
Methods of Modern Mathematical Physics. II. Fourier
  Analysis, Self-Adjointness.
Academic Press [Harcourt Brace Jovanovich Publishers], New York,
  1975.

[SSV12]Schillingetal2012
René L. Schilling, Renming Song, and Zoran Vondraček.
Bernstein Functions, volume 37 of De Gruyter Studies in
  Mathematics.
Walter de Gruyter & Co., Berlin, second edition, 2012.
Theory and Applications.

[SWW22]Songetal2022
Renming Song, Peixue Wu, and Shukun Wu.
Heat kernel estimates for non-local operator with multisingular
  critical killing.
arXiv e-prints, page arXiv:2203.03891, March 2022.

[V1́8]Vazquez2018
Juan Luis Vázquez.
Asymptotic behaviour for the fractional heat equation in the
  Euclidean space.
Complex Var. Elliptic Equ., 63(7-8):1216–1231, 2018.




