#  RealTimeTransport: An open-source C++ library for quantum transport
simulations in the strong coupling regime

Konstantin Nestmann<EMAIL_ADDRESS>Martin Leijnse  Division
of Solid State Physics and NanoLund, Lund University, Lund, Sweden  Maarten R.
Wegewijs  Institute for Theory of Statistical Physics, RWTH Aachen, Aachen,
Germany  JARA-FIT, Aachen, Germany  Peter Grünberg Institut, Forschungszentrum
Jülich, Jülich, Germany

###### Abstract

The description of quantum transport in the strong system-reservoir coupling
regime poses a significant theoretical and computational challenge that
demands specialized tools for accurate analysis. RealTimeTransport is a new
open-source C++ library that enables the computation of both stationary and
transient transport observables for generic quantum systems connected to
metallic reservoirs. It computes the Nakajima-Zwanzig memory kernels for both
dynamics and transport in real-time going beyond traditional expansions in the
bare system-reservoir couplings. Currently, several methods are available: (i)
A renormalized perturbation theory in leading and next-to-leading order which
avoids the low-temperature breakdown that limits the traditional theory. (ii)
Starting from this well-behaved reference solution a 2- and 3-loop self-
consistent renormalization-group transformation of the memory kernels is
implemented. This allows refined quantitative predictions even in the presence
of many body resonances, such as the Kondo enhancement of cotunneling. This
paper provides an overview of the theory, the architecture of
RealTimeTransport and practical demonstrations of the currently implemented
methods. In particular, we analyze the stationary transport through a serial
double quantum dot and showcase for the $T=0$ interacting Anderson model the
complete time-development of single-electron tunneling (SET), cotunneling-
assisted SET (CO-SET) and inelastic cotunneling resonances throughout the
entire gate-bias stability diagram. We discuss the range of applicability of
the implemented methods and benchmark them against other advanced approaches.

##  I Introduction

The accurate description of quantum transport phenomena is fundamental to the
understanding of quantum dot devices. These are interesting platforms to study
fundamental aspects of open quantum systems, featuring a rich interplay of
non-equilibrium physics and many-body effects, as well as for various
(quantum) technologies. CODE(0x560026ecf290)hough their stationary transport
properties are by now fairly well-understood, their dynamical aspects remain
challenging.

From a methodological perspective an accurate description in extended
parameter regimes is challenging. Non-equilibrium Green’s function methods are
well-suited to study systems at low temperature, but struggle to treat strong
Coulomb interactions which are typically incorporated perturbatively or in
mean-field approximation. CODE(0x560026ecf290)s has been improved by various
means, in particular by applying renormalization group (RG) methods within
this framework. CODE(0x560026ecf290)he other hand, a plethora of quantum
master equation (QME) approaches exist, CODE(0x560026ecf290)ch treat the
quantum dot energy scales non-perturbatively. These instead assume that the
system-reservoir coupling is weak, which prohibits their usage at temperatures
below the energy scale set by the tunnel rate. Just like Green’s function
approaches, master equation-like methods can be improved by applying RG ideas,
but in this case to approach stronger system-reservoir couplings and lower
temperatures. CODE(0x560026ecf290) above mentioned methods are semi-
analytical, which often aids the understanding of the results, and the work
reported here follows this line. Fully numerical approaches such as scattering
state numerical renormalization groups, CODE(0x560026ecf290)me-dependent)
density matrix renormalization groups CODE(0x560026ecf290)uantum monte carlo
methods CODE(0x560026ecf290)st to name a few, have also been successfully
applied.

In this paper we present RealTimeTransport, a C++ library which implements
several methods to analyze both transient and stationary quantum transport
quantities based on the Nakajima-Zwanzig memory kernel of the open system.
111The library is available at https://github.com/kn-code/RealTimeTransport As
the name suggests, the package is consistently based on the time-domain
representation of the dynamics and transport. It thereby maintains a useful
connection to the intuition of traditional quantum master equations (QMEs):
inside the “black box” one can still identify rates for state transitions and
Bloch coherence-vector dynamics. At the same time it fixes shortcomings of
these QMEs by fully accounting for renormalization and memory effects. Key
technical steps of our approach were heavily inspired by the so-called
“$E$-flow” RG formulated in frequency domain ($E$), which was successfully
used to analytically and numerically describe the non-equilibrium Kondo effect
in the stationary and time-dependent case. CODE(0x560026ecf290)nstead remain
in the time domain noting that there the memory kernel is typically well-
behaved and can be efficiently resolved using well-developed Chebyshev
interpolation techniques. CODE(0x560026ecf290)

Roughly speaking our approach works as follows: In ordinary bare perturbation
theory one initially eliminates all dissipation from the problem by setting
the coupling to the reservoir to zero. One then computes the corrections to
this “too small dissipation” to increase it. In our approach we do exactly the
opposite: We initially analytically solve the problem in the physical limit of
infinite temperature leading to “too large dissipation” and then compute
corrections to this. In prior work this was done by a “$T$-flow” from
infinitely high temperature to the low temperature of interest.
CODE(0x560026ecf290)s involved solving a differential RG equation with respect
to the physical temperature $T$, which at low $T$ approaches a fixed point (in
the space of time-dependent memory kernels and vertex functions).

In the present work we perform this computation in a new way. We exploit that
the “too large dissipation” computed at $T=\infty$ allows to construct an
ansatz which is well-behaved even at $T=0$. This reference solution
incorporates all $T=\infty$ effects and is the leading term of a much
simplified renormalized perturbation series. In particular, this ansatz is
very similar to the bare perturbation expression whose virtual intermediate
evolutions are regularized as $i0^{+}\to\Sigma_{\infty}$, i.e., by replacing
the problematic zero-dissipation by the $T=\infty$ Lindblad generator.
Importantly, this superoperator-valued regularization is systematically
derived from a physical limit and not put in “by hand”. Our main result is
that from this ansatz one can initiate a discrete RG iteration which converges
stepwise to the same fixed point as the previous continuous $T$-flow, thus
reducing the “too large dissipation” to the desired correct value.
Importantly, for the low temperature of interest this discrete iteration turns
out to be more efficient then the continous $T$-flow.

Although the RG aspects of the outlined method make it stand out, it is firmly
routed in the well-developed line of quantum master equation approaches for
which various numerical packages have been published. For example, the QuTiP
(Quantum Toolbox in Python) package CODE(0x560026ecf290)lements
phenomenological and leading order approaches in the system-reservoir
coupling, such as the Lindblad, Bloch-Redfield and Floquet-Markov QMEs. QmeQ
(Quantum master equations for Quantum dot transport calculations)
CODE(0x560026ecf290)o implements several leading order methods, but
additionally provides next-to-leading order approximations of the stationary
state and particle/energy current through a consistent expansion of the memory
kernel or the second-order von Neumann approach.

The paper is structured as follows: In Sec. II we first present the time-space
perturbation theory for the memory kernel of the density operator and then
develop the above mentioned renormalized version of this. Based on this we set
up the new discrete iterative RG transformation and comment on some
implementation details in Sec. III. In Sec. IV we then showcase the practical
library usage for several examples. In particular, we analyze the stationary
transport through a serial double quantum dot and present the time evolution
of the non-equilibrium interacting Anderson model in a magnetic field at zero
temperature. This showcases the emergence of different strongly-correlated
phenomena at distinct time-scales which is mapped out for all applied bias and
gate voltages. We note that the _stationary_ limit of this problem has by now
been experimentally probed in detail in many types of quantum dots, ranging
from semi-conductor heterostructures, carbon-based materials to molecular and
atomic quantum dots. CODE(0x560026ecf290)ally, we benchmark our methods
against the density matrix renormalization group and quantum Monte Carlo
approach.

##  II Computation of the memory kernel

In the following we set $\hbar=k_{B}=e=1$.

###  II.1 Model and notation

We consider a quantum dot system connected to several non-interacting electron
reservoirs in the thermodynamic limit. Each of the reservoirs labeled by $r$
is assumed to be initially in a grand canonical equilibrium with temperature
$T_{r}$ and chemical potential $\mu_{r}$. Thus the Hamiltonian describing the
isolated reservoirs is given by

$\displaystyle H_{R}=\sum_{r\nu}\int
d\omega(\omega+\mu_{r})a^{\dagger}_{r\nu}(\omega)a_{r\nu}(\omega),$ (1)

and the initial state of the reservoirs is

$\displaystyle\rho_{R}=\prod_{r}\frac{1}{Z_{r}}e^{-\beta_{r}\sum_{\nu}\int
d\omega\omega a^{\dagger}_{r\nu}(\omega)a_{r\nu}(\omega)}.$ (2)

Here $\nu$ denotes a generic channel index of the reservoirs, which is often
just the spin $\sigma$. The field operators $a^{\dagger}_{r\nu}$ ($a_{r\nu}$)
create (destroy) an electron with channel index $\nu$ in reservoir $r$ and
anti-commute as

$\displaystyle\big{\\{}a_{r\nu}(\omega),a^{\dagger}_{r^{\prime}\nu^{\prime}}(\omega^{\prime})\big{\\}}$
$\displaystyle=\delta_{rr^{\prime}}\delta_{\nu\nu^{\prime}}\delta(\omega-\omega^{\prime}),$
(3a)
$\displaystyle\big{\\{}a_{r\nu}(\omega),a_{r^{\prime}\nu^{\prime}}(\omega^{\prime})\big{\\}}$
$\displaystyle=0.$ (3b)

It is useful to introduce an additional particle-hole index $\eta$ and define

$\displaystyle a_{\eta r\nu\omega}$
$\displaystyle\coloneqq\left\\{\begin{array}[]{lrl}a^{\dagger}_{r\nu}(\omega)&\text{for}&\eta=+\\\
a_{r\nu}(\omega)&\text{for}&\eta=-\end{array}\right..$ (6)

For notational convenience we collect all indices into a multi-index written
as a number,

$\displaystyle 1\equiv\eta,r,\nu,\omega;\quad\bar{1}\equiv-\eta,r,\nu,\omega,$
(7)

where the overbar indicates the inversion of the particle-hole index. This
allows us to rewrite the anti-commutation relations (3) compactly as

$\displaystyle\\{a_{1},a_{2}\\}=\delta_{1\bar{2}}.$ (8)

Electrons can tunnel between the reservoirs and quantum dot system, which is
captured by the tunneling Hamiltonian

$\displaystyle H_{T}=\sum_{r\nu l}\int d\omega\sqrt{\lambda_{r\nu}}t_{r\nu
l}d^{\dagger}_{l}a_{r\nu}(\omega)+\text{h.c.}$ (9)

Here $d_{l}$ denotes a field operator of the quantum dot system and $l$ a
(multi)index labeling the single particle states, including the spin and
orbital quantum numbers. Importantly, the Fock space can always be constructed
such that the field operators $d_{l}$ of the system and $a_{r\nu}(\omega)$ of
the reservoir commute, rather than anti-commute. This is explained in Sec. II.
A. of Ref. [Saptsov and Wegewijs, 2014]: one transforms from the conventional
anti-commutation to commutation using the fermion parity operators of the
system and reservoir, respectively. This is crucial for tracing out the
fermionic reservoir degrees of freedom later on, because it avoids keeping
track of various minus signs.

We assume that the density of states $\lambda_{r\nu}$ and tunneling amplitudes
$t_{r\nu l}$ are energy independent (wideband limit) and real-valued. They
will later enter the transport equations via the spectral density, which we
define as

$\displaystyle\Gamma_{rll^{\prime}}\coloneqq
2\pi\sum_{\nu}\lambda_{r\nu}t_{r\nu l}t_{r\nu l^{\prime}}.$ (10)

The Hamiltonian $H$ describing the isolated quantum dot system is only
constrained by the assumption that it should commute with the fermion parity
(as any physical Hamiltonian should). CODE(0x560026ecf290)s is because the
fermion parity superselection principle is explicitly incorporated into the
superfermionic perturbation theory introduced later. Thus, the total
Hamiltonian we consider is given by

$\displaystyle H_{\text{tot}}=H+H_{R}+H_{T}.$ (11)

###  II.2 Superfermionic perturbation theory

For initial total states that factorize,
$\rho^{0}_{\text{tot}}=\rho_{0}\otimes\rho_{R}$, the dynamics of the reduced
density operator $\rho(t)\coloneqq\operatorname{Tr}_{R}\rho_{\text{tot}}(t)$
can be described using the superoperator-valued propagator $\Pi(t)$ as
$\rho(t)=\Pi(t)\rho_{0}$. The dynamics of $\Pi(t)$ follows the time-nonlocal
Nakajima-Zwanzig quantum master equation, CODE(0x560026ecf290)Π(t) = -i LΠ(t)
- i ∫_0^t ds K(t-s) Π(s), where $\mathcal{L}\bullet=[H,\bullet]$ denotes the
bare Liouvillian, $\mathcal{K}$ the memory kernel and the bullet an arbitrary
operator argument. Our main goal is to compute the memory kernel
$\mathcal{K}$. The perturbative expansion of $\mathcal{K}$ in the system-
environment coupling is well-established, CODE(0x560026ecf290) most commonly
done in frequency space and in particular done at $\omega=0$, where the
stationary state can be extracted. However, in order to implement the
renormalization group transformation presented in Sec. II.4, it is necessary
to resolve the memory kernel either for all frequencies or alternatively for
all times. Because the memory kernel as a function of time is typically well
behaved and (at finite temperatures) exponentially decaying, it can be
efficiently interpolated using Chebyshev polynomials, CODE(0x560026ecf290)ch
is done within RealTimeTransport. We will therefore outline the main steps to
derive a time-space perturbation theory for $\mathcal{K}$ using the
superfermionic formulation as presented in detail in Refs. [Saptsov and
Wegewijs, 2012, 2014]. This formulation is also at the heart of the more
powerful methods presented in the subsequent sections. The crucial starting
point is the definition of superfield operators as

$\displaystyle
J^{p}_{1}\bullet\coloneqq\frac{1}{\sqrt{2}}\,\Big{[}a_{1}\bullet+p(-\mathds{1})^{n_{R}}\bullet(-\mathds{1})^{n_{R}}a_{1}\Big{]},$
(12)

where $(-\mathds{1})^{n_{R}}$ denotes the fermion parity of the reservoirs.
These superfields are called _superfermions_ , because they act analogously to
ordinary field operators $a_{1}$ in Liouville space. For example, they
anticommute as

$\displaystyle\left\\{J_{1}^{p_{1}},J_{2}^{p_{2}}\right\\}=\delta_{p_{1}\bar{p}_{2}}\delta_{1\bar{2}}$
(13)

and obey a super-Pauli principle $(J^{p}_{1})^{2}=0$. CODE(0x560026ecf290)s,
$J^{+}_{1}$ ($J^{-}_{1}$) can be thought of as creating (annihilating) an
excitation in Liouville space. By defining similar superfermions in the
quantum dot systems Liouville space,

$\displaystyle
D^{p}_{1}\bullet\coloneqq\frac{1}{\sqrt{2}}\,\Big{[}d_{1}\bullet+p(-\mathds{1})^{n}\bullet(-\mathds{1})^{n}d_{1}\Big{]},$
(14)

we can express the tunneling Liouvillian
$\mathcal{L}_{T}\coloneqq[H_{T},\bullet]$ as

$\displaystyle\mathcal{L}_{T}=\sqrt{\lambda_{1}}t_{1}D_{1}^{p}J_{\bar{1}}^{\bar{p}}.$
(15)

Here $\bar{p}\coloneqq-p$ and we implicitly sum/integrate over $p$ and
$1=\eta,l,r,\nu,\omega$. Note that in writing Eq. (15) we explicitly used the
fermion parity superselection principle
$(-\mathds{1})^{n_{\text{tot}}}\rho_{\text{tot}}(-\mathds{1})^{n_{\text{tot}}}=\rho_{\text{tot}}$,
where $(-\mathds{1})^{n_{\text{tot}}}$ denotes the fermion parity of the total
system and $\rho_{\text{tot}}$ any physical state. CODE(0x560026ecf290)fining
the unperturbed propagators of the dot and reservoir as $\Pi_{0}(t)\coloneqq
e^{-i\mathcal{L}t}$ and $\Pi_{R}(t)\coloneqq e^{-i\mathcal{L}_{R}t}$
respectively, we can expand the full propagator $\Pi$ in orders of the
tunneling Liouvillian $\mathcal{L}_{T}$. By implicitly summing over
$n\in\mathbbm{N}_{0}$ and integrating over all intermediate times $\tau_{i}$
in a time-ordered fashion $t\geq\tau_{1}\geq\dots\geq\tau_{n}\geq 0$ we obtain

$\displaystyle\Pi(t)=$
$\displaystyle(-i)^{n}\big{\langle}\Pi_{0R}(t-\tau_{1})\mathcal{L}_{T}\dots\mathcal{L}_{T}\Pi_{0R}(\tau_{n})\big{\rangle}_{R}$
(16) $\displaystyle=$
$\displaystyle(-i)^{n}\sqrt{\lambda_{1}}t_{1}\dots\sqrt{\lambda_{n}}t_{n}\Pi_{0}(t-\tau_{1})D_{1}^{p_{1}}\dots
D_{n}^{p_{n}}\Pi_{0}(\tau_{n})$
$\displaystyle\times\big{\langle}\Pi_{R}(t-\tau_{1})J_{\bar{1}}^{\bar{p}_{1}}\Pi_{R}(\tau_{1}-\tau_{2})\dots
J_{\bar{n}}^{\bar{p}_{n}}\Pi_{R}(\tau_{n})\big{\rangle}_{R}.$ (17)

In (16) we used the shorthand notation $\Pi_{0R}=\Pi_{0}\Pi_{R}$ and
$\langle\dots\rangle_{R}\coloneqq\operatorname{Tr}_{R}(\dots)\rho_{R}$. In
(17) we inserted (15) and used that by construction all system superoperators
commute with those of the reservoir, allowing us to pull all of them to the
left. The last line of Eq. (17) only involves reservoir fields and can be
evaluated using a superoperator version of Wick’s theorem,
CODE(0x560026ecf290)ch decomposes the $n$-point expectation value $\langle
J_{1}^{p_{1}}\dots J_{n}^{p_{n}}\rangle_{R}$ into sums of products of two-
point expectation values $\langle J_{1}^{p_{1}}J_{2}^{p_{2}}\rangle_{R}$. A
key simplification achieved by the superfermions (12) is that out of the four
two-point expectation values one obtains by setting $p_{1/2}=\pm$, only two
are non-zero. These are given by

$\displaystyle\langle J_{1}^{-}J_{2}^{+}\rangle_{R}$
$\displaystyle=\delta_{1\bar{2}},$ (18a) $\displaystyle\langle
J_{1}^{-}J_{2}^{-}\rangle_{R}$
$\displaystyle=\delta_{1\bar{2}}\tanh\left(\frac{\eta_{2}\beta_{r_{2}}\omega_{2}}{2}\right).$
(18b)

The internal $\omega_{i}$ integrations within Eq. (17) can then be
incorporated into time-space contraction functions defined by

$\displaystyle\gamma_{\eta_{1}l_{1},\eta_{2}l_{2}}^{p}(\tau)\coloneqq\sqrt{\lambda_{1}\lambda_{2}}t_{1}t_{2}\langle
J_{\bar{1}}^{-}J_{\bar{2}}^{p}\rangle_{R}e^{-i\eta_{1}(\omega_{1}+\mu_{r_{1}})\tau},$
(19)

where we again implicitly sum/integrate over all indices on the right hand
side which do not occur on the left hand side. Evaluating the frequency
integrals we obtain

$\displaystyle\gamma_{\eta_{1}l_{1},\eta_{2}l_{2}}^{+}(\tau)$
$\displaystyle=\delta_{\eta_{1}\bar{\eta}_{2}}\tfrac{1}{2}\sum_{r}\Gamma_{rl_{1}l_{2}}\bar{\delta}(\tau),$
(20a) $\displaystyle\gamma_{\eta_{1}l_{1},\eta_{2}l_{2}}^{-}(\tau)$
$\displaystyle=\delta_{\eta_{1}\bar{\eta}_{2}}\sum_{r}\frac{-i\Gamma_{rl_{1}l_{2}}T_{r}}{\sinh(\pi
T_{r}\tau)}e^{-i\eta_{1}\mu_{r}\tau}.$ (20b)

The $\bar{\delta}$ distribution occurring in $\gamma^{+}$ – defined as
$\int_{0}^{t}\bar{\delta}(t-s)f(s)ds=f(t)$ – stems from the wideband limit
taken in the definition of the tunneling Hamiltonian (9). Thus, we obtain the
series expansion of the propagator in orders of the system-environment
coupling, in a form where all reservoir degrees of freedom have been
integrated out:

$\displaystyle\Pi(t)=$
$\displaystyle(-i)^{n}\sum_{\text{contr.}}(-1)^{N_{p}}\big{(}\prod_{\text{pairs
$i<j$}}\delta_{p_{i}+}\gamma_{ij}^{\bar{p}_{j}}(\tau_{i}-\tau_{j})\big{)}$
$\displaystyle\times\Pi_{0}(t-\tau_{1})D_{1}^{p_{1}}\Pi_{0}(\tau_{1}-\tau_{2})\dots
D_{n}^{p_{n}}\Pi_{0}(\tau_{n}).$ (21)

Each term in Eq. (21) can be represented by a diagram, where vertices
(superfermions), drawn as bullets, interrupt the bare evolutions $\Pi_{0}$ and
are contracted in pairs. The memory kernel $\mathcal{K}$ is then given by the
sum of all connected diagrams without legs, i.e., without the leftmost and
rightmost bare propagator,

$\displaystyle\mathcal{K}(t)=$
$\displaystyle(-i)^{n}\sum_{\text{connected}}(-1)^{N_{p}}\big{(}\prod_{\text{pairs
$i<j$}}\delta_{p_{i}+}\gamma_{ij}^{\bar{p}_{j}}(\tau_{i}-\tau_{j})\big{)}$
$\displaystyle\times
D_{1}^{p_{1}}\Pi_{0}(t-\tau_{1})D_{2}^{p_{2}}\dots\Pi_{0}(\tau_{n-2})D_{n}^{p_{n}},$
(22)

see the next section for examples. The factor $(-1)^{N_{p}}$ refers to the
Wick-contraction sign, which is given by the number of crossings of
contraction functions.

###  II.3 Renormalized perturbation theory around the infinite-temperature
limit

The $\bar{\delta}$-function inside the memory kernel leads to a time-local
contribution,

$\displaystyle\mathcal{K}(t-s)=\Sigma_{\infty}\bar{\delta}(t-s)+\Sigma(t-s),$
(23)

where

$\displaystyle\Sigma_{\infty}\coloneqq-\frac{i}{2}\sum_{\eta
rll^{\prime}}\Gamma_{rll^{\prime}}D^{+}_{\eta l}D^{-}_{\bar{\eta}l^{\prime}}.$
(24)

Since the $\gamma^{-}$ contraction function (20b) vanishes if all
$T_{r}\rightarrow\infty$, we can infer that the Liouvillian

$\displaystyle L_{\infty}\coloneqq L+\Sigma_{\infty}$ (25)

generates the exact semigroup dynamics at infinite temperature:

$\displaystyle\Pi_{\infty}(t)\coloneqq\lim_{T_{r}\rightarrow\infty}\Pi(t)=e^{-iL_{\infty}t}.$
(26)

This infinite temperature solution can be used as a new reference in a
_renormalized_ perturbation series. The key technical step is to realize that
two vertices connected by a $\gamma^{+}$ contraction can never have any other
vertices between them, due to vanishing support of the intermediate time
integrations, which enables the resummation of all $\gamma^{+}$ contractions
into $\Pi_{\infty}$ propagators (see Ref. [Saptsov and Wegewijs, 2014] for
details). It is now easy to obtain the rules of the renormalized perturbation
theory, for which one has to adapt those of the bare one as follows: (i) only
creation superfermions $D^{+}_{1}$ and $\gamma^{-}_{1}$ contractions are
allowed and (ii) every bare Liouvillian $\mathcal{L}$ must be renormalized as
$\mathcal{L}\rightarrow\mathcal{L}_{\infty}$, which includes replacing bare
propagators as $\Pi_{0}\rightarrow\Pi_{\infty}$. RealTimeTransport implements
the leading and next-to-leading order of this renormalized perturbation
theory. In leading order, the renormalized memory kernel is explicitly given
by

$\displaystyle-i\Sigma^{(1)}(t)=\raisebox{-0.3pt}{\includegraphics[height=12.80365pt]{diagrams/memory_kernel_1.pdf}}=-\gamma^{-}_{12}(t)D_{1}^{+}\Pi_{\infty}(t)D_{2}^{+}.$
(27)

It is instructive to analyze the worst case $T\rightarrow 0$, where the bare
perturbation theory suffers from a complete breakdown due to the very slow and
oscillatory decay of the contraction function $\gamma^{-}(t)\propto
e^{-i\eta_{1}\mu_{r}t}/t$. CODE(0x560026ecf290)ontrast, the limit
$T\rightarrow 0$ is unproblematic in Eq. (27) because the intermediate
renormalized propagator is still exponentially decaying towards the infinte-
temperature stationary state
$\Pi_{\infty}(t)\rightarrow\tfrac{1}{d}\bm{|}\mathds{1}\bm{)}\bm{(}\mathds{1}\bm{|}$,
which is a left zero eigenvector of the $D^{+}_{2}$ creation superfermion on
the right. CODE(0x560026ecf290)e next-to-leading order corrections are

$\displaystyle-i\Sigma^{(2)}=\raisebox{-0.3pt}{\includegraphics[height=12.80365pt]{diagrams/memory_kernel_2.pdf}}+\raisebox{-0.3pt}{\includegraphics[height=12.80365pt]{diagrams/memory_kernel_3.pdf}},$
(28)

where the diagrams are given by

$\displaystyle\raisebox{-0.3pt}{\includegraphics[height=12.80365pt]{diagrams/memory_kernel_2.pdf}}=$
$\displaystyle\gamma^{-}_{14}(t)\gamma^{-}_{23}(\tau_{1}-\tau_{2})D^{+}_{1}\Pi_{\infty}(t-\tau_{1})D^{+}_{2}$
$\displaystyle\times\Pi_{\infty}(\tau_{1}-\tau_{2})D^{+}_{3}\Pi_{\infty}(\tau_{2})D^{+}_{4},$
(29)
$\displaystyle\raisebox{-0.3pt}{\includegraphics[height=12.80365pt]{diagrams/memory_kernel_3.pdf}}=$
$\displaystyle-\gamma^{-}_{13}(t-\tau_{2})\gamma^{-}_{24}(\tau_{1})D^{+}_{1}\Pi_{\infty}(t-\tau_{1})D^{+}_{2}$
$\displaystyle\times\Pi_{\infty}(\tau_{1}-\tau_{2})D^{+}_{3}\Pi_{\infty}(\tau_{2})D^{+}_{4}.$
(30)

Here we again integrate in a time-ordered fashion over all intermediate times
$\tau_{i}$, i.e., $t\geq\tau_{1}\geq\tau_{2}\geq 0$, and sum over all occuring
multi-indices. Importantly, the singularities of $\gamma^{-}_{ij}(\tau)$ at
$\tau=0$ within Eqs. (27)–(30) never contribute due to the algebra of the
superfermions, which was shown in detail in App. A of Ref. [Nestmann and
Wegewijs, 2021]. However, they need to be implemented carefully in order to
not cause numerical problems. Finally, we mention that for non-interacting
systems it can be shown that the renormalized perturbation theory terminates
due to the super-Pauli principle. CODE(0x560026ecf290)cifically, for these
systems the renormalized perturbation theory of order $n_{l}$ is exact, where
$n_{l}$ denotes the number of single particle states. For example, this means
that for a single spinful level (without interactions) Eqs. (27)–(30) give the
exact memory kernel.

###  II.4 Renormalization group transformation

We now follow Ref. [Nestmann and Wegewijs, 2022] to derive a renormalization
group transformation with a fixed point given by the time-dependent memory
kernel $\Sigma$ [Eq. (23)]. Firstly, in each diagram we resum all connected
subblocks without uncontracted lines to full propagators $\Pi$, which we draw
as double lines. The remaining diagrams are called irreducible:

$\displaystyle-i\Sigma=\raisebox{-0.3pt}{\includegraphics[height=12.80365pt]{diagrams/memory_kernel_1_resummed.pdf}}+\raisebox{-0.3pt}{\includegraphics[height=12.80365pt]{diagrams/memory_kernel_3_resummed.pdf}}+\cdots.$
(31)

For example, this leads to the diagram (29) being already contained within the
first diagram of Eq. (31). Secondly, we define the effective 1-point vertex
$\bar{D}_{1}$ as sum over all irreducible diagrams with 1 uncontracted line,

$\displaystyle\bar{D}_{1}\equiv\raisebox{-0.3pt}{\includegraphics[height=12.80365pt]{diagrams/effective_vertex_1.pdf}}=\raisebox{-0.3pt}{\includegraphics[height=12.80365pt]{diagrams/bare_vertex_1.pdf}}+\raisebox{-0.3pt}{\includegraphics[height=12.80365pt]{diagrams/bare_vertex_2.pdf}}+\cdots,$
(32)

Because the effective vertex $\bar{D}_{1}$ represents a block of time, we need
to distinguish between the time arguments of the leftmost vertex $t$, the
uncontracted vertex $\tau$, and the rightmost vertex $s$, i.e.,
$\bar{D}_{1}\equiv\bar{D}_{1}(t,\tau,s)$. However, since we are only dealing
with time-translation invariant Hamiltonians, the effective vertex only
depends on time differences, $\bar{D}_{1}\equiv\bar{D}_{1}(t-\tau,\tau-s)$. We
can now utilize the effective vertex to express the memory kernel as
CODE(0x560026ecf290)iΣ(t) &=
= - γ^-_12(t-τ_2) D_1^+ Π(t-τ_1) ¯D_2(τ_1-τ_2,τ_2). This result is formally
exact and represents the memory kernel as a functional of the full propagator
$\Pi$ and the effective vertex $\bar{D}_{1}$. However, since the propagator
itself is determined by the memory kernel $\Sigma$ via Eq. (II.2), we see that
Eq. (II.4) sets up a self-consistent RG transformation of the form

$\displaystyle\Sigma=\mathcal{F}\big{[}\Sigma,\bar{D}_{1}\big{]}.$ (33)

Defining in an analogous way the effective two-point vertex as the sum over
all irreducible diagrams with two free legs,

$\displaystyle\bar{D}_{12}\equiv\raisebox{-0.3pt}{\includegraphics[height=12.80365pt]{diagrams/effective_vertex_2.pdf}}=\raisebox{-0.3pt}{\includegraphics[height=12.80365pt]{diagrams/effective_vertex_2_1.pdf}}+\cdots$
(34)

we find that the effective one-point vertex can be expressed as,
CODE(0x560026ecf290)diagrams/effective_vertex_1.pdf =  \+  \+  \+ , which
represents a self-consistent equation enslaved to Eq. (II.4) of the form

$\displaystyle\bar{D}_{1}=\mathcal{F}_{1}\big{[}\Sigma,\bar{D}_{1},\bar{D}_{12}\big{]}.$
(35)

This indicates that the complete series (31) can be rewritten as an infinite
hierarchy of self-consistent equations for $n$-point effective vertices, where
the $n$-point vertex is determined by the $n+1$-point vertex. This is similar
in spirit to the flow equations of the functional renormalization group for
Green’s functions, where the self-energy is determined through an infinite
hierarchy of $n$-particle vertices. CODE(0x560026ecf290) Ref. [Nestmann and
Wegewijs, 2022] Eqs. (II.4)–(II.4) were approximately solved by transforming
them into a continuous renormalization group flow using the environment
temperature $T$ as flow parameter (dubbed $T$-flow). The main idea was to
compute memory kernel corrections $\delta\Sigma$ which occur when the
temperature is lowered in many small steps $\delta T$, thereby generating
higher order coupling effects. In RealTimeTransport we instead solve these
equations iteratively in a 2- or 3-loop scheme as described below. We
numerically checked that the solution generated this way is equivalent to the
corresponding solution obtained from the continuous $T$-flow, but find that
this discrete fixed-point iteration is more efficient than solving the
differential $T$-flow RG equations. In the 2-loop scheme we set the effective
two point vertex $\bar{D}_{12}$ to zero and approximate the one-point vertex
$\bar{D}_{1}$ via the first two terms of Eq. (32). This gives a self-
consistent scheme for the memory kernel of order $\mathcal{O}({D^{+}}^{4})$.
Here we count the number of bare creation superfermions in each term, which
are proportional to the number of contraction functions $\gamma^{-}$. As
initial approximation for the memory kernel we use the renormalized next-to-
leading order memory kernel [Eqs. (27)–(30)] and solve Eq. (II.2) to obtain
the corresponding propagator. We use these to obtain an initial approximation
of $\bar{D}_{1}$, which we can then use to recompute the memory kernel and
propagator. This scheme is iterated until the error between subsequent
iterations is below some threshold. In the 3-loop scheme we instead
approximate the effective two point vertex $\bar{D}_{12}$ via the first term
of Eq. (34) and expand the right hand side of Eq. (II.4) such that all terms
of order $\mathcal{O}({D^{+}}^{5})$ are contained. This gives a self-
consistent scheme for the memory kernel of order $\mathcal{O}({D^{+}}^{6})$.
The 3-loop iteration is instead started with the two-loop memory kernel. The
idea is that this starting point should be closer to the three loop memory
kernel, thus requiring fewer iterations to achieve convergence.

##  III Library architecture

RealTimeTransport currently implements the leading and next-to-leading order
renormalized perturbation theory as well as the discussed 2- and 3-loop
renormalization group transformations of the previous section. Thus, the
memory kernel has to be resolved in time-space. Internally this is achieved
through a Chebyshev interpolation with an adaptively chosen degree, based on a
user-provided error goal. This approach leverages the fact that the memory
kernel is typically a well-behaved, exponentially decaying function of time.
For such smooth functions Chebyshev interpolations are nearly optimal in a
well-defined sense CODE(0x560026ecf290) vastly superior to interpolations on
equidistant grids. The primary challenge is that the occurring superoperators
are large matrices of dimension $d^{2}\times d^{2}$, where $d$ denotes the
Hilbert-Fock space dimension, which need to be multiplied with each other many
times when diagrams are evaluated by the numerical integration routines. To
reduce the size of matrices which need to be handled we use that the memory
kernel has a block diagonal structure. Specifically, if an observable of the
total system $N_{\text{tot}}=N+N_{R}$ is conserved and furthermore
$[H,N]=[H_{R},N_{R}]=0$, then

$\displaystyle\bm{(}s_{1}s_{2}\bm{|}\Sigma\bm{|}s_{1^{\prime}}s_{2^{\prime}}\bm{)}\neq
0\,\Longrightarrow\,N_{s_{1}}-N_{s_{2}}=N_{s_{1^{\prime}}}-N_{s_{2^{\prime}}}$
(36)

where $\bm{(}A\bm{|}B\bm{)}\coloneqq\operatorname{Tr}A^{\dagger}B$ denotes the
Hilbert-Schmidt scalar product, $\bm{|}s_{1}s_{2}\bm{)}=|s_{1}\rangle\langle
s_{2}|$ and $N|s\rangle=N_{s}|s\rangle$, see also App. F of Ref. [Saptsov and
Wegewijs, 2012]. This means that each block of the memory kernel is
characterized by a difference $\Delta N=N_{s_{1}}-N_{s_{2}}$, provided the
basis vectors are correctly organized. Each block of the memory kernel is thus
internally represented by a matrix-valued Chebyshev interpolation object. We
furthermore use that the superfermions $D^{+}_{1}$ have a sparse structure
inherited from the field operators $d_{1}$, i.e., almost all matrix elements
are zero. This allows us to separate each superfermion into a block structure
defined by the memory kernel and only save the non-zero blocks as dense
matrices. The remaining linear algebra with dense matrices is handled using
the well-established Eigen library. CODE(0x560026ecf290)e memory kernel
computation is implemented flexibly in a model-independent way. To create a
new model one has to define a class that inherits from an abstract Model class
and implement several methods, which provide, for example, the Hamiltonian,
the field operators and information about the block structure of the memory
kernel. Several important models have already been implemented in the library,
such as the single and double quantum dot systems discussed in the next
section.

##  IV Examples

In this section we exemplify the practical usage of the RealTimeTransport
library, showcase each of the four implemented methods and discuss their
scopes of applicability.

###  IV.1 Double quantum dot – first order renormalized perturbation theory

As a first example, we consider a spinless double quantum dot connected to two
biased reservoirs. In the notation of Sec. II.1 the index $l\in\\{1,2\\}$
refers to the dot index and the reservoir index $\nu$ can be suppressed. Thus
$d_{l}^{\dagger}$ ($d_{l}$) creates (annihilates) an electron in the dot $l$.
We take the Hamiltonian of the double dot to be

$\displaystyle
H=\epsilon_{1}n_{1}+\epsilon_{2}n_{2}+Un_{1}n_{2}+\Omega(d_{1}^{\dagger}d_{2}+d_{2}^{\dagger}d_{1}),$
(37)

where $n_{l}\coloneqq d_{l}^{\dagger}d_{l}$ is the occupation of dot $l$ with
energy $\epsilon_{l}$, $U$ the Coulomb interaction and $\Omega$ the
hybridization between the dots. The coupling between each dot and its adjacent
leads can be parameterized by $\Gamma_{rl}\coloneqq
2\pi\lambda_{r}|t_{rl}|^{2}$ such that
$\Gamma_{rll^{\prime}}=\sqrt{\Gamma_{rl}\Gamma_{rl^{\prime}}}$ [Eq. (10)]. In
the following we focus on a serial setup with
$\Gamma_{L1}=\Gamma_{R2}\equiv\Gamma$, $\Gamma_{L2}=\Gamma_{R1}=0$ and
furthermore set $\epsilon_{1}=\epsilon_{2}\equiv\epsilon$. An interesting
feature of this model is that coherences, described by non-diagonal elements
of the density matrix, do not vanish in the stationary state, for which we
illustrate the computation using the first order renormalized perturbation
theory in Listing 1. The complete code for this can be found in the examples
folder of the repository.CODE(0x560026ecf290)ce the double dot model is
already implemented in the library, we only need to select parameters (line 1–
11) and create the model using the createModel function (line 14). The cost of
computing the memory kernel is mainly determined by the maximum time up to
which it should be computed and an accuracy goal for the interpolation of the
memory kernel in time space. After the memory kernel is computed (line 22) the
stationary state can be directly accessed (line 26).

⬇

1// Double dot parameters

2Real e = -1; // Dot energies

3Real U = 5; // Coulomb interaction

4Real Omega = 2; // Hybridisation

5Real V = 0.5; // Bias voltage

6RealVector T{{1, 1}}; // Temperatures leads

7RealVector mu{{V/2, -V/2}}; // Chem. potentials

8

9// Serial setup: lead L - dot 1 - dot 2 - lead R

10RealVector Gamma1{{1, 0}}; // Dot 1 - lead L

11RealVector Gamma2{{0, 1}}; // Dot 2 - lead R

12

13// Create model

14auto model = createModel<DoubleDot>(

15 e, e, U, Omega, T, mu, Gamma1, Gamma2);

16

17// Computation parameters

18Real tMax = 10; // Maximum simulation time

19Real errGoal = 1e-6; // Interpolation error goal

20

21// Memory kernel computation

22auto K = computeMemoryKernel(model,

23 RenormalizedPT::Order::_1, tMax, errGoal);

24

25// Compute stationary state and print it

26Matrix rho = K.stationaryState();

27std::cout << "rho =\n" << rho << "\n";

Listing 1: Computation of the stationary state using the first order
renormalized perturbation theory.

To access the particle currents into the reservoirs one has to compute an
additional current kernel, which is shown in Listing 2. After the current
kernel is computed (line 3) the stationary current can be directly accessed by
passing the previously obtained stationary state as a parameter (line 7).

⬇

1// Current kernel computation

2int r = 0; // Left lead

3auto K_I = computeCurrentKernel(model, r,

4 RenormalizedPT::Order::_1, tMax, errGoal);

5

6// Compute current and print it

7Real I = K_I.stationaryCurrent(rho);

8std::cout << "I = " << I << "\n";

Listing 2: Computation of the stationary current. Figure 1: Leading order
renormalized perturbation theory at $T=2\Gamma$, $\Gamma$ and $0$: Stationary
current as function of interdot hybridization $\Omega$ and dot energies
$\epsilon$ for a double quantum dot with interaction $U=5\Gamma$ at bias
voltage $V=\Gamma/2$.

We illustrate the dependence of the stationary current on the hybridization
$\Omega$ and dot energies $\epsilon$ by sweeping these parameters for three
different temperatures $T\in\\{0,1,2\\}\Gamma$ in Fig. 1 and plot the ratio
$I/I_{\text{max}}$, where $I_{\text{max}}$ denotes the maximum current within
the scanned parameter range. As expected one can clearly see that the
broadening of the transport resonances is decreasing with temperature, but
unlike in bare perturbation theory, the width remains finite for $T\rightarrow
0$. Furthermore, as $\Omega\rightarrow 0$ the current vanishes because no
tunneling between the dots is possible. This basic feature is completely
missed by the popular Pauli master equation, which is even at high
temperatures unable to describe the current correctly for
$\Omega\lesssim\Gamma$ due to the role of coherences. Other commonly used
first order approaches, such as the von Neumann or Redfield approach, which do
account for coherences break down when $T\lesssim\Gamma$, see Ref. [Kiršanskas
_et al._ , 2017] for a detailed discussion. The two blue lines of maximum
current of width $V$ in Fig. 1(c) are due to single electron tunneling (SET)
and mark the transition of the stationary total dot occupation $n=n_{1}+n_{2}$
between $2\leftrightarrow 1$ and $1\leftrightarrow 0$ as indicated. We note
that cotunneling effects are important away from the resonances and not
included consistently through the leading order renormalized perturbation
theory. Here the next-to-leading order renormalized perturbation theory or the
self-consistent RG methods should be utilized, which are illustrated in the
next section.

###  IV.2 Transient Coulomb diamonds – second order renormalized perturbation
theory

We now consider the single impurity Anderson model with energy $\epsilon$,
magnetic field $B$ and Coulomb repulsion $U$ described by the Hamiltonian

$\displaystyle
H=\epsilon\left(n_{\uparrow}+n_{\downarrow}\right)+\tfrac{1}{2}B\left(n_{\uparrow}-n_{\downarrow}\right)+Un_{\uparrow}n_{\downarrow}.$
(38)

The indices $l$ and $\nu$ of Sec. II.1 are spin indices $\sigma$. We take the
tunnel rate $\Gamma$ to be spin- and reservoir-independent, i.e.,
$\Gamma_{r\sigma\sigma^{\prime}}\equiv\Gamma\delta_{\sigma\sigma^{\prime}}$
[Eq. (10)]. The computation of transient occupations using the next-to-leading
order renormalized perturbation theory is shown in Listing 3 and the complete
code can also be found in the examples folder of the repository.
CODE(0x560026ecf290)er the model is set up similarly to before (line 1–8) we
use an additional block argument in the computation of the memory kernel and
propagator (line 18 and 20). This uses the fact that the memory kernel and
propagator are block diagonal, since the particle number and spin component
along the magnetic field are conserved. For the considered model the
stationary state is encoded in a single block, which can also be used to
compute transients as long as the initial state is described by a diagonal
density matrix. Thus, only computing a single block accelerates the
computation. Omitting the block argument as in Listing 1 defaults to computing
the full memory kernel and propagator.

⬇

1// Create Anderson dot with two reservoirs

2Real epsilon = -4; // Dot energy

3Real B = -1; // Magnetic field

4Real U = 10; // Coulomb repulsion

5RealVector T{{0, 0}}; // Temperatures leads

6RealVector mu{{2, -2}}; // Chem. potentials

7RealVector Gamma{{1, 1}}; // Tunnel rate to leads

8auto model = createModel<AndersonDot>(

9 epsilon, B, U, T, mu, Gamma);

10

11// Computation parameters

12Real tMax = 5; // Maximum simulation time

13Real errGoal = 1e-4; // Interpolation error goal

14int block = 0; // Only compute first block

15

16// Compute memory kernel & propagator

17auto method = RenormalizedPT::Order::_2;

18auto K = computeMemoryKernel(

19 model, method, tMax, errGoal, block);

20auto Pi = computePropagator(K, block);

21

22// Set initial state: Basis 0, Up, Down, UpDown

23Matrix rho0 = Matrix::Zero(4, 4);

24rho0.diagonal() << 0, 1, 0, 0;

25

26// Operator nUp

27Matrix nUp = Matrix::Zero(4, 4);

28nUp.diagonal() << 0, 1, 0, 1;

29

30// Compute occupations Tr{nUp*rho(t)}

31auto t = RealVector::LinSpaced(100, 0, tMax);

32RealVector n(t.size());

33for (int i = 0; i < t.size(); ++i)

34 n[i] = (Pi(t[i], rho0) * nUp).trace().real();

35

36// Print

37std::cout << "t = " << t.transpose() << "\n"

38 << "n = " << n.transpose() << "\n";

Listing 3: Computation of transient occupations $\langle
n_{\uparrow}\rangle(t)$.

The computation of transient currents is shown in Listing 4. The current
kernel is computed (line 10) and used together with the previously obtained
propagator in the function computeCurrent (line 14), which returns an object
containing a Chebyshev interpolation of the transient current. This allows to
efficiently evaluate the current at arbitrary (vectors of) time arguments
(line 15).

⬇

8// Compute current kernel

9int r = 0; // Left reservoir

10auto K_I = computeCurrentKernel(

11 model, r, method, tMax, errGoal, block);

12

13// Compute transient current and print

14auto I = computeCurrent(K_I, Pi, rho0);

15RealVector I_t = I(t);

16std::cout << "I = " << I_t.transpose() << "\n";

Listing 4: Computation of the transient currents. Figure 2: Next-to-leading
order renormalized perturbation theory at $T=0$: Transient occupation of the
excited state $\langle n_{\uparrow}\rangle(t)$ as function of bias voltage
$V=2\mu_{L}=-2\mu_{R}$ and gate voltage $V_{g}=-\epsilon$ for an initially
empty Anderson quantum dot $\rho_{0}=|0\rangle\\!\langle 0|$ with strong
interaction $U=30\Gamma$ and magnetic field $B=5\Gamma$ .

In Fig. 2 we show the transient occupation $\langle n_{\uparrow}\rangle(t)$
for an initially empty dot, $\rho_{0}=|0\rangle\\!\langle 0|$, as function of
the bias voltage $V$ and gate voltage $V_{g}=-\epsilon$ for several times $t$.
These results can be rationalized in detail based on well-established
understanding of the stationary limit based on the transition rates.
CODE(0x560026ecf290) Fig. 2(a)–(b) fast SET processes fill the dot on a time
scale $t\Gamma\sim 1$ to a total number of $n\coloneqq
n_{\uparrow}+n_{\downarrow}\approx 1$ and $2$ electrons as indicated in (b)
(except in the $n\approx 0$ region where the dot is already stationary). In
these regions, the rate of filling is essentially independent of the spin
value, allowing the excited level $|\\!\\!\uparrow\rangle$ to be rapidly
occupied with the same probability as the ground level
$|\\!\\!\downarrow\rangle$. In the $n\approx 1$ Coulomb diamond [white center
in Fig. 2(b)] this value $n_{\uparrow}\approx 1/2$ exceeds the stationary
probability [Fig. 2(f)] and multi-stage decay sets in. At first this decay
only occurs along a band of width $\propto B$ _inside_ the Coulomb diamond
[dark blue in Fig. 2(c)]. Here the excess energy $B$ of the excited state
$|\\!\\!\uparrow\rangle$ enables SET relaxation indirectly via the empty state
$|0\rangle$ to the ground state $|\\!\\!\downarrow\rangle$ at rate
$\propto\Gamma$. We point out that this band is known as the cotunneling-
assisted SET (CO-SET) regime in _stationary_ transport, since there the
excited level can only be populated by inelastic cotunneling, see Refs.
[Leijnse and Wegewijs, 2008; Gaudenzi  _et al._ , 2017] and references
therein. Here, on the other hand, this band is instead populated by fast SET
into the empty dot due to our choice of initial state. Indeed, at this point
in time ($t\Gamma\sim 1$) cotunneling has had no time yet to make an impact on
the _time-dependent_ transport, as its rate scales as $\propto\Gamma^{2}/U$.
This impact is seen only in Fig. 2(d)–(e) where the excitation probability has
been reduced uniformly throughout the $n\approx 1$ diamond where SET is
energetically forbidden. Here
$|\\!\\!\uparrow\rangle\to|\\!\\!\downarrow\rangle$ relaxation occurs directly
by inelastic cotunneling processes involving _both_ reservoirs and releasing
the excess energy $B$. In Fig. 2(f) this reduction has continued at low bias
$|V|<B$ whereas it stalls at bias $|V|>B$ in the two visible triangles: only
in these regimes _excitation_
$|\\!\\!\downarrow\rangle\to|\\!\\!\uparrow\rangle$ by inelastic cotunneling
is energetically possible while relaxation by CO-SET is forbidden. This
cotunneling re-excitation of the dot effectively counteracts the decay and
establishes a sizeable stationary value in competition with cotunneling
relaxation. CODE(0x560026ecf290)ortantly, RealTimeTransport allows this
rationalization to be further quantified by making use of the time-dependent
memory-kernel matrix accessible internally. For example, one interesting
avenue would be to construct an effective time-dependent rate picture by
converting this memory kernel to a time-local generator. CODE(0x560026ecf290)

###  IV.3 2-loop RG corrections

We now investigate when it is important to complement the renormalized
perturbation theory with the self-consistency imposed by the 2-loop
renormalization group transformation. We plot in Fig. 3 the (stationary)
nonlinear conductance $dI/dV$ predicted by both methods across the lines
indicated in Fig. 2(f). The conductance is computed using a central difference
$dI/dV\approx(I(V+\delta V/2)-I(V-\delta V/2))/\delta V$ with $\delta
V=0.01\Gamma$. Since we have a substantial magnetic field $B$ the Kondo
resonance is not strongly developed here, even at the particle hole symmetric
point $\epsilon=-U/2$ (blue curve). Consistent with this, the RG corrections
are moderate and occur where they are expected: namely inside the Coulomb
diamond they suppress (enhance) the conductance below (above) the onset of ICT
at $V\approx B=5\Gamma$ (blue curve). This has the effect of sharpening the
ICT onset of the conductance relative to the renormalized perturbation theory.
Likewise, the SET signature of the same spin-flip resonance is sharpened. This
is in line with the general spirit of our approach as outlined in the
introduction: We first introduce too much dissipation to subsequently compute
corrections which reduce it in a non-trivial way.

Figure 3:  Nonlinear conductance $dI/dV$ at $T=0$ for three gate voltages
$V_{g}=-\epsilon$ along the indicated cuts in Fig. 2(f) as predicted by the
next-to-leading order renormalized perturbation theory (full lines) and 2-loop
renormalization group (crosses) for the same parameters as in Fig. 2.

From a user perspective, performing a calculation with the renormalization
group methods is very similar to the perturbation theory methods and is shown
in Listing 5. Because the computations are more time consuming it can be
useful to parallelize them. We implemented multi-threading capabilities using
the Taskflow library. CODE(0x560026ecf290)s allows to define an Executor with
a given number of threads (line 3) that can be passed as an additional
argument to compute the memory and current kernels (lines 5 and 10).

⬇

1// Set up model and define parameter as before

2// Compute memory kernel with 8 threads

3tf::Executor executor(8);

4auto method = IteratedRG::Order::_2;

5auto K = computeMemoryKernel(model,

6 method, tMax, errGoal, executor);

7auto Pi = computePropagator(K);

8

9// Compute one block of the current kernel

10auto K_I = computeCurrentKernel(K, Pi, r,

11 method, tMax, errGoal, executor, block);

Listing 5: Multi-threaded computation of the memory and current kernel using
the RG methods. Figure 4:  Comparison of the stationary current at $T=0$ as
function of the bias $V$ at the symmetry point $\epsilon=-U/2$ and magnetic
field $B=0$: next-to-leading order renormalized perturbation theory, 2- and
3-loop RG, tDMRG and the tQMC method.

###  IV.4 3-loop RG corrections

Finally, we turn to the 3-loop renormalization group corrections at $T=0$. We
focus on the particle-hole symmetric point $\epsilon=-U/2$ and the most
challenging case of zero magnetic field, $B=0$. In this case the Kondo effect
becomes important for small bias $V$, which is challenging to describe for
many theoretical methods. CODE(0x560026ecf290)fig:method-comparison-2 we
compare the stationary currents obtained from our next-to-leading order
renormalized perturbation theory and our 2- and 3-loop renormalization group
with the time-dependent density matrix renormalization group (tDMRG) and the
real-time quantum Monte Carlo method (tQMC). The data for the latter two were
taken from Refs. [Werner, Oka, and Millis, 2009; Heidrich-Meisner, Feiguin,
and Dagotto, 2009; Werner  _et al._ , 2010], but note also Ref. [Eckel  _et
al._ , 2010] for further comparisons with the functional renormalization group
and the iterative real-time path integral approach. For small interaction
$U=2\Gamma$ we find that all methods are largely consistent with each other.
The renormalized perturbation theory consistently predicts a current that is
too small at low bias $V$, especially for larger interactions shown up to
$U=10\Gamma$. This can still be seen in our 2-loop RG, however to a much
smaller extent. The 3-loop corrections only become important if $V$ is very
small and $U$ is large. Even there we find good overall agreement between the
currents predicted by our 3-loop RG and the tQMC method. We note that our
discrete RG iteration reproduces the data from the continuous $T$-flow,
CODE(0x560026ecf290)firming the expectation that they converge to the same
fixed point.

##  V Summary

We have presented the open source C++ RealTimeTransport library which is
available at https://github.com/kn-code/RealTimeTransport to model stationary
and transient quantum transport phenomena. The library currently implements
four methods: the leading and next-to-leading order renormalized perturbation
theory around the infinite temperature limit as well as more advanced 2- and
3-loop renormalization group transformations. We exemplified the practical
usage of each of these methods and compared their ranges of applicability. In
contrast to the ordinary bare perturbation theory in the tunnel coupling, our
renormalized perturbation theory can also be applied all the way down to zero
temperature. We have shown that it captures sequential and cotunneling effects
in leading and next-to-leading order, respectively. However, especially at low
temperatures even higher order effects become important. These can be
systematically incorporated using our self-consistent perturbative
renormalization group transformation of the memory kernel [Eq. (33)] of
increasing loop orders. For an Anderson dot we found that the most advanced
3-loop corrections contribute noticeably only if the system is close to the
Kondo resonance. In this case we found overall good agreement between our
3-loop iteration and the tQMC method at finite bias. The version of the
library described in this paper is RealTimeTransport 1.0.

## Acknowledgments

K.N. acknowledges funding from the European Union under the Horizon Europe’s
Marie Skłodowska-Curie Project 101104590. M.L. acknowledges funding from
NanoLund, the Swedish Research Council (Grant Agreement No. 2020-03412) and
the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme under the Grant Agreement No. 856526.

## Data Availability

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

## References

## References

  * Cronenwett, Oosterkamp, and Kouwenhoven (1998) S. M. Cronenwett, T. H. Oosterkamp, and L. P. Kouwenhoven, A Tunable Kondo Effect in Quantum Dots, Science 281, 540–544 (1998).
  * Sasaki  _et al._ (2000) S. Sasaki, S. De Franceschi, J. M. Elzerman, W. G. van der Wiel, M. Eto, S. Tarucha, and L. P. Kouwenhoven, Kondo effect in an integer-spin quantum dot, Nature 405, 764–767 (2000).
  * van der Wiel  _et al._ (2002) W. G. van der Wiel, S. De Franceschi, J. M. Elzerman, T. Fujisawa, S. Tarucha, and L. P. Kouwenhoven, Electron transport through double quantum dots, Rev. Mod. Phys. 75, 1–22 (2002).
  * Su  _et al._ (2016) T. A. Su, M. Neupane, M. L. Steigerwald, L. Venkataraman, and C. Nuckolls, Chemical principles of single-molecule electronics, Nat. Rev. Mater. 1, 1–15 (2016).
  * Gaudenzi  _et al._ (2017) R. Gaudenzi, M. Misiorny, E. Burzurí, M. R. Wegewijs, and H. S. J. van der Zant, Transport mirages in single-molecule devices, J. Chem. Phys. 146 (2017), 10.1063/1.4975767.
  * Josefsson  _et al._ (2018) M. Josefsson, A. Svilans, A. M. Burke, E. A. Hoffmann, S. Fahlvik, C. Thelander, M. Leijnse, and H. Linke, A quantum-dot heat engine operating close to the thermodynamic efficiency limits, Nat. Nanotechnol. 13, 920–924 (2018).
  * Rammer (2007) J. Rammer,  _Quantum Field Theory of Non-equilibrium States_ (Cambridge University Press, Cambridge, England, UK, 2007).
  * Metzner  _et al._ (2012) W. Metzner, M. Salmhofer, C. Honerkamp, V. Meden, and K. Schönhammer, Functional renormalization group approach to correlated fermion systems, Rev. Mod. Phys. 84, 299–352 (2012).
  * Schoeller and Schön (1994) H. Schoeller and G. Schön, Mesoscopic quantum transport: Resonant tunneling in the presence of a strong coulomb interaction, Phys. Rev. B 50, 18436–18452 (1994).
  * Gurvitz and Prager (1996) S. A. Gurvitz and Y. S. Prager, Microscopic derivation of rate equations for quantum transport, Phys. Rev. B 53, 15932–15943 (1996).
  * König, Schoeller, and Schön (1997) J. König, H. Schoeller, and G. Schön, Cotunneling at resonance for the single-electron transistor, Phys. Rev. Lett. 78, 4482–4485 (1997).
  * Timm (2008) C. Timm, Tunneling through molecules and quantum dots: Master-equation approaches, Phys. Rev. B 77, 195416 (2008).
  * Koller  _et al._ (2010) S. Koller, M. Grifoni, M. Leijnse, and M. R. Wegewijs, Density-operator approaches to transport through interacting quantum dots: Simplifications in fourth-order perturbation theory, Phys. Rev. B 82, 235307 (2010).
  * Schoeller (2009) H. Schoeller, A perturbative nonequilibrium renormalization group method for dissipative quantum mechanics, Eur. Phys. Journ. Special Topics 168, 179–266 (2009).
  * Schoeller (2018) H. Schoeller, Dynamics of open quantum systems, (2018), arXiv:1802.10014 [cond-mat.stat-mech] .
  * Anders (2008) F. B. Anders, Steady-state currents through nanodevices: A scattering-states numerical renormalization-group approach to open quantum systems, Phys. Rev. Lett. 101, 066804 (2008).
  * Nghiem and Costi (2014) H. T. M. Nghiem and T. A. Costi, Generalization of the time-dependent numerical renormalization group method to finite temperatures and general pulses, Phys. Rev. B 89, 075118 (2014).
  * Daley  _et al._ (2004) A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, Time-dependent density-matrix renormalization-group using adaptive effective Hilbertspaces, J. Stat. Mech.: Theory Exp. 2004, P04005 (2004).
  * White and Feiguin (2004) S. R. White and A. E. Feiguin, Real-time evolution using the density matrix renormalization group, Phys. Rev. Lett. 93, 076401 (2004).
  * Schollwöck (2005) U. Schollwöck, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259–315 (2005).
  * Rubtsov, Savkin, and Lichtenstein (2005) A. N. Rubtsov, V. V. Savkin, and A. I. Lichtenstein, Continuous-time quantum monte carlo method for fermions, Phys. Rev. B 72, 035122 (2005).
  * Mühlbacher and Rabani (2008) L. Mühlbacher and E. Rabani, Real-time path integral approach to nonequilibrium many-body quantum systems, Phys. Rev. Lett. 100, 176403 (2008).
  * Werner, Oka, and Millis (2009) P. Werner, T. Oka, and A. J. Millis, Diagrammatic monte carlo simulation of nonequilibrium systems, Phys. Rev. B 79, 035320 (2009).
  * Cohen  _et al._ (2015) G. Cohen, E. Gull, D. R. Reichman, and A. J. Millis, Taming the dynamical sign problem in real-time evolution of quantum many-body problems, Phys. Rev. Lett. 115, 266802 (2015).
  * Note (1) The library is available at https://github.com/kn-code/RealTimeTransport.
  * Pletyukhov and Schoeller (2012) M. Pletyukhov and H. Schoeller, Nonequilibrium kondo model: Crossover from weak to strong coupling, Phys. Rev. Lett. 108, 260601 (2012).
  * Bruch  _et al._ (2022) V. Bruch, M. Pletyukhov, H. Schoeller, and D. M. Kennes, Floquet renormalization group approach to the periodically driven kondo model, Phys. Rev. B 106, 115440 (2022).
  * Trefethen (2019) L. N. Trefethen,  _Approximation theory and approximation practice_ (SIAM, 2019).
  * Nestmann and Wegewijs (2022) K. Nestmann and M. R. Wegewijs, Renormalization group for open quantum systems using environment temperature as flow parameter, SciPost Phys. 12, 121 (2022).
  * Johansson, Nation, and Nori (2012) J. R. Johansson, P. D. Nation, and F. Nori, QuTiP: An open-source Python framework for the dynamics of open quantum systems, Comput. Phys. Commun. 183, 1760–1772 (2012).
  * Johansson, Nation, and Nori (2013) J. R. Johansson, P. D. Nation, and F. Nori, QuTiP 2: A Python framework for the dynamics of open quantum systems, Comput. Phys. Commun. 184, 1234–1240 (2013).
  * Kiršanskas  _et al._ (2017) G. Kiršanskas, J. N. Pedersen, O. Karlström, M. Leijnse, and A. Wacker, Qmeq 1.0: An open-source python package for calculations of transport through quantum dot devices, Computer Physics Communications 221, 317–342 (2017).
  * Bockrath  _et al._ (1997) M. Bockrath, D. H. Cobden, P. L. McEuen, N. G. Chopra, A. Zettl, A. Thess, and R. E. Smalley, Single-Electron Transport in Ropes of Carbon Nanotubes, Science 275, 1922–1925 (1997).
  * Nygård, Cobden, and Lindelof (2000) J. Nygård, D. H. Cobden, and P. E. Lindelof, Kondo physics in carbon nanotubes, Nature 408, 342–346 (2000).
  * De Franceschi  _et al._ (2001) S. De Franceschi, S. Sasaki, J. M. Elzerman, W. G. van der Wiel, S. Tarucha, and L. P. Kouwenhoven, Electron cotunneling in a semiconductor quantum dot, Phys. Rev. Lett. 86, 878–881 (2001).
  * Schleser  _et al._ (2005) R. Schleser, T. Ihn, E. Ruh, K. Ensslin, M. Tews, D. Pfannkuche, D. C. Driscoll, and A. C. Gossard, Cotunneling-mediated transport through excited states in the coulomb-blockade regime, Phys. Rev. Lett. 94, 206805 (2005).
  * Saptsov and Wegewijs (2014) R. B. Saptsov and M. R. Wegewijs, Time-dependent quantum transport: Causal superfermions, exact fermion-parity protected decay modes, and pauli exclusion principle for mixed quantum states, Phys. Rev. B 90, 045407 (2014).
  * Wick, Wightman, and Wigner (1952) G. C. Wick, A. S. Wightman, and E. P. Wigner, The intrinsic parity of elementary particles, Phys. Rev. 88, 101–105 (1952).
  * Aharonov and Susskind (1967) Y. Aharonov and L. Susskind, Charge superselection rule, Phys. Rev. 155, 1428–1431 (1967).
  * Nakajima (1958) S. Nakajima, On quantum theory of transport phenomena, Prog. Theor. Phys. 20, 948 (1958).
  * Zwanzig (1960) R. Zwanzig, Ensemble method in the theory of irreversibility, J. Chem. Phys. 33, 1338 (1960).
  * Leijnse and Wegewijs (2008) M. Leijnse and M. R. Wegewijs, Kinetic equations for transport through single-molecule transistors, Phys. Rev. B 78, 235424 (2008).
  * Saptsov and Wegewijs (2012) R. B. Saptsov and M. R. Wegewijs, Fermionic superoperators for zero-temperature nonlinear transport: Real-time perturbation theory and renormalization group for anderson quantum dots, Phys. Rev. B 86, 235432 (2012).
  * Nestmann and Wegewijs (2021) K. Nestmann and M. R. Wegewijs, General connection between time-local and time-nonlocal perturbation expansions, Phys. Rev. B 104, 155407 (2021).
  * Guennebaud, Jacob  _et al._ (2010) G. Guennebaud, B. Jacob, _et al._ , Eigen v3, (2010).
  * Wegewijs and Nazarov (2001) M. R. Wegewijs and Yu. V. Nazarov, Inelastic co-tunneling through an excited state of a quantum dot, arXiv (2001), 10.48550/arXiv.cond-mat/0103579, cond-mat/0103579 .
  * Megier, Smirne, and Vacchini (2020) N. Megier, A. Smirne, and B. Vacchini, The interplay between local and non-local master equations: exact and approximated dynamics, New J. Phys. 22, 083011 (2020).
  * Nestmann, Bruch, and Wegewijs (2021) K. Nestmann, V. Bruch, and M. R. Wegewijs, How quantum evolution with memory is generated in a time-local way, Phys. Rev. X 11, 021041 (2021).
  * Huang  _et al._ (2022) T.-W. Huang, D.-L. Lin, C.-X. Lin, and Y. Lin, Taskflow: A lightweight parallel and heterogeneous task graph computing system, IEEE Trans. Parallel Distrib. Syst. 33, 1303–1320 (2022).
  * Eckel  _et al._ (2010) J. Eckel, F. Heidrich-Meisner, S. G. Jakobs, M. Thorwart, M. Pletyukhov, and R. Egger, Comparative study of theoretical methods for non-equilibrium quantum transport, New Journal of Physics 12, 043042 (2010).
  * Heidrich-Meisner, Feiguin, and Dagotto (2009) F. Heidrich-Meisner, A. E. Feiguin, and E. Dagotto, Real-time simulations of nonequilibrium transport in the single-impurity anderson model, Phys. Rev. B 79, 235336 (2009).
  * Werner  _et al._ (2010) P. Werner, T. Oka, M. Eckstein, and A. J. Millis, Weak-coupling quantum monte carlo calculations on the keldysh contour: Theory and application to the current-voltage characteristics of the anderson model, Phys. Rev. B 81, 035108 (2010).


# Positivity of oscillatory integrals
and Hankel transforms

Yong-Kum Cho111corresponding author<EMAIL_ADDRESS>Department of
Mathematics, College of Natural Sciences, Chung-Ang University, 84 Heukseok-
Ro, Dongjak-Gu, Seoul 06974, Korea. Seok-Young
<EMAIL_ADDRESS>Department of Mathematics, University of
Central Florida, 4393 Andromeda Loop N., Orlando, FL 32816, USA.
Young Woong<EMAIL_ADDRESS>Department of Mathematics, Chung-Ang
University.

Abstract. In consideration of the integral transform whose kernel arises as an
oscillatory solution of certain second-order linear differential equation, its
positivity is investigated on the basis of Sturm’s theory. As applications,
positivity criteria are obtained for Hankel transforms as well as
trigonometric integrals defined on the positive real line.

Keywords. Bessel functions; Fourier transforms; Hankel transforms;

oscillatory; positivity; Sturm’s oscillation theory.

AMS Subject Classifications. 34C10; 42A38; 44A20; 33C10.

##  1 Introduction

In 1823, C. Sturm [1] founded an oscillation theory for the second-order
linear differential equation expressible in the normal form

$u^{\prime\prime}+\phi(t)u=0,\quad t>0.$ (1.1)

In this paper we consider the integral transform of the type

$(Uf)(x)=\int_{0}^{\infty}f(t)u(xt)\,dt\qquad(x>0),$ (1.2)

where $f(t)$ is a nonnegative function for $\,t>0\,$ such that the integral
exists and the kernel $u(t)$ is an oscillatory solution of the differential
equation (1.1), and aim to investigate if the inequality $\,(Uf)(x)>0\,$ holds
true for all $\,x>0\,$ on the basis of Sturm’s oscillation theory.

A great deal of functions encountered in mathematical physics, which are
oscillatory in nature, arise as solutions of (1.1) or images of (1.2). In the
special case $\,\phi(t)=1,\,t>0,$ for instance, (1.2) gives rise to Fourier
sine or cosine transforms on the positive real line. Concerning positivity, a
classical theorem due to Pólya [2] states that if $f(t)$ is a continuous,
positive and decreasing function such that $\,\int_{0}^{\infty}f(t)dt\,$
exists, then

$\displaystyle\left\\{\begin{aligned}
&{\,{\rm(i)}\,\quad\int_{0}^{\infty}f(t)\sin(xt)\,dt>0,}\\\
&{{\rm(ii)}\quad\int_{0}^{\infty}f(t)\cos(xt)\,dt>0,}\end{aligned}\right.$
(1.3)

for all $\,x>0,$ where (ii) holds true under the additional assumption that
$f(t)$ is convex for $\,t>0\,$ (see Tuck [3] for related issues).

Another special case is the Hankel transform defined by

$(H_{\nu}f)(x)=\int_{0}^{\infty}f(t)J_{\nu}(xt)\sqrt{xt}\,dt\qquad(x>0),$
(1.4)

where $J_{\nu}(z)$ stands for the Bessel function of the first kind and order
$\nu$ is assumed to be real with $\,\nu>-1.$ Owing to identities [4, 3.4],

$\displaystyle J_{-1/2}(z)=\sqrt{\frac{2}{\pi z}\,}\cos
z,\,\,J_{1/2}(z)=\sqrt{\frac{2}{\pi z}\,}\sin z,$

it reduces to Fourier transforms of type (1.3) when $\,\nu=\pm 1/2.$

As it will become clearer in the sequel, we deduce positivity of $(Uf)(x)$ by
decomposing it into sums of integrals over intervals between consecutive zeros
of $u(t)$ and by comparing arches of graph of integrands. An advantage of this
approach is that sufficient conditions on $f(t)$ for positivity are mild and
analogous to aforementioned Pólya’s conditions. Our method, however, requires
that $\,u(0+)=0\,$ and zeros of $u(t)$ must be enumerable in ascending order
of magnitude, which could be a drawback in practice.

Concerning positivity of Hankel transforms, we remark that Zhang [5] has
recently obtained a set of positivity criteria applicable to certain class of
completely monotone functions. Unlike our method, his approaches are based on
Laplace transforms of Bessel functions.

The positivity arises often as a critical issue in various disguises. For
example, it is known that the radial function $\,F({\bf x})=f(\|{\bf
x}\|),\,{\bf x}\in\mathbb{R}^{d},$ is positive definite on $\mathbb{R}^{d}$ if
its Fourier transform is nonnegative. In view of the well-known Fourier
transform formula

$\displaystyle\widehat{F}(\xi)=(2\pi)^{\frac{d}{2}}\|\xi\|^{1-\frac{d}{2}}\int_{0}^{\infty}t^{\frac{d}{2}}f(t)J_{\frac{d-2}{2}}(\|\xi\|t)\,dt\qquad(\xi\in\mathbb{R}^{d})$

and the fact that any Hankel transform of order $\,\nu>\mu>-1\,$ can be
written as the Hankel transform of order $\mu$, it is positivity of Hankel
transforms of order $\,\nu\geq(d-2)/2\,$ that plays a decisive role in
determining the positive definiteness of radial functions on $\mathbb{R}^{d}$
(see [6, 7] for details).

The present paper is organized as follows. We revisit Sturm’s oscillation
theory in section 2 with our focuses on the distribution of zeros and the wave
of arches of oscillatory solutions. We establish Pólya-type positivity
criteria for the oscillatory integral (1.2) in section 3, one of which will be
illustrated with Fourier transforms in section 4. We present main applications
to Hankel transforms in section 5 and discuss $Y$-transforms in the last
section.

##  2 Sturm’s oscillation theory

This section aims to revisit some of Sturm’s oscillation theory, original and
extended, which are relevant to the present work. As standard, given an
interval $I$, a function is said to be _oscillatory_ on $I$ when it changes
sign infinitely often. If the sign change occurs at most finitely many times,
it is said to be _non-oscillatory_ on $I$. In what follows, it will be assumed
that the underlying function $\phi(t)$ in (1.1) is continuous for $\,t>0.$

###  2.1 Oscillations with scattered zeros

One of classical Sturm’s theorems states that if there exist $\,a>0\,$ and
$\,m>0\,$ such that $\,\phi(t)\geq m\,$ for all $\,t\geq a,$ then any non-
trivial solution of (1.1) is oscillatory on the interval $[a,\infty)$. As an
alternative, it is Kneser [8], Fite [9] and Hille [10] who proposed the
following weaker conditions:

$\displaystyle{\rm(A1)}\quad\liminf_{t\,\to\,\infty}\big{[}t^{2}\phi(t)\big{]}>\frac{1}{4}.$
$\displaystyle{\rm(A2)}\quad\int_{a}^{\infty}\phi(t)dt=+\infty.$
$\displaystyle{\rm(A3)}\quad\liminf_{t\,\to\,\infty}\left[t\int_{t}^{\infty}\phi(s)ds\right]>\frac{1}{4}.$

######  Lemma 2.1.

Suppose that $\,\phi(t)\geq 0\,$ for all $\,t\geq a,$ where $\,a>0,$ and one
of assumptions (A1), (A2), (A3) holds true. Then any non-trivial solution
$u(t)$ of (1.1) is oscillatory on $[a,\infty).$ Moreover, $u(t)$ has an
infinity of zeros which are all simple, not concentrated on a finite
subinterval of $[a,\infty)$ but scattered over $[a,\infty)$. In particular,
the zeros of $u(t)$ are enumerable in ascending order of magnitude on
$[a,\infty).$

The proof for oscillation is simple and available in many textbooks or
expository articles (see [11, 12], for instance). The simplicity of each zero
is a consequence of the existence and uniqueness theorem of Picard. If the
zeros of $u(t)$ are concentrated on a finite subinterval, then it follows by
the Bolzano-Weierstrass theorem that the zeros have an accumulation point at
which both $u(t),u^{\prime}(t)$ vanish. By Picard’s theorem, $\,u(t)\equiv 0,$
which leads to a contradiction. It is thus proved that the zeros are not
concentrated on a finite interval but scattered over the interval
$[a,\infty)$.

######  Remark 2.1.

If $u(t)$ is oscillatory on $(0,b],\,b>0,$ then the zeros of $u(t)$ has an
accumulation point $0$ and hence may not be enumerable in ascending order of
magnitude. In order to avoid such an occasion, we recall that Kneser and Hille
also proved that if $\,\phi(t)\geq 0\,$ for all $\,t\geq a>0\,$ and

$\limsup_{t\,\to\,\infty}\big{[}t^{2}\phi(t)\big{]}<\frac{1}{4}\quad\text{or}\quad\limsup_{t\,\to\,\infty}\left[t\int_{t}^{\infty}\phi(s)ds\right]<\frac{1}{4},$

then any non-trivial solution of (1.1) is non-oscillatory on $[a,\infty).$

By considering $\,v(t)=tu(1/t)\,$ which obviously satisfies

$v^{\prime\prime}+\frac{1}{t^{4}}\phi\left(\frac{1}{t}\right)v=0,\quad t>0,$

it is simple to infer from the above non-oscillation criteria that any non-
trivial solution of (1.1) is non-oscillatory on the interval $\,(0,b],\,b>0,$
when $\phi(t)$ satisfies either of the following.

$\displaystyle{\rm(B1)}\quad\limsup_{t\,\to\,0+}\big{[}t^{2}\phi(t)\big{]}<\frac{1}{4}.$
$\displaystyle{\rm(B2)}\quad\limsup_{t\,\to\,0+}\left[\frac{1}{t}\int_{0}^{t}s^{2}\phi(s)\,ds\right]<\frac{1}{4}.$

###  2.2 Oscillation magnitudes

As for the distribution of zeros and _relative magnitudes of oscillation_ ,
Sturm’s comparison theorems state as follows (see [13, 14, 1, 4]).

Let $\big{(}v(t),w(t)\big{)}$ be a $C^{1}$ solution of the system of equations

$\left\\{\begin{aligned} &{\,v^{\prime\prime}+\phi(t)v=0,}\\\
&{\,w^{\prime\prime}+\psi(t)w=0,}\end{aligned}\right.\quad t\in I,$

where $I$ is an open interval and $\phi(t),\psi(t)$ are continuous for $\,t\in
I.$

  * (I)

If $\,\phi(t)>\psi(t)\,$ for $\,t\in I,$ then $v(t)$ has at least one zero
between two consecutive zeros of $w(t)$. As a consequence, if
$\,0<m<\phi(t)<M\,$ for $\,t\in I\,$ and $t_{1},t_{2}$ are consecutive zeros
of $v(t)$ with $\,t_{1}<t_{2},$ then

$\frac{\pi}{\,\sqrt{M}\,}<t_{2}-t_{1}<\frac{\pi}{\,\sqrt{m}\,}.$

  * (II)

For $\,a\in I,$ suppose that $\,v(a)=w(a)=0\,$ and
$\,v^{\prime}(a+)=w^{\prime}(a+)>0.$ If $\,\phi(t)<\psi(t)\,$ for
$\,t>a,\,t\in I,$ then $\,v(t)>w(t)\,$ for all $t$ between $a$ and the first
zero of $w(t)$ and hence the first zero of $w(t)$ for $\,t>a\,$ is on the left
of the first zero of $v(t)$ for $\,t>a.$

######  Lemma 2.2.

Suppose that $u(t)$ is a nontrivial solution of (1.1) satisfying $\,u(0+)=0$,
where $\phi(t)$ is continuous, strictly increasing and subject to one of
assumptions (A1), (A2), (A3). Then the following holds true.

  * (i)

$u(t)$ is oscillatory on $(0,\infty)$ with an infinity of simple zeros which
are enumerable in ascending order of magnitude.

  * (ii)

The wave of arches of graph $\,s=u(t)\,$ decreases in the sense

$\left\\{\begin{aligned}
&{\quad\zeta_{k}-\zeta_{k-1}>\zeta_{k+1}-\zeta_{k},}\\\
&{\,\big{|}u(\zeta_{k}-t)\big{|}>\big{|}u(\zeta_{k}+t)\big{|},\quad
0<t<\zeta_{k+1}-\zeta_{k},}\end{aligned}\right.$

where $\big{(}\zeta_{k}\big{)}$ denotes all zeros of
$u(t),\,0\equiv\zeta_{0}<\zeta_{1}<\zeta_{2}<\cdots.$

###### Proof.

As $\phi(t)$ is strictly increasing, either $\phi(t)$ has a unique zero
$\,a>0\,$ or $\,\phi(t)\geq 0\,$ for all $\,t>0.$ In the former case, since it
is known that $u(t)$ has at most one zero for $\,0<t<a\,$ due to the sign
condition $\,\phi(t)<0\,$ (see, for instance, [11]), the stated result of part
(i) follows from Lemma 2.1.

In the latter case, we note that $\,\phi(0+)\,$ exists and hence

$\limsup_{t\,\to\,0+}\big{[}t^{2}\phi(t)\big{]}=0<\frac{1}{4},$

that is, condition (B1) is verified. As a consequence, for any fixed $\,b>0,$
$u(t)$ is non-oscillatory on the interval $\,(0,b]\,$ and oscillatory on the
interval $\,[b,\infty)\,$ so that part (i) follows once again by Lemma 2.1.

By considering $-u(t)$ otherwise, we may assume $\,u(t)>0$ for $0<t<\zeta_{1}$
and hence $u(t)$ is positive on $\,\big{(}\zeta_{2k},\,\zeta_{2k+1}\big{)}\,$
and negative on $\,\big{(}\zeta_{2k+1},\,\zeta_{2k+2}\big{)}\,$ for each
$\,k=0,1,\cdots.$ We fix $\,k\geq 1\,$ and consider auxiliary functions

$\displaystyle v(t)$ $\displaystyle=(-1)^{k+1}u\left(\zeta_{k}-t\right),\quad
0\leq t<\zeta_{k}-\zeta_{k-1},$ $\displaystyle w(t)$
$\displaystyle=(-1)^{k}u\left(\zeta_{k}+t\right),\quad 0\leq
t<\zeta_{k+1}-\zeta_{k},$

both of which are positive. Evidently, $v(t)$ and $w(t)$ satisfy

$\displaystyle\quad v^{\prime\prime}+\phi\left(\zeta_{k}-t\right)v=0,\quad
w^{\prime\prime}+\phi\left(\zeta_{k}+t\right)w=0,$ $\displaystyle
v(0)=w(0)=0,\quad
v^{\prime}(0)=w^{\prime}(0)=(-1)^{k}u^{\prime}(\zeta_{k})>0,$

whence part (ii) follows at once from Sturm’s theorem (II). ∎

######  Remark 2.2.

When $\phi(t)$ is increasing, not necessarily in the strict sense, an obvious
modification indicates that part (ii) must be replaced by

$\left\\{\begin{aligned}
&{\quad\zeta_{k}-\zeta_{k-1}\geq\zeta_{k+1}-\zeta_{k},}\\\
&{\,\big{|}u(\zeta_{k}-t)\big{|}\geq\big{|}u(\zeta_{k}+t)\big{|},\,\,0<t<\zeta_{k+1}-\zeta_{k}.}\end{aligned}\right.$
(2.1)

The inequalities of part (ii) as well as (2.1) are reversed in the case when
$\phi(t)$ is decreasing. Although these results, often referred to as Sturm’s
convexity theorem, are known and available in the aforementioned references,
we gave a proof for part (ii) for the sake of completeness.

##  3 Oscillatory integrals

Our main theorem reads as follows.

######  Theorem 3.1.

Suppose that $u(t)$ is an oscillatory solution of (1.1), where $\phi(t)$ is
continuous, strictly increasing and subject to one of assumptions (A1), (A2),
(A3). Let $\big{(}\zeta_{k}\big{)}$ denote all zeros of $u(t)$ arranged in
ascending order and assume that $\,u(0+)=0$ and $\,u(t)>0\,$ for
$\,0<t<\zeta_{1}.$ Put

$(Uf)(x)=\int_{0}^{\infty}f(t)u(xt)\,dt\qquad(x>0).$

If $f(t)$ is a nonnegative, sectionally continuous and decreasing function
such that it does not vanish identically and satisfies

$\displaystyle{\rm(P1)}\quad\int_{0}^{1}f(t)u(t)\,dt<\infty,$
$\displaystyle{\rm(P2)}\quad\lim_{k\to\infty}\int_{\zeta_{k}}^{\zeta_{k+1}}f(t/x)|u(t)|\,dt=0\quad\text{for
each $\,x>0,$}$

then $(Uf)(x)$ exists for each $\,x>0\,$ with $\,(Uf)(x)>0.$

###### Proof.

The sign assumption on $u(t)$ for $\,0<t<\zeta_{1}\,$ implies

$\left\\{\begin{aligned} &{u(t)>0,\quad\zeta_{2k}<t<\zeta_{2k+1},}\\\
&{u(t)<0,\quad\zeta_{2k+1}<t<\zeta_{2k+2},}\end{aligned}\right.\,\,\,\,k=0,1,2,\cdots,$
(3.1)

where $\,\zeta_{0}\equiv 0.$ By Lemma 2.2,
$\,\zeta_{2k+1}-\zeta_{2k}>\zeta_{2k+2}-\zeta_{2k+1}\,$ and

$\displaystyle u(\zeta_{2k+1}-t)>\big{|}u(\zeta_{2k+1}+t)\big{|},\quad
0<t<\zeta_{2k+2}-\zeta_{2k+1}.$ (3.2)

Let us fix $\,x>0\,$ and put $\,g(t)=f(t/x),\,t>0.$ We decompose

$\displaystyle x(Uf)(x)$
$\displaystyle=\sum_{k=0}^{\infty}\int_{\zeta_{k}}^{\zeta_{k+1}}g(t)u(t)dt=\sum_{k=0}^{\infty}B_{k},$
$\displaystyle B_{k}$
$\displaystyle:=\left(\int_{\zeta_{2k}}^{\zeta_{2k+1}}+\int_{\zeta_{2k+1}}^{\zeta_{2k+2}}\right)g(t)u(t)dt,$

which is legitimate if $(Uf)(x)$ exists. By changing variables, we write

$\displaystyle B_{k}$
$\displaystyle=\int_{0}^{\zeta_{2k+1}-\zeta_{2k}}g(\zeta_{2k+1}-t)u(\zeta_{2k+1}-t)\,dt$
$\displaystyle+\int_{0}^{\zeta_{2k+2}-\zeta_{2k+1}}g(\zeta_{2k+1}+t)u(\zeta_{2k+1}+t)\,dt$
$\displaystyle=\int_{0}^{\zeta_{2k+2}-\zeta_{2k+1}}\Big{[}g(\zeta_{2k+1}-t)u(\zeta_{2k+1}-t)+g(\zeta_{2k+1}+t)u(\zeta_{2k+1}+t)\Big{]}dt$
$\displaystyle+\int_{\zeta_{2k+2}-\zeta_{2k+1}}^{\zeta_{2k+1}-\zeta_{2k}}g(\zeta_{2k+1}-t)u(\zeta_{2k+1}-t)\,dt.$
(3.3)

Since $g(t)$ is nonnegative, decreasing and sectionally continuous, (3.1),
(3.2) clearly indicate that both integrands of (3) are positive when the
interval $\big{[}\zeta_{2k},\,\zeta_{2k+2}\big{]}$ contains the support of
$g(t)$ and nonnegative in the general case. Thus it is evident that
$\,B_{k}>0\,$ in the former case and $\,B_{k}\geq 0\,$ for all $k$. By
modifying the support condition in an obvious way, we may conclude that
$\,(Uf)(x)>0\,$ unless $g(t)$ vanishes identically.

It remains to verify that $(Uf)(x)$ exists, that is, the integral defining
$(Uf)(x)$ converges. For this purpose, we decompose

$x(Uf)(x)=\sum_{k=0}^{\infty}(-1)^{k}A_{k},\quad
A_{k}:=\int_{\zeta_{k}}^{\zeta_{k+1}}g(t)|u(t)|dt.$ (3.4)

Since $f(t)$ is monotone and $u(t)$ is continuous for $\,t>0,$ the condition
(P1) implies that the integral of $f(t)|u(xt)|$, possibly improper, over any
right-neighborhood of $0$ exists. It is thus evident that

$A_{0}=x\int_{0}^{\zeta_{1}/x}f(t)|u(xt)|dt<\infty.$

We may assume that $\,A_{k}>0\,$ for all $k$ for the convergence is obvious
otherwise. Proceeding as above, we observe that

$\displaystyle A_{k}-A_{k+1}=$
$\displaystyle(-1)^{k}\int_{0}^{\zeta_{k+2}-\zeta_{k+1}}\Big{[}g(\zeta_{k+1}-t)u(\zeta_{k+1}-t)$
$\displaystyle\qquad\qquad\qquad\quad+\,g(\zeta_{k+1}+t)u(\zeta_{k+1}+t)\Big{]}dt$
$\displaystyle+(-1)^{k}\int^{\zeta_{k+1}-\zeta_{k}}_{\zeta_{k+2}-\zeta_{k+1}}g(\zeta_{k+1}-t)u(\zeta_{k+1}-t)\,dt>0,$

where $\,k=0,1,2,\cdots,$ and so $(A_{k})$ is strictly decreasing. Since
$\,A_{k}\to 0\,$ as $\,k\to\infty\,$ by condition (P2), it follows by the
alternating series test that the series on the right of (3.4) converges and
thus $(Uf)(x)$ exists. ∎

######  Remark 3.1.

If we assume instead that $\,u(t)<0\,$ for all $\,0<t<\zeta_{1},$ without
altering other assumptions, then the same proof verifies

$\,(Uf)(x)<0\quad\text{for all}\quad x>0.$

Suppose that $\phi(t)$ is merely increasing so that (2.1) holds true. In this
case, a close inspection on (3) shows that both integrands are nonnegative so
that $\,B_{k}\geq 0\,$ for all $k$ under the same assumption on $g(t)$. If the
interval $\big{[}\zeta_{2k},\,\zeta_{2k+2}\big{]}$ contains the support of
$g(t)$, where $g(t)$ is assumed to be strictly decreasing, then the first
integrand is positive and hence $\,B_{k}>0.$

######  Theorem 3.2.

Let $u(t)$ be an oscillatory solution of (1.1), where $\phi(t)$ is continuous,
increasing, not necessarily in the strict sense, and subject to one of
assumptions (A1), (A2), (A3), such that $\,u(0+)=0$ and $\,u(t)>0\,$ for all
$\,0<t<\zeta_{1}.$ If $f(t)$ is a nonnegative and sectionally continuous
function such that it is strictly decreasing on its support and subject to
conditions (P1), (P2), then $(Uf)(x)$ exists for each $\,x>0\,$ with
$\,(Uf)(x)>0.$

##  4 Fourier transforms

This section aims to refine Pólya’s results concerning positivity of Fourier
transforms by revisiting in the framework of oscillatory integrals.

As for the Fourier sine transform on the positive real line

$(Uf)(x)=\int_{0}^{\infty}f(t)\sin\,(xt)\,dt\qquad(x>0),$

the kernel $\,u(t)=\sin t\,$ is an oscillatory solution of the differential
equation (1.1) with $\,\phi(t)\equiv 1\,$ such that $\,u(0+)=0\,$ and
$\,u(t)>0\,$ for $\,0<t<\pi.$ Thus the Fourier sine transform falls under the
scope of Theorem 3.2.

  * •

Due to the asymtotic behavior $\,\sin t\,\sim\,t\,$ as $\,t\to 0,$ it is
evident that condition (P1) holds when $\,f(t)t\,$ is integrable over the
interval $(0,1).$

  * •

We claim that condition (P2) holds if $f(t)$ is a nonnegative decreasing
function such that $\,f(t)\to 0\,$ as $\,t\to\infty.$ Indeed, given $\,x>0\,$
and $\,\epsilon>0,\,$ if we choose $N$ so that $\,f(t/x)<\epsilon/2\,$ for all
$\,t\geq N\pi,$ then

$\int_{k\pi}^{(k+1)\pi}f(t/x)|\sin
t|\,dt\leq\frac{\epsilon}{2}\int_{k\pi}^{(k+1)\pi}|\sin t|\,dt=\epsilon,\quad
k\geq N.$

To sum up, Theorem 3.2 yields the following result.

######  Theorem 4.1.

The inequality

$\int_{0}^{\infty}f(t)\sin(xt)\,dt>0$

holds true for all $\,x>0\,$ when $f(t)$ is nonnegative, sectionally
continuous, strictly decreasing on its support,

$\int_{0}^{1}f(t)t\,dt<\infty\quad\text{and}\quad\lim_{t\to\infty}f(t)=0.$

We illustrate Theorem 4.1 with

$\int_{0}^{\infty}\frac{\sin(xt)}{\,t^{\gamma}(t^{2}+a^{2})^{\delta}}\,dt>0,\quad
x>0,\,a>0,$

valid for $\,0\leq\gamma<2,\,\delta>0\,$ or $\,0<\gamma<2,\,\delta=0,$ which
includes the following special cases of interest ([15, (20), (28), 2.2])

$\displaystyle\int_{0}^{\infty}\frac{\sin(xt)}{t(t^{2}+a^{2})}\,dt$
$\displaystyle=\frac{\pi}{2a^{2}}\left(1-e^{-ax}\right),$
$\displaystyle\int_{0}^{\infty}\frac{\sin(2xt)}{\sqrt{t(t^{2}+a^{2})}}\,dt$
$\displaystyle=\sqrt{\pi x}\,I_{1/4}(ax)K_{1/4}(ax),$

where $\,I_{\nu}(x),K_{\nu}(x)\,$ stand for the modified Bessel functions of
the first and second kind, respectively (see, for instance, [4, 3.7]).

An integration by parts gives

$\int_{0}^{\infty}f(t)\cos(xt)\,dt=\frac{1}{x}\int_{0}^{\infty}\left[-f^{\prime}(t)\right]\sin(xt)\,dt,$

provided that $\,tf(t)\to 0\,$ as $\,t\to 0+\,$ and $\,f(t)\to 0\,$ as
$\,t\to\infty.$ By applying Theorem 3.2 to $\,-f^{\prime}(t),$ we obtain the
following result. As customary, a function $h$ will be said to be sectionally
smooth on an interval $I$ if both $\,h,h^{\prime}\,$ are sectionally or
piecewise continuous on $I$.

######  Corollary 4.1.

The inequality

$\int_{0}^{\infty}f(t)\cos(xt)\,dt>0$

holds true for all $\,x>0\,$ when $f(t)$ is nonnegative, sectionally smooth,
strictly decreasing and convex on its support,

$\int_{0}^{1}f(t)\,dt<\infty\quad\text{and}\quad\lim_{t\to\infty}f(t)=0.$

We illustrate Corollary 4.1 with

$\int_{0}^{\infty}\frac{\cos(xt)}{(t+a)^{\lambda}}\,dt>0,\quad x>0,\,a>0,$

valid for any $\,\lambda>0.$ As this example may indicate, it should be noted
that Theorem 4.1 and Corollary 4.1 extend Pólya’s results to the case of non-
integrable functions. Additionally relevant examples are

$\displaystyle\int_{0}^{\infty}t^{-\nu}\sin(xt)\,dt$
$\displaystyle=\frac{\pi}{2\Gamma(\nu)}\csc\left(\frac{\nu\pi}{2}\right)x^{\nu-1}>0,$
$\displaystyle\int_{0}^{\infty}t^{-\nu}\cos(xt)\,dt$
$\displaystyle=\frac{\pi}{2\Gamma(\nu)}\sec\left(\frac{\nu\pi}{2}\right)x^{\nu-1}>0,$

valid for all $\,x>0,$ where the first inequality holds true for $\,0<\nu<2\,$
and the second inequality for $\,0<\nu<1\,$ (see [15, (1), 1.3, (1), 2.3] ).
As obvious, the function $\,f(t)=t^{-\nu}\,$ is not integrable for any
$\,\nu\in\mathbb{R}.$

##  5 Hankel transforms

In the framework of Sturm’s oscillation theory, the kernel
$\,u(t)=\sqrt{t}J_{\nu}(t)\,$ of Hankel transform (1.4) arises as a solution
of (1.1) with

$\phi(t)=1+\frac{1-4\nu^{2}}{4t^{2}},\quad t>0.$ (5.1)

As a smooth function, $\phi(t)$ satisfies all oscillation assumptions (A1),
(A2), (A3). Moreover, it decreases strictly to the value $1$ when
$\,|\nu|<1/2\,$ and increases strictly to the value $1$ when $\,|\nu|>1/2.$ As
customary, $\big{(}j_{\nu,k}\big{)}$ will denote all positive zeros of
$J_{\nu}(z)$ arranged in ascending order.

Our principal theorem states as follows.

######  Theorem 5.1.

For $\,\nu>1/2,$ we have $\,(H_{\nu}f)(x)>0\,$ for all $\,x>0\,$ if $f(t)$ is
a nonnegative, sectionally continuous and decreasing function for $\,t>0\,$
such that it does not vanish identically,

$\int_{0}^{1}t^{\nu+1/2}f(t)dt<\infty\quad\textit{and}\quad\lim_{t\to\infty}f(t)=0.$

###### Proof.

On considering the asymptotic behavior

$\displaystyle u(t)=\sqrt{t}J_{\nu}(t)$
$\displaystyle=\sqrt{t}\,\sum_{m=0}^{\infty}\frac{(-1)^{m}}{\,m!\,\Gamma(\nu+1+m)\,}\left(\frac{t}{2}\right)^{2m+\nu}$
$\displaystyle\,\sim\,\frac{\sqrt{2}}{\Gamma(\nu+1)}\left(\frac{t}{2}\right)^{\nu+1/2}\quad\text{as
$\,t\to 0,$}$

it is clear that $\,u(0+)=0,\,u(t)>0\,$ for $\,0<t<j_{\nu,1}\,$ and condition
(P1) holds true. By Theorem 3.1, hence, it suffices to show that condition
(P2) holds true under the stated assumptions on $f(t)$.

Since $\phi(t)$ increases strictly to the value $1$, if we apply Sturm’s
theorem (I) with any $m$ satisfying $\,0<m<\phi(t)<1\,$ and then letting
$\,m\to 1,$ it is straightforward to deduce that

$j_{\nu,2}-j_{\nu,1}>j_{\nu,3}-j_{\nu,2}>\cdots\,\to\,\pi.$ (5.2)

Hankel’s asymptotic expansion ([4, 7.21]) indicates

$u(t)=\sqrt{\frac{2}{\pi}}\,\cos\left(t-\frac{\nu\pi}{2}-\frac{\pi}{4}\right)+O\big{(}t^{-1}\big{)}$
(5.3)

as $\,t\to\infty.$ By combining (5.2) with McMahon’s asymptotic expansion for
large zeros of Bessel functions ([16], [4, 15.53])

$j_{\nu,k}=\left(k+\frac{\nu}{2}-\frac{1}{4}\right)\pi+O\big{(}k^{-1}\big{)}$
(5.4)

as $\,k\to\infty,$ we may choose an integer $N$ so that

$j_{\nu,k}\geq\frac{k\pi}{2},\quad\pi<j_{\nu,k+1}-j_{\nu,k}<2\pi\quad\text{for
all $\,k\geq N.$}$

For $\,k\geq N,$ these inequalities clearly imply that

$\displaystyle\int_{j_{\nu,k}}^{j_{\nu,k+1}}\left|\cos\left(t-\frac{\nu\pi}{2}-\frac{\pi}{4}\right)\right|\,dt\leq
2\pi,\quad\int_{j_{\nu,k}}^{j_{\nu,k+1}}\frac{\,dt\,}{t}\leq\frac{4}{k},$

whence it follows from Hankel’s formula (5.3) that

$\lim_{k\to\infty}\int_{j_{\nu,k}}^{j_{\nu,k+1}}f(t/x)|u(t)|\,dt=0,\quad x>0,$

under the assumption $\,f(t)\to 0\,$ as $\,t\to\infty.$ It is thus verified
that (P2) holds true and Theorem 3.1 gives the desired positivity. ∎

As readily observed, Theorem 5.1 gives

$\int_{0}^{\infty}\frac{J_{\nu}(xt)\sqrt{xt}\,}{\,t^{\gamma}(t^{2}+a^{2})^{\delta}}\,dt\,>0,\quad
x>0,\,a>0,$

valid for $\,\nu\geq 1/2\,$ and $\,\gamma=0,\,\delta>0\,$ or
$\,0<\gamma<\nu+3/2,\,\delta\geq 0,$ which includes the following special
cases [15, (3), (11), (17), 8.5] :

$\displaystyle\int_{0}^{\infty}J_{\nu}(xt)\,dt$ $\displaystyle=\frac{1}{x},$
$\displaystyle\int_{0}^{\infty}\frac{J_{\nu}(2xt)\,dt\,}{\sqrt{t^{2}+a^{2}}}\,$
$\displaystyle=I_{\nu/2}(ax)K_{\nu/2}(ax),$
$\displaystyle\int_{0}^{\infty}\frac{J_{\nu}(2xt)\,dt\,}{\,t^{\nu}(t^{2}+a^{2})^{\nu+1/2}}\,$
$\displaystyle=\frac{(2/a)^{2\nu}\Gamma(\nu+1)}{\Gamma(2\nu+1)}\,x^{\nu}I_{\nu}(ax)K_{\nu}(ax).$

By applying instead Theorem 5.3, to be established in the sequel, it is
possible to extend the stated range of positivity to the case $\,\nu>-1.$

###  5.1 Alternative: The case $\,0<\nu<1/2$

While we have investigated positivity of Hankel transforms in the cases
$\,\nu>1/2\,$ or $\,\nu=\pm 1/2\,$ on the basis of comparing the wave of
arches, such a method is no longer suitable in dealing with other cases due to
our restrictive assumptions on the kernel and $\phi(t)$. As an alternative, we
consider the kernel of type $\,u(t)=\sqrt{t}J_{\nu}\left(t^{\alpha}\right)\,$
which arises as a solution of (1.1) with

$\displaystyle\phi(t)=\frac{1}{t^{2}}\left(\alpha^{2}t^{2\alpha}+\frac{1-4\alpha^{2}\nu^{2}}{\,4\,}\right),\quad
t>0.$ (5.5)

When $\,\alpha>1,\,\alpha|\nu|\geq 1/2,$ $\phi(t)$ is strictly increasing.
Since

$u(t)=\sqrt{t}J_{\nu}\left(t^{\alpha}\right)\,\sim\,\frac{t^{\alpha\nu+1/2}\,}{2^{\nu}\Gamma(\nu+1)}\quad\text{as
$\,t\to 0+,\,$}$

if $\,\alpha\nu+1/2>0,$ then $\,u(0+)=0,\,u(t)>0\,$ for
$\,0<t<j_{\nu,1}^{1/\alpha}\,$ and (P1) is equivalent to the integrability of
$\,f(t)t^{\alpha\nu+1/2}\,$ for $\,0<t<1.$

In view of Hankel’s asymptotic formula

$u(t)=\sqrt{\frac{2}{\pi}}\,t^{\frac{1-\alpha}{2}}\left[\cos\left(t^{\alpha}-\frac{\nu\pi}{2}-\frac{\pi}{4}\right)+O\big{(}t^{-\alpha}\big{)}\right]$

as $\,t\to\infty\,$ and a modification of (5.2) in the case $\,|\nu|\leq 1/2,$

$j_{\nu,2}-j_{\nu,1}\leq j_{\nu,3}-j_{\nu,2}\leq\cdots\,\to\,\pi,$

it is easy to confirm that (P2) holds with no further assumptions on $f(t)$.

It is thus shown by Theorem 3.1 that

$(Uf)(x)=\int_{0}^{\infty}f(t)J_{\nu}\big{[}(xt)^{\alpha}\big{]}\sqrt{xt}\,dt>0\qquad(x>0),$

valid for $\,\nu>0,\,\alpha>1\,$ with $\,\alpha\nu\geq 1/2,$ provided that
$f(t)$ is a nonnegative, sectionally continuous and decreasing function
satisfying

$\int_{0}^{1}f(t)\,t^{\alpha\nu+1/2}\,dt<\infty.$

By changing variables appropriately, we find that

$\alpha
x^{\frac{\alpha-1}{2\alpha}}(Uf)(x^{\frac{1}{\alpha}})=\int_{0}^{\infty}f(t^{\frac{1}{\alpha}})t^{\frac{3(1-\alpha)}{2\alpha}}J_{\nu}(xt)\sqrt{xt}\,dt>0.$

In the case $\,0<\nu<1/2,$ if we choose $\,\alpha=1/2\nu,$ then this
alternative method yields an extension of Theorem 5.1 which reads as follows.

######  Theorem 5.2.

For $\,0<\nu<1/2,$ we have $\,(H_{\nu}f)(x)>0\,$ for all $\,x>0\,$ if $f(t)$
is a nonnegative sectionally continuous function such that it does not vanish
identically, $\,t^{3/2-3\nu}f(t)\,$ is decreasing for $\,t>0\,$ and

$\int_{0}^{1}t^{\nu+1/2}f(t)\,dt<\infty.$

######  Remark 5.1.

In [17], Makai exploited the scaled function
$\,u(t)=\sqrt{t}J_{\nu}\left(\beta t^{1/2\nu}\right)\,$ in proving that
$\,j_{\nu,k}/\nu\,$ is a strictly decreasing function of $\,\nu>0.$

We illustrate with Gegenbauer’s formula

$\int_{0}^{\infty}e^{-bt}J_{\nu}(xt)\frac{dt}{t}=\frac{\big{[}\sqrt{x^{2}+b^{2}}-b\big{]}^{\nu}\,}{\nu
x^{\nu}}>0,$ (5.6)

valid for all $\,x>0\,$ and $\,b>0,\,\nu>0\,$ (see [15, (2), 8.6], [4, 13.2]).

  * •

In the case $\,\nu\geq 1/2,$ we identify

$\displaystyle\int_{0}^{\infty}e^{-bt}J_{\nu}(xt)\frac{dt}{t}=\frac{1}{\sqrt{x}}\big{(}H_{\nu}f\big{)}(x),\,\,f(t)=t^{-\frac{3}{2}}e^{-bt}.$

As elementary, $f(t)$ is positive, continuous, strictly decreasing and
$\,f(t)\to 0\,$ as $\,t\to\infty.$ In addition, it is trivial to see that

$\int_{0}^{1}t^{\nu+1/2}f(t)dt=\int_{0}^{1}t^{\nu-1}e^{-bt}dt<\infty.$

Thus it follows by theorems 4.1, 5.1 that $\,\big{(}H_{\nu}f\big{)}(x)>0.$

  * •

In the case $\,0<\nu<1/2,$ we note that the associated function

$t^{3/2-3\nu}f(t)=t^{-3\nu}e^{-bt},\quad t>0,$

is decreasing and hence $\,\big{(}H_{\nu}f\big{)}(x)>0\,$ by Theorem 5.2.

###  5.2 The general case $\,\nu>-1$

By exploiting the recurrence relation of Bessel functions, it is possible to
deduce positivity of Hankel transforms in the general case $\,\nu>-1\,$ on the
basis of the preceding two theorems.

######  Theorem 5.3.

For $\,\nu>-1,$ suppose that $f(t)$ is a nonnegative, decreasing and
sectionally smooth function for $\,t>0\,$ satisfying

$\int_{0}^{1}t^{\nu+1/2}f(t)dt<\infty\quad\textit{and}\quad\lim_{t\to\infty}f(t)=0.$
(5.7)

Put

$\displaystyle g(t)$
$\displaystyle=-t^{\nu+1/2}\frac{d}{dt}\left[t^{-\nu-1/2}f(t)\right]$
$\displaystyle=\left(\nu+\frac{1}{2}\right)\frac{f(t)}{t}-f^{\prime}(t),\quad
t>0.$

If both $f(t),g(t)$ do not vanish identically for $\,t>0,$ then
$\,(H_{\nu}f)(x)>0\,$ for all $\,x>0\,$ under each of the following additional
case-assumptions.

  * (i)

For $\,\nu>-1/2,$ $g(t)$ is nonnegative, decreasing and
$\,\lim_{t\to\infty}g(t)=0.$

  * (ii)

For $\,\nu=-1/2,$ $f(t)$ is strictly decreasing and convex on its support.

  * (iii)

For $\,\nu<-1/2,$ $g(t)$ is nonnegative and $t^{-3/2-3\nu}g(t)$ is decreasing.

###### Proof.

By using the differentiation formula [4, 3.2]

$\frac{d}{dt}\left[t^{\nu+1}J_{\nu+1}(t)\right]=t^{\nu+1}J_{\nu}(t),\quad
t>0,$

we fix $\,x>0\,$ and integrate by parts to obtain

$\displaystyle\left(H_{\nu}f\right)(x)$
$\displaystyle=\frac{1}{\sqrt{x}}\int_{0}^{\infty}t^{-\nu-1/2}f(t)\frac{d}{dt}\left[t^{\nu+1}J_{\nu+1}(xt)\right]dt$
$\displaystyle=\frac{1}{\sqrt{x}}\,\Big{[}t^{1/2}f(t)J_{\nu+1}(xt)\Big{]}_{t=0}^{\infty}$
$\displaystyle\quad\qquad-\frac{1}{x}\int_{0}^{\infty}t^{\nu+1/2}\frac{d}{dt}\left[t^{-\nu-1/2}f(t)\right]J_{\nu+1}(xt)\sqrt{xt}\,dt$
$\displaystyle=\frac{1}{x}\int_{0}^{\infty}g(t)J_{\nu+1}(xt)\sqrt{xt}\,dt=\frac{1}{x}\left(H_{\nu+1}g\right)(x),$
(5.8)

provided that $\,G(t):=t^{1/2}f(t)J_{\nu+1}(xt)\,$ tends to zero as $t$
approaches zero or infinity. To justify, we first note that
$\,t^{1/2}J_{\nu+1}(xt)\,$ is bounded for large $t$, due to Hankel’s formula
(5.3), and hence the second assumption of (5.7) indicates $\,G(t)\to 0\,$ as
$\,t\to\infty.$ Next, by inspecting the first assumption of (5.7) and
monotonicity of $f(t)$, it is not difficult to see that $\,t^{\nu+3/2}f(t)\to
0\,$ as $\,t\to 0+.$ By combining with the initial behavior

$t^{1/2}J_{\nu+1}(xt)\,\sim\,\frac{x^{\nu+1}t^{\nu+3/2}}{\Gamma(\nu+1)2^{\nu+1}\,}\quad\text{as
$\,t\to 0+,$}$

we find that $\,G(t)\to 0\,$ as $\,t\to 0+\,$ and thus (5.2) is legitimate.

Concerning part (i), Theorem 5.1 asserts that $\,(H_{\nu+1}g)(x)>0\,$ when
$g(t)$ is nonnegative, sectionally continuous, decreasing with

$\int_{0}^{1}t^{\nu+3/2}g(t)dt<\infty\quad\text{and}\quad\lim_{t\to\infty}g(t)=0.$

On integrating by parts, it is easy to observe that

$\int_{0}^{1}t^{\nu+3/2}g(t)dt=-f(1)+(2\nu+2)\int_{0}^{1}t^{\nu+1/2}f(t)dt<\infty$
(5.9)

and the result follows. Part (ii) is equivalent to Corollary 4.1. In the case
$\,\nu<-1/2,$ that is, $\,0<\nu+1<1/2,$ since Theorem 5.2 indicates that
$\,(H_{\nu+1}g)(x)>0\,$ when $g(t)$ is nonnegative, sectionally continuous,
subject to condition (5.9) and $\,t^{-3/2-3\nu}g(t)\,$ is decreasing, the
result follows. ∎

As an illustration, we take $\,f(t)=t^{-\beta},\,t>0,$ and consider

$\displaystyle(H_{\nu}f)(x)$
$\displaystyle=\int_{0}^{\infty}t^{-\beta}J_{\nu}(xt)\sqrt{xt}\,dt$
$\displaystyle=2^{-\beta+1/2}x^{\beta-1}\frac{\,\Gamma\big{(}(2\nu-2\beta+3)/4\big{)}\,}{\,\Gamma\big{(}(2\nu+2\beta+1)/4\big{)}\,},$
(5.10)

valid for $\,\nu>-1,\,0<\beta<\nu+3/2\,$ ([15, 8.5, (7)], [18, 1.4, (4.6)]).
Since

$g(t)=(\nu+\beta+1/2)t^{-\beta-1},\quad t>0,$

it is elementary to find by applying Theorem 5.3 that the Hankel transform
$(H_{\nu}f)(x)$ given by (5.2) is positive for all $\,x>0\,$ if

$\left\\{\begin{aligned} &{\,\,\,0<\beta<\nu+3/2\quad\text{for
$\,\,\,\nu\geq-1/2,$}}\\\ &{-\nu-1/2<\beta<\nu+3/2\quad\text{for
$\,\,-1<\nu<-1/2,$}}\end{aligned}\right.$

the optimal range of $\,\nu,\beta\,$ for positivity in this example.

##  6 Y-transforms

As another oscillatory solution of the differential equation (1.1) with (5.1),
if we take $\,u(t)=\sqrt{t}\,Y_{\nu}(t),$ where $Y_{\nu}(t)$ denotes a Neumann
function or the Bessel function of the second kind defined by

$Y_{\nu}(t)=\frac{\,\cos{\nu\pi}J_{\nu}(t)-J_{-\nu}(t)\,}{\sin{\nu\pi}},$

then the corresponding oscillatory integral takes the form

$\displaystyle(N_{\nu}f)(x)$
$\displaystyle=\int_{0}^{\infty}f(t)Y_{\nu}(xt)\sqrt{xt}\,dt$
$\displaystyle=\cot\nu\pi\,(H_{\nu}f)(x)-\csc\nu\pi\,\left(H_{-\nu}f\right)(x),$
(6.1)

often referred to as the $Y$-transform of order $\nu$ of $f(t)$ (see [15],
[18]).

Concerning positivity, we remark that Theorem 3.1 is not applicable to
$Y$-transforms. Indeed, while $\phi(t)$ is increasing only when $\,|\nu|\geq
1/2,$ it is easy to see that $\,|u(t)|=\sqrt{t}\,|Y_{\nu}(t)|\,$ tends to
infinity as $t$ approaches zero when $\,|\nu|>1/2.$ Nevertheless, it is
possible to deduce positivity or negativity of $(N_{\nu}f)(x)$ by making use
of the relation (6) in the case $\,1/2<|\nu|<1,$ provided that both
$\,(H_{\nu}f)(x),\,(H_{-\nu}f)(x)\,$ are positive at the same time: It is
positive for $\,-1<\nu<-1/2\,$ and negative for $\,1/2<\nu<1.$

In practice, it is straightforward to obtain the following positivity result
for $Y$-transforms on applying Theorem 3.1 and Theorem 5.3.

######  Theorem 6.1.

For $\,1/2<|\nu|<1,$ let $f(t)$ be a nonnegative, decreasing and sectionally
continuous function such that it does not vanish identically,

$\int_{0}^{1}t^{-|\nu|+1/2}f(t)dt<\infty\quad\textit{and}\quad\lim_{t\to\infty}f(t)=0.$

If $g(t)$ is nonnegative and $\,t^{-3/2+3|\nu|}g(t)\,$ is decreasing for
$\,t>0,$ where

$g(t)=\left(-|\nu|+\frac{1}{2}\right)\frac{f(t)}{t}-f^{\prime}(t),\quad t>0,$

then the following inequalities hold true for each $\,x>0:$

$\left\\{\begin{aligned} &{\,\,\,(N_{\nu}f)(x)>0\quad\text{for
$\,\,-1<\nu<-1/2,$}}\\\ &{\,\,\,(N_{\nu}f)(x)<0\quad\text{for
$\,\,1/2<\nu<1.$}}\end{aligned}\right.$

Acknowledgements. We are grateful to anonymous referees for their suggestions
which helped greatly improve the original version of the present paper.

Funding. The research by Yong-Kum Cho was supported by the National Research
Foundation of Korea funded by the Ministry of Science and ICT
(2021R1A2C1007437). The research by Young Woong Park was supported by the
Chung-Ang University Graduate Research Scholarship in 2021.

Disclosure Statement. No potential conflict of interest was reported by the
authors.

## References

  * [1] Sturm C. Mémoire sur les équations différentielles du second order. J Math Pures Appl. 1834;1:106-186.
  * [2] Pólya G. Über die nullstellen gewisser ganzer funktionen. Math Z. 1918;2:352-383.
  * [3] Tuck EO. On positivity of Fourier transforms. Bull Austral Math Soc. 2006;74:133–138.
  * [4] Watson GN. A Treatise on the Theory of Bessel Functions. London: Cambridge University Press; 1922.
  * [5] Zhang R. Positivity and Hankel transforms. Proc Amer Math Soc. 2022;150:3025-3035.
  * [6] Schaback R, Wu Z. Operators on radial functions. J Comput Appl Math. 1996;73:257-270.
  * [7] Wendland H. Scattered Data Approximation. Cambridge (UK): Cambridge University Press; 2005.
  * [8] Kneser A. Untersuchungen über die reellen nullstellen der integrale linearer differentialgleichungen. Math Ann. 1893;42:409-435.
  * [9] Fite WB. Concerning the zeros of the solutions of certain differential equations. Trans Amer Math Soc. 1918;19:755-759.
  * [10] Hille E. Non-oscillation theorems. Trans Amer Math Soc. 1948;64:234-252.
  * [11] Simmons GF. Differential Equations with Applications and Historical Notes. CRC Press: 3rd edition; 2017.
  * [12] Walker BK. Oscillation criteria of the Hille type for ordinary differential equations. Master thesis: Texas Tech University; 1970.
  * [13] Laforgia A, Muldoon M. Some consequences of the Sturm comparison theorem. Amer Math Monthly. 1986;93:89–94
  * [14] Makai E. On a monotonic property of certain Sturm-Liouville functions. Acta Math Acad Sci Hungar. 1952;3:165-172.
  * [15] Erdélyi A, Magnus W, Oberhettinger F, Tricomi FG. Tables of Integral Transforms, Vol I, II. McGraw-Hill; 1954.
  * [16] McMahon J. On the roots of the Bessel and certain related functions. Annals of Math. 1895;9:23-30.
  * [17] Makai E. On zeros of Bessel functions. Univ Beograd Publ Elektrotehn Fak Ser Mat Fiz No 602-633. 1978; 109-110.
  * [18] Oberhettinger F. Tables of Bessel Transforms. Springer-Verlag; 1972.