# Darboux transformation of diffusion processes

Alexey Kuznetsov, Minjian Yuan 111Dept. of Mathematics and Statistics, York
University, 4700 Keele Street, Toronto, ON, M3J 1P3, Canada.
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Darboux transformation of a second order linear differential operator is a
well-known technique with many applications in mathematics and physics. We
study Darboux transformation from the point of view of Markov semigroups of
diffusion processes. We construct a Darboux transform of a diffusion process
through a combination of Doob’s $h$-transform and a version of Siegmund
duality. Our main result is a simple formula that connects transition
probability densities of the two processes. We provide several examples of
Darboux transformed diffusion processes related to Brownian motion and
Ornstein-Uhlenbeck process. For these examples we compute explicitly the
transition probability density and derive its spectral representation.

Keywords: diffusion process, killed Brownian motion, Darboux transform, Sturm-
Liouville theory, Markov semigroup, Doob’s transform, Siegmund duality 2020
Mathematics Subject Classification : Primary 60J60, Secondary 60J35

##  1 Introduction

Darboux transformation is a well-known technique for studying second order
linear differential operators that goes back to the late 19th century [7]. To
introduce this transformation, we start with a second order linear
differential operator

${\mathcal{L}}=\frac{1}{2}\partial_{y}^{2}-c(y),$ (1)

where $y\in(l,r)\subseteq{\mathbb{R}}$. We assume that $c(y)$ is continuous on
$(l,r)$ and that ${\mathcal{L}}$ acts on $C^{2}$ functions on $(l,r)$. Let $h$
be a positive function on $(l,r)$ that satisfies ${\mathcal{L}}h=\lambda h$
for some $\lambda\in{\mathbb{R}}$. It is easy to check that the second order
differential operator ${\mathcal{L}}-\lambda$ can be factorized as a product
of two first order differential operators

${\mathcal{L}}-\lambda=\frac{1}{2}{\mathcal{D}}_{\frac{1}{h}}{\mathcal{D}}_{h},$
(2)

where we denoted

${\mathcal{D}}_{h}:=\partial_{y}-\frac{h^{\prime}(y)}{h(y)}=h(y)\partial_{y}\frac{1}{h(y)},\;\;\;{\mathcal{D}}_{\frac{1}{h}}:=\partial_{y}+\frac{h^{\prime}(y)}{h(y)}=\frac{1}{h(y)}\partial_{y}h(y).$
(3)

Now we introduce the Darboux transform of ${\mathcal{L}}$ as a second order
linear differential operator $\widetilde{\mathcal{L}}$ defined via

$\widetilde{\mathcal{L}}-\lambda=\frac{1}{2}{\mathcal{D}}_{h}{\mathcal{D}}_{\frac{1}{h}}.$
(4)

Thus, the Darboux transform of a second order linear operator is obtained by
factorizing it (or, more precisely, its shift by $\lambda$) as a product of
two first order differential operators and then interchanging the order of the
two factors. It is an easy exercise to verify that $\widetilde{\mathcal{L}}$
is given by

$\widetilde{\mathcal{L}}=\frac{1}{2}\partial_{y}^{2}-\tilde{c}(y),$ (5)

where

$\tilde{c}(y):=c(y)-\partial_{y}^{2}\ln
h(y)=\Big{(}\frac{h^{\prime}(y)}{h(y)}\Big{)}^{2}-c(y)-2\lambda.$ (6)

The key property of the Darboux transformation, which makes it so useful in
applications, is the following intertwining relation

$\widetilde{\mathcal{L}}{\mathcal{D}}_{h}={\mathcal{D}}_{h}{\mathcal{L}},$ (7)

which follows from (2) and (4). In particular, if $f$ is a solution to
${\mathcal{L}}f=\mu f$ then $g={\mathcal{D}}_{h}f$ solves
$\widetilde{\mathcal{L}}g=\mu g$. Note that
$\widetilde{\mathcal{L}}(1/h)=\lambda(1/h)$ (this follows immediately from
(4)), therefore the Darboux transformation is invertible: starting from the
operator $\widetilde{\mathcal{L}}$ and applying the Darboux transformation
with the function $\tilde{h}(y)=1/h(y)$ we will recover the original operator
${\mathcal{L}}$. The Darboux transformation can be iterated several times –
these are usually called Darboux-Crum transformations (see [6] and
[12][Section 2.4]).

A vast literature exists on the applications of the Darboux transformation in
applied mathematics and physics. For example, the Darboux transformation is
used to construct exactly solvable potentials for Schrödinger equation [1] and
it plays a crucial role in the study of exceptional orthogonal polynomials [9,
10]. Many other applications in integrable systems are described in [13] and
[18]. However, it seems that the Darboux transformation has not yet been
studied in the context of diffusion processes, and this is what we intend to
do in the present paper.

Assuming that the functions $c(y)$ and $\tilde{c}(y)$ in (1) and (5) are
positive, we can consider the operators ${\mathcal{L}}$ and
$\widetilde{\mathcal{L}}$ as the infinitesimal generators of killed Brownian
motion processes $Y$ and $\widetilde{Y}$ with state space $(l,r)$. When a
boundary point $l$ or $r$ is non-singular for one of these processes (see
Section 3 for a review of Feller’s boundary classification), we have many
different ways to construct a diffusion process with a given infinitesimal
generator by specifying different boundary conditions
(reflecting/killing/elastic/etc.). Instead of first specifying the boundary
conditions and then studying the resulting processes, we will proceed in the
opposite direction. We will first set a goal for what kind of process we want
to construct, and then we will see what boundary conditions and other
ingredients are needed for this construction to work. What follows next is an
informal description of what we plan to do. Let $Y$ be a killed Brownian
motion on $(l,r)$ with the Markov semigroup $\\{{\mathcal{P}}_{t}\\}_{t\geq
0}$ (where we denoted ${\mathcal{P}}_{t}f(y):={\mathbb{E}}_{y}[f(Y_{t})]$) and
with the infinitesimal generator ${\mathcal{L}}$ given by (1). Assume that we
found another killed Brownian motion process $\widetilde{Y}$ with the Markov
semigroup $\\{\widetilde{\mathcal{P}}_{t}\\}_{t\geq 0}$ and the infinitesimal
generator $\widetilde{\mathcal{L}}$ given by (5). Moreover, we assume that the
intertwining relationship for generators (7) carries over to a similar
relationship for semigroups, so that
$\widetilde{\mathcal{P}}_{t}{\mathcal{D}}_{h}={\mathcal{D}}_{h}{\mathcal{P}}_{t}$.
We want to use this intertwining relationship to express the transition
operators $\widetilde{\mathcal{P}}_{t}$ in terms of ${\mathcal{P}}_{t}$. One
way to define operators $\widetilde{\mathcal{P}}_{t}$ which satisfy such
intertwining relationship is via

$\widetilde{\mathcal{P}}_{t}={\mathcal{D}}_{h}{\mathcal{P}}_{t}{\mathcal{I}}^{(z)}_{h},$
(8)

where $z$ is some fixed point in $(l,r)$ and the integral operator
${\mathcal{I}}^{(z)}_{h}$ is defined as follows

${\mathcal{I}}^{(z)}_{h}f(y):=h(y)\int_{z}^{y}\frac{f(u)}{h(u)}{\textnormal{d}}u.$

Assume that $h$ is a positive $\lambda$-invariant function for
${\mathcal{P}}_{t}$, meaning that it satisfies ${\mathcal{P}}_{t}h=e^{\lambda
t}h$. This condition implies ${\mathcal{L}}h=\lambda h$. Then the operators
$\widetilde{\mathcal{P}}_{t}$ defined in (8) would indeed satisfy the
intertwining relationship
$\widetilde{\mathcal{P}}_{t}{\mathcal{D}}_{h}={\mathcal{D}}_{h}{\mathcal{P}}_{t}$,
since ${\mathcal{I}}^{(z)}_{h}{\mathcal{D}}_{h}f(y)=f(y)-f(z)h(y)/h(z)$ for
$C^{1}$ functions $f$ and ${\mathcal{D}}_{h}{\mathcal{P}}_{t}h=0$ (due to
$\lambda$-invariance of $h$). A similar argument shows that the operators
defined in (8) would not depend on the choice of $z\in(l,r)$. If it is indeed
true that the Markov semigroups of $Y$ and $\widetilde{Y}$ satisfy (8), then
it is easy to verify that the transition probability densities of
$\widetilde{Y}$ and $Y$ should satisfy:

$p_{t}^{\widetilde{Y}}(x,y)\stackrel{{\scriptstyle?}}{{=}}\frac{h(x)}{h(y)}\partial_{x}\Big{[}\frac{1}{h(x)}\int_{y}^{r}p_{t}^{Y}(x,u)h(u){\textnormal{d}}u\Big{]}.$
(9)

We added a question mark in the above equation to emphasize that this is only
a hypothesis at this stage. Indeed, the above discussion is non-rigorous, and
this construction of the process $\widetilde{Y}$ has many gaps. For example,
it is not clear why the operators (8) would indeed define a Markov semigroup,
in particular, why the right-hand side of (9) should be positive. It is
unclear whether we should impose the boundary conditions on $Y$ and
$\widetilde{Y}$ in some specific way for (9) to be true. However, this
informal reasoning is essentially correct, as we will demonstrate in this
paper. Our main result is the construction of a killed Brownian motion process
$\widetilde{Y}$ such the transition probabilities of $Y$ and $\widetilde{Y}$
satisfy an analogue of (9) (the difference is that we will allow adding a
constant to $\tilde{c}(y)$ that would make the resulting function positive,
and this will introduce an extra exponential factor in (9)). We will show that
the process $\widetilde{Y}$ can be constructed in three steps, starting from a
killed Brownian motion $Y$ and a $\lambda$-invariant function $h$: Doob’s
$h$-transform followed by a version of Siegmund duality [23] followed by
another Doob’s $h$-transform.

The primary motivation for constructing Darboux transformed diffusion
processes is that this would greatly expand the number of diffusions for which
the transition probability density can be computed in closed form. There are
only a few examples (Brownian motion, Ornstein-Uhlenbeck process, Bessel
process, square root diffusion process, and several others) for which the
transition probability is known in closed form, and all of these processes are
very useful in applications, particularly in Mathematical Finance. Expanding
the number of such examples is, therefore, a worthwhile task.

The paper is organized as follows. In Section 2 we recall some basic facts
about diffusion processes, including Feller’s classification of boundary
points for diffusion processes. In Section 3 we discuss Siegmund duality [3,
23] and present a new version of it (which can be applied to processes with
reflecting instead of absorbing boundary conditions). In Section 4 we
construct the Darboux transform of a killed Brownian motion and present an
example that illustrates the necessity of $\lambda$-invariance of $h$: we show
that if $h$ solves ${\mathcal{L}}h=\lambda h$ but is not $\lambda$-invariant
for the process $Y$, then the right-hand side of (9) may be negative. Section
5 provides five examples of the Darboux transformation applied to various
diffusion processes. In the first four examples we start with the Brownian
motion process on an interval with various boundary conditions, and in the
last example we apply the Darboux transformation to the Brownian motion killed
at rate $y^{2}/2$ (this latter process is closely related to the Ornstein-
Uhlenbeck process). In all examples we describe the boundary behavior and
derive explicit formulas for the transition probability density of the
transformed process. In four of these examples we also provide the spectral
representation of the transition probability density. At the end of Section 5
we discuss connections between our results and those obtained in the physics
literature on propagators for the one-dimensional Schrödinger equation.

##  2 Preliminaries

In this section we introduce notation and recall several basic facts about
diffusion processes, which will be useful later. We will mostly follow
[2][Chapter II]; other good references are [14, 17]. Let $X$ be a regular
diffusion process on an interval $(l,r)\subseteq{\mathbb{R}}$. We assume that
the infinitesimal generator of $X$ has the form

${\mathcal{L}}_{X}=\frac{1}{2}\sigma^{2}(x)\partial_{x}^{2}+b(x)\partial_{x}-c(x),$
(10)

where the functions $b$, $c$ and $\sigma$ are continuous on $(l,r)$ and
$\sigma(x)>0$ and $c(x)\geq 0$ for $x\in(l,r)$. The operator
${\mathcal{L}}_{X}$ is acting on the class of $C^{2}$ functions on $(l,r)$
with appropriate boundary conditions, which will be described below. We denote
by $s(x)$ the scale function of $X$ and by ${\textnormal{d}}m(x)$ and
${\textnormal{d}}k(x)$ the speed measure and the killing measure of $X$. The
derivative of $s$ and the densities of measures ${\textnormal{d}}m(x)$ and
${\textnormal{d}}k(x)$ can be found via

$s^{\prime}(x)=e^{-B(x)},\;\;m^{\prime}(x)=2\sigma^{-2}(x)e^{B(x)},\;\;\;\;k^{\prime}(x)=c(x)m^{\prime}(x),$
(11)

where $B^{\prime}(x):=2\sigma^{-2}(x)b(x)$, see [2][Section II.9]. We denote

$f^{\pm}(x)=\lim\limits_{h\to 0^{\pm}}\frac{f(x+h)-f(x)}{s(x+h)-s(x)}.$

Let us now recall Feller’s boundary classification (see [2][Section II.6]). We
fix a point $z\in(l,r)$ and denote

$R(x):=(m(z)-m(x)+k(z)-k(x))s^{\prime}(x),\;\;\;Q(x):=(s(z)-s(x))(m^{\prime}(x)+k^{\prime}(x)).$
(12)

The left boundary $l$ is called exit if $R\in L_{1}((l,z))$ and entrance if
$Q\in L_{1}((l,z))$. The conditions for the right boundary $r$ are the same,
except for an obvious change of the interval $(l,z)\mapsto(z,r)$. A boundary
point that is both not-entrance and not-exit is called natural while a
boundary point that is both entrance and exit is called non-singular. If $l$
(or $r$) is a non-singular boundary point for $X$, we impose one of the
following boundary conditions

  * (i)

elastic boundary condition: $f(l+)=\alpha f^{+}(l+)$ for $\alpha>0$ (or
$f(r-)=-\beta f^{-}(r-)$ for $\beta>0$).

  * (ii)

reflecting boundary condition: $f^{+}(l+)=0$ (or $f^{-}(r-)=0$);

  * (iii)

killing boundary condition: $f(l+)=0$ (or $f(r-)=0$).

In this paper we do not consider sticky boundary conditions or traps (in the
terminology of [2][Section II.7]) – these are boundary conditions involving
the second derivative at a point.

Let $\lambda>0$ and denote by $\psi_{\lambda}(x)$ and $\varphi_{\lambda}(x)$
the fundamental solutions to the second order ODE ${\mathcal{L}}_{X}f=\lambda
f$. We recall that the fundamental solutions are unique (up to multiplication
by a constant) positive solutions to ${\mathcal{L}}_{X}f=\lambda f$ such that
$\psi_{\lambda}$ is increasing and $\varphi_{\lambda}$ is decreasing (see
[2][Section II.10]). Also, if $l$ is a non-singular boundary point, then
$\psi_{\lambda}$ must satisfy one of the boundary conditions stated above (and
the same applies for $\varphi_{\lambda}$ at the right boundary point $r$).
Their Wronskian is defined as follows

${\textnormal{Wr}}[\varphi_{\lambda};\psi_{\lambda}]:=\varphi_{\lambda}(x)\psi^{\prime}_{\lambda}(x)-\varphi^{\prime}_{\lambda}(x)\psi_{\lambda}(x),$
(13)

and it is known that the quantity
$\omega_{\lambda}:={\textnormal{Wr}}[\varphi_{\lambda};\psi_{\lambda}]/s^{\prime}(x)$
is independent of $x$. Let $p^{X}_{t}(x,y)$ be the transition probability
density of the process $X$ with respect to the speed measure, that is

${\mathbb{P}}_{x}(X_{t}\in A)=\int_{A}p^{X}_{t}(x,y){\textnormal{d}}m(y),$

for all Borel sets $A\subset(l,r)$ The function $p^{X}_{t}(x,y)$ is symmetric
in $x$ and $y$, (see [2][Section II.4]). The Green function of the process $X$
has the following expression in terms of the fundamental solutions
$\psi_{\lambda}$ and $\varphi_{\lambda}$:

$G^{X}_{\lambda}(x,y):=\int_{0}^{\infty}e^{-\lambda
t}p^{X}_{t}(x,y){\textnormal{d}}y=\begin{cases}\omega_{\lambda}^{-1}\psi_{\lambda}(x)\varphi_{\lambda}(y),\;\;\;&{\textnormal{
if }}x\leq y,\\\
\omega_{\lambda}^{-1}\psi_{\lambda}(y)\varphi_{\lambda}(x),\;\;\;&{\textnormal{
if }}y\leq x,\end{cases}$ (14)

as stated in [2][Section II.11].

We also mention the following fact: if a process $X=\\{X_{t}\\}_{t\geq 0}$ is
mapped into a process $Y=\\{y(X_{t})\\}_{t\geq 0}$ via an increasing
differentiable function $y=y(x)$, then the scale function and speed measure
are transformed as follows:

$s_{Y}(y(x))=s_{X}(x),\;\;\;m_{Y}^{\prime}(y(x))=m_{X}^{\prime}(x)/y^{\prime}(x).$

The state space for process $Y$ is $(y(l),y(r))$, and the boundary behavior of
$Y$ at points $y(l)$ and $y(r)$ is the same as the boundary behavior of $X$ at
$l$ and $r$.

##  3 Siegmund duality

Let $Z$ and $\widetilde{Z}$ be two Markov processes on $[0,\infty)$. We say
that $\widetilde{Z}$ is a Siegmund dual of $Z$ if

${\mathbb{P}}_{y}(Z_{t}\geq x)={\mathbb{P}}_{x}(\widetilde{Z}_{t}\leq y),$
(15)

for all $x,y,t\geq 0$. This concept was introduced by Siegmund [23] in 1976
and was later studied in [3, 5, 15]. Diffusion processes satisfying (15) were
called conjugate diffusions in [25]. The existence of a dual process
$\widetilde{Z}$ requires the process $Z$ to be stochastically monotone, see
[3, 23]. The latter property can be stated as follows: for any $t>0$ and $x>0$
the function $y\mapsto{\mathbb{P}}_{y}(Z_{t}\geq x)$ is non-decreasing.

Siegmund duality is applied when one of the processes has an absorbing
boundary. For example, assuming that $Z$ is conservative and taking the limit
of (15) as $x\to 0^{+}$, we see that the process $\widetilde{Z}$ must have an
absorbing boundary at zero. Thus, in this case, the distribution of
$\widetilde{Z}_{t}$ will have an atom at zero, and assuming that this
distribution has a density on $(0,\infty)$, it must be given by

$p_{t}^{\widetilde{Z}}(x,y)=\partial_{y}{\mathbb{P}}_{y}(Z_{t}\geq x),$ (16)

which follows from (15) by taking derivative in $y$. Our goal in this section
is to construct a version of Siegmund duality, which works for processes with
reflecting/killing boundary conditions and where we would still have a result
similar to (16).

For the rest of this section, we assume that $X$ is a conservative diffusion
process on the interval $(l,r)\subseteq{\mathbb{R}}$ having the infinitesimal
generator

${\mathcal{L}}_{X}=\frac{1}{2}\sigma^{2}(x)\partial_{x}^{2}+b(x)\partial_{x}=\frac{{\textnormal{d}}}{{\textnormal{d}}m(x)}\frac{{\textnormal{d}}}{{\textnormal{d}}s(x)},$
(17)

where the scale function $s(x)$ and the speed measure ${\textnormal{d}}m(x)$
are determined via (11). We assume that $\sigma(x)$ is a positive and
continuously differentiable function of $x\in(l,r)$. If some boundary points
are non-singular for $X$, we impose a reflecting boundary condition at that
point.

We want to construct a new diffusion process $\widetilde{X}$ on $(l,r)$, which
has the speed measure/scale function equal to the scale function/speed measure
of $X$ (see Section 4 in [5]). In other words, the infinitesimal generator of
$\widetilde{X}$ should have the form

${\mathcal{L}}_{\widetilde{X}}=\frac{{\textnormal{d}}}{{\textnormal{d}}s(x)}\frac{{\textnormal{d}}}{{\textnormal{d}}m(x)}=\frac{1}{2}\sigma^{2}(x)\partial_{x}^{2}+\big{(}\sigma(x)\sigma^{\prime}(x)-b(x)\big{)}\partial_{x}.$
(18)

Note that the operators ${\mathcal{L}}_{X}$ and
${\mathcal{L}}_{\widetilde{X}}$ satisfy the following intertwining relation

$\frac{{\textnormal{d}}}{{\textnormal{d}}s(x)}{\mathcal{L}}_{X}={\mathcal{L}}_{\widetilde{X}}\;\frac{{\textnormal{d}}}{{\textnormal{d}}s(x)},$
(19)

which will play an important role later.

To complete our construction of the process $\widetilde{X}$, we need to
discuss boundary behavior at $l$ and $r$. The type of the boundary behavior
depends on the functions $R$ and $Q$ defined via (12), and when we compute
their analogues $\widetilde{R}$ and $\widetilde{Q}$ for the operator
${\mathcal{L}}_{\widetilde{X}}$ we find that $\widetilde{Q}\equiv R$ and
$\widetilde{R}\equiv Q$ (because the killing measure is zero and
$m_{X}(s)=s_{\widetilde{X}}(s)$ and $s_{X}(x)=m_{\widetilde{X}}(x)$). Note
that $l$ can not be exit-not-entrance for $X$ since we assumed that $X$ is
conservative. From the boundary classification presented on page 12, it
follows that

  * (i)

if $l$ is a non-singular reflecting boundary for $X$, then it is also a non-
singular boundary for $\widetilde{X}$, and in this case we impose a killing
boundary condition $\tilde{f}(l+)=0$ for $\widetilde{X}$;

  * (ii)

if $l$ is a natural boundary for $X$, then it stays a natural boundary for
$\widetilde{X}$;

  * (iii)

if $l$ is an entrance-not-exit, it becomes an exit-not-entrance boundary for
$\widetilde{X}$.

The same applies to the right boundary point $r$. Note that in cases (ii) and
(iii) we do not need to impose boundary conditions (as they are already
determined by the character of the boundary point, see [2][Section II.8]).
This discussion of the boundary behavior completes our construction of the
process $\widetilde{X}$.

The following is our main result in this section, which connects the
transition probabilities of $X$ and $\widetilde{X}$.

######  Theorem 1.

For $t>0$ and $x_{1},x_{2},y\in(l,r)$ such that $x_{1}<x_{2}$ we have

${\mathbb{P}}_{x_{1}}(X_{t}\leq y)={\mathbb{P}}_{x_{2}}(X_{t}\leq
y)+{\mathbb{P}}_{y}(x_{1}<\widetilde{X}_{t}\leq x_{2}).$ (20)

###### Proof.

We will first establish this result in the special case when $b(x)\equiv 0$.
In this case the process $X$ is in the natural scale, that is $s_{X}(x)\equiv
x$. Let $\lambda>0$ and $\psi_{\lambda}$ and $\varphi_{\lambda}$ be the
fundamental solutions to ${\mathcal{L}}_{X}f=\lambda f$. We claim that the
fundamental solutions to ${\mathcal{L}}_{\widetilde{X}}f=\lambda f$ are given
by

$\tilde{\psi}_{\lambda}(x)=\psi^{\prime}_{\lambda}(x),\;\;\;\tilde{\varphi}_{\lambda}(x)=-\varphi^{\prime}_{\lambda}(x).$
(21)

To prove this, we note first that if $f$ is a solution to
${\mathcal{L}}_{X}f=\lambda f$, then $\tilde{f}=f^{\prime}$ is a solution to
${\mathcal{L}}_{\widetilde{X}}\tilde{f}=\lambda\tilde{f}$, as follows from the
intertwining relation (19). It is clear that the functions
$\tilde{\psi}_{\lambda}$ and $\tilde{\varphi}_{\lambda}$ are positive on
$(l,r)$ (since $\psi_{\lambda}$ is increasing and $\varphi_{\lambda}$ is
decreasing). Next we note that $\lambda>0$ and $\psi_{\lambda}$ is a positive
solution to $\tfrac{1}{2}\sigma^{2}(x)f^{\prime\prime}(x)=\lambda f(x)$ ,
which implies that $\psi^{\prime\prime}_{\lambda}(x)>0$, so that
$\tilde{\psi}_{\lambda}=\psi^{\prime}_{\lambda}$ is increasing. The same
argument shows that $\tilde{\varphi}_{\lambda}$ is decreasing. If $l$ is a
natural or entrance-not-exit boundary for $X$, it becomes a natural or exit-
not-entrance boundary for $\widetilde{X}$. In these cases we do not need to
impose a boundary condition for $\psi_{\lambda}$ at $l$ (the increasing
solution to $\widetilde{\mathcal{L}}f=\lambda f$ is unique, up to a
multiplicative constant, as follows from see [2][Section II.10]). If $l$ is a
non-singular reflecting boundary for $X$, then necessarily
$\psi^{\prime}_{\lambda}(l+)=0$, which implies the correct killing boundary
condition $\tilde{\psi}_{\lambda}(l+)=0$. The same considerations apply to the
boundary condition of $\varphi_{\lambda}$ at $r$. This ends the proof of (21).

Next, we check that if
$\omega_{\lambda}={\textnormal{Wr}}[\varphi_{\lambda};\psi_{\lambda}]$ then

${\textnormal{Wr}}[\tilde{\varphi}_{\lambda};\tilde{\psi}_{\lambda}]=-\varphi^{\prime}_{\lambda}(x)\psi^{\prime\prime}_{\lambda}(x)+\varphi^{\prime\prime}_{\lambda}(x)\psi^{\prime}_{\lambda}(x)=\frac{2\lambda}{\sigma^{2}(x)}{\textnormal{Wr}}[\psi_{\lambda};\phi_{\lambda}]=\lambda\omega_{\lambda}s_{\widetilde{X}}^{\prime}(x),$
(22)

where in the second step we used the fact that $\psi_{\lambda}$ and
$\varphi_{\lambda}$ are solutions to
$\tfrac{1}{2}\sigma^{2}(x)f^{\prime\prime}(x)=\lambda f(x)$. Using (14) and
(22) we can write the Green function of the process $\widetilde{X}$ as

$G^{\widetilde{X}}_{\lambda}(x,y)=\begin{cases}-(\lambda\omega_{\lambda})^{-1}\psi^{\prime}_{\lambda}(x)\varphi^{\prime}_{\lambda}(y),\;\;\;&{\textnormal{
if }}x\leq y,\\\
-(\lambda\omega_{\lambda})^{-1}\psi^{\prime}_{\lambda}(y)\tilde{\varphi}^{\prime}_{\lambda}(x),\;\;\;&{\textnormal{
if }}y\leq x.\end{cases}$ (23)

Now we have all the ingredients to complete the proof of (20). We rewrite (20)
in the following form:

$\int_{l}^{y}p^{X}_{t}(x_{1},u){\textnormal{d}}u=\int_{l}^{y}p^{X}_{t}(x_{2},u){\textnormal{d}}u+\int_{x_{1}}^{x_{2}}p^{\widetilde{X}}_{t}(y,v){\textnormal{d}}v.$
(24)

We need to prove that (24) holds for all $t>0$ and $x_{1},x_{2},y$ in $(l,r)$
with $x_{1}<x_{2}$. We take $\lambda>0$, multiply both sides of (29) by
$\exp(-\lambda t)$ and integrate in $t$ over $(0,\infty)$. After applying
Fubini’s theorem, we obtain

$2\int_{l}^{y}G^{X}_{\lambda}(x_{1},u)\sigma^{-2}(u){\textnormal{d}}u=2\int_{l}^{y}G^{X}_{\lambda}(x_{2},u)\sigma^{-2}(u){\textnormal{d}}u+\int_{x_{1}}^{x_{2}}G^{\widetilde{X}}_{\lambda}(y,v){\textnormal{d}}v.$
(25)

By the uniqueness of Laplace transform, to prove (24), it is enough to show
that (25) holds for all $\lambda>0$ and $x_{1},x_{2},y\in(l,r)$ such that
$x_{1}<x_{2}$. Identity (25) is true if

$2\partial_{x}\int_{l}^{y}G^{X}_{\lambda}(x,u)\sigma^{-2}(u){\textnormal{d}}u+G^{\widetilde{X}}_{\lambda}(x,y)=0,$
(26)

for $x,y\in(l,r)$. We will prove first that (26) holds when $l<x<y$. Using
(14) and (23) we rewrite (26) in the form

$\frac{2}{\omega_{\lambda}}\partial_{x}\Big{[}\varphi_{\lambda}(x)\int_{l}^{x}\psi_{\lambda}(u)\sigma^{-2}(u){\textnormal{d}}u+\psi_{\lambda}(x)\int_{x}^{y}\varphi_{\lambda}(u)\sigma^{-2}(u){\textnormal{d}}u\Big{]}-\frac{1}{\lambda\omega_{\lambda}}\psi^{\prime}_{\lambda}(x)\varphi^{\prime}_{\lambda}(y)=0.$
(27)

Since $\psi_{\lambda}$ is a solution to
$\tfrac{1}{2}\sigma^{2}(x)f^{\prime\prime}(x)=\lambda f(x)$, we have

$2\int_{l}^{x}\psi_{\lambda}(u)\sigma^{-2}(u){\textnormal{d}}u=\frac{1}{\lambda}\int_{l}^{x}\psi^{\prime\prime}_{\lambda}(u){\textnormal{d}}u=\frac{1}{\lambda}\psi^{\prime}_{\lambda}(x),$

since $\psi^{\prime}_{\lambda}(l+)=0$ due to boundary conditions on the
process $X$ (see [2][Section II.10]). Similarly

$2\int_{x}^{y}\varphi_{\lambda}(u)\sigma^{-2}(u){\textnormal{d}}u=\frac{1}{\lambda}\big{(}\varphi^{\prime}_{\lambda}(y)-\varphi^{\prime}_{\lambda}(x)\big{)}.$

With the help of the above two equations, we can rewrite (27) as follows

$\frac{1}{\lambda\omega_{\lambda}}\partial_{x}\big{(}{\textnormal{Wr}}[\varphi_{\lambda};\psi_{\lambda}]+\psi_{\lambda}(x)\varphi^{\prime}_{\lambda}(y)\big{)}-\frac{1}{\lambda\omega_{\lambda}}\psi^{\prime}_{\lambda}(x)\varphi^{\prime}_{\lambda}(y)=0,$

and this identity is true since
${\textnormal{Wr}}[\varphi_{\lambda};\psi_{\lambda}]=\omega_{\lambda}$ does
not depend on $x$. Thus, we have established (26) when $l<x<y$. The proof that
(26) holds when $l<y<x$ is done in the same way, and the details are omitted.

This ends the proof of (26) in the driftless case, when $b(x)\equiv 0$. The
general case follows by considering the driftless process
$Y=\\{s(X_{t})\\}_{t\geq 0}$ and constructing the transformed process
$\widetilde{Y}$ and defining $\widetilde{X}_{t}=s^{-1}(\widetilde{Y}_{t})$.
The fact that identity (26) holds for processes $Y$ and $\widetilde{Y}$
implies that it also holds for $X$ and $\widetilde{X}$. We leave it as an
exercise to check that the process $\widetilde{X}$ thus constructed has the
infinitesimal generator as in (18). $\sqcap\kern-8.0pt\hbox{$\sqcup$}$

Condition (20) is equivalent to saying that for any fixed $t>0$ and
$y\in(l,r)$ the function

$x\mapsto{\mathbb{P}}_{x}(X_{t}\leq y)-{\mathbb{P}}_{y}(\widetilde{X}_{t}>x)$
(28)

is constant (does not depend on $x$). Taking derivatives in $x$ of this
constant function (28) and relabelling variables $x\leftrightarrow y$ we
obtain the following result:

######  Corollary 1.

The transition probability densities of $X$ and $\widetilde{X}$ satisfy

$p^{\widetilde{X}}_{t}(x,y)=\partial_{y}{\mathbb{P}}_{y}(X_{t}>x)=\partial_{y}\int_{x}^{r}p^{X}_{t}(y,u){\textnormal{d}}u.$
(29)

##  4 Darboux transform of killed Browian motion

We start with a diffusion process $Y$ on the interval
$(l,r)\subseteq{\mathbb{R}}$, which is a Brownian motion killed at rate
$c(y)$, where $c(y)$ is a positive continuous function for $y\in(l,r)$. The
infinitesimal generator is
${\mathcal{L}}_{Y}=\frac{1}{2}\partial_{y}^{2}-c(y)$. We take the speed
measure to be the Lebesgue measure ${\textnormal{d}}m(x)={\textnormal{d}}x$,
in this case the scale function becomes $s(x)=2x$ (see (11)). If $l=-\infty$,
then $l$ is a natural boundary. If $l$ is finite, then the boundary
classification conditions $R\in L_{1}((l,z))$ and $Q\in L_{1}((l,z))$ from 12
are equivalent to

$R\in
L_{1}((l,z))\Leftrightarrow\int_{l}^{z}c(y)(y-l){\textnormal{d}}y<\infty\;{\textnormal{
and }}\;Q\in
L_{1}((l,z))\Leftrightarrow\int_{l}^{z}c(y){\textnormal{d}}y<\infty$

Thus, if a boundary is an entrance, it must also be an exit, and we have the
following three possibilities for boundary behavior of the process $Y$:
natural, exit-not-entrance and non-singular (entrance-not-exit boundary is not
possible for a killed Brownian motion). We will only consider killed Brownian
motion processes that have one of elastic/reflecting/killing boundary
conditions at each non-singular boundary point (see page 12).

Next, we assume that for some $\lambda\in{\mathbb{R}}$ we have found a
positive function $h(y)$ that is $\lambda$-invariant for the process $Y$, that
is

${\mathbb{E}}_{x}[h(Y_{t})]=e^{\lambda t}h(x),\;\;\;x\in(l,r),\;t>0.$ (30)

This $\lambda$-invariant function $h$ must satisfy the following properties:

(i) $h$ is a solution to ${\mathcal{L}}_{Y}h=\lambda h$ on $(l,r)$; (31) (ii)
$h$ satisfies appropriate boundary conditions at each non-singular boundary
point. (32)

The first property follows from the fact that the process $\exp(-\lambda
t)h(Y_{t})$ is a martingale. To check the second property, we take
$\mu>\max(0,\lambda)$, multiply both sides of (30) by $e^{-\mu t}$, integrate
in $t\in(0,\infty)$ and obtain

$(\mu-\lambda)\int_{l}^{r}G^{Y}_{\mu}(x,y)h(y){\textnormal{d}}y=h(x),$ (33)

where $G^{Y}_{\mu}(x,y)$ is the Green function of the process $Y$. Thus, if
$h$ is bounded on $(l,r)$ (which happens, for example, if both boundaries are
non-singular), it must be in the domain of the infinitesimal generator
${\mathcal{L}}_{Y}$ (because it lies in the image of the resolvent operator).
Therefore, $h$ must satisfy the appropriate boundary conditions at each non-
singular boundary (see [2][Section II.7]). The proof of (32) in the general
case, when $h$ can be unbounded, is given in Appendix B.

In what follows, we will require one additional assumption on $h$:

${\mathfrak{m}}_{h}:=\sup\limits_{y\in(l,r)}\Big{[}c(y)-\Big{(}\frac{h^{\prime}(y)}{h(y)}\Big{)}^{2}\Big{]}<\infty.$
(34)

Remark 1. Note that ${\mathfrak{m}}_{h}$ is finite when $c(y)$ is bounded on
$(l,r)$. In all our examples in Section 5, we actually have
${\mathfrak{m}}_{h}=0$, though it is easy to construct examples with
${\mathfrak{m}}_{h}$ equal to any positive number (one could take Brownian
motion on ${\mathbb{R}}$ killed at a positive constant rate and a
$\lambda$-invariant function $h(y)=\cosh(y)$, see Section 5.1). It could be
true that ${\mathfrak{m}}_{h}$ is always finite, but we could not prove this
or find a counterexample.

As the first step of constructing the Darboux transform of $Y$, we introduce a
new Markov process $X$ via Doob’s $h$-transform:

${\mathbb{P}}_{x}(X_{t}\in A)=\frac{e^{-\lambda
t}}{h(x)}{\mathbb{E}}_{x}[{\mathbf{1}}_{\\{Y_{t}\in A\\}}h(Y_{t})].$ (35)

The process $X$ is a conservative Markov process, since $h$ is
$\lambda$-invariant for $Y$. Moreover, it is a diffusion process on $(l,r)$
with the infinitesimal generator given by

${\mathcal{L}}_{X}=\frac{1}{h}{\mathcal{L}}_{Y}h-\lambda=\frac{1}{2}\partial_{y}^{2}+\frac{h^{\prime}(y)}{h(y)}\partial_{y}.$
(36)

The derivatives of the scale function and the speed measure of the process $X$
are $s^{\prime}(x)=2h^{-2}(x)$ and $m^{\prime}(x)=h^{2}(x)$ (this follows from
(11)), thus the functions $R$ and $Q$ in (12) are given by

$R_{X}(x)=2h^{-2}(x)\int_{x}^{z}h^{2}(y){\textnormal{d}}y,\;\;\;Q_{X}(x)=2h^{2}(x)\int_{x}^{z}h^{-2}(y){\textnormal{d}}y.$
(37)

The above two functions $R_{X}$ and $Q_{X}$ can be used to find the type of
the boundary points $l$ and $r$ for the process $X$ (see Section 2).

For the rest of this paper we assume that if $l$ or $r$ is a non-singular
boundary point for the process $X$, then it is a reflecting boundary. The
following proposition will be useful for our examples in Section 5.

######  Proposition 1.

If a non-singular boundary point of $Y$ has a reflecting or elastic boundary
condition, then this boundary point is also non-singular for $X$ with a
reflecting boundary condition.

###### Proof.

The Green function of $X$ (with respect to the speed measure
${\textnormal{d}}m(x)=h^{2}(x){\textnormal{d}}x$) is given by

$G^{X}_{\mu}(x,y)=\frac{G^{Y}_{\lambda+\mu}(x,y)}{h(x)h(y)},$

where $\mu>\max(0,-\lambda)$ and $x,y\in(l,r)$. Comparing the above formula
and (14), we see that the fundamental increasing and decreasing solutions
$\psi^{X}_{\mu}$ and $\varphi^{X}_{\mu}$ for the operator ${\mathcal{L}}_{X}$
must be given by

$\psi^{X}_{\mu}(x)=\frac{\psi^{Y}_{\lambda+\mu}(x)}{h(x)},\;\;\;\varphi^{X}_{\mu}(x)=\frac{\varphi^{Y}_{\lambda+\mu}(x)}{h(x)}.$

Assume that $l$ is a non-singular boundary point for $Y$ with reflecting or
elastic boundary condition, which can be written as $g^{\prime}(l+)=\gamma
g(l+)$ for some $\gamma\geq 0$ (here we used the fact that $Y$ is in the
natural scale, so that $g^{+}(y)=g^{\prime}(y)$). The function $h$, being
$\lambda$-invariant for $Y$, must satisfy equation ${\mathcal{L}}_{Y}h=\lambda
h$ and the corresponding boundary condition at $l$. The fundamental solution
$\psi^{Y}_{\lambda+\mu}(x)$ must satisfy the same boundary condition at $l$.
From here, it is easy to check that

$\lim\limits_{x\to
l+}h^{2}(x)\frac{{\textnormal{d}}}{{\textnormal{d}}x}\psi^{X}_{\mu}(x)=\frac{1}{2}\lim\limits_{x\to
l+}\Big{[}h(x)\frac{{\textnormal{d}}}{{\textnormal{d}}x}\psi^{Y}_{\lambda+\mu}(x)-\psi^{Y}_{\lambda+\mu}(x)\frac{{\textnormal{d}}}{{\textnormal{d}}x}h(x)\Big{]}=0.$

Thus, for $\mu>\max(0,-\lambda)$ the function $f(x)=\psi^{X}_{\mu}(x)$
satisfies $f^{+}(l+)=0$, which is a reflecting boundary condition. The values
of the Green function $G^{X}_{\mu}(x,y)$ for $\mu>\max(0,-\lambda)$ uniquely
determine the transition probability density of the process $X$ (by the
uniqueness of the Laplace transform), and we conclude that $l$ is indeed a
reflecting boundary for $X$.

When $r$ is a non-singular boundary point for $Y$, the proof proceeds along
the same lines, except that now we focus on the boundary behavior of the
decreasing fundamental solution $\varphi^{X}_{\mu}(x)$.
$\sqcap\kern-8.0pt\hbox{$\sqcup$}$

The process $X$ thus constructed is a conservative diffusion process with
reflecting boundary condition at every non-singular boundary point. We
construct its Siegmund transform $\widetilde{X}$ as described in Section 3.
According to (18), the infinitesimal generator of $\widetilde{X}$ is

${\mathcal{L}}_{\widetilde{X}}=\frac{1}{2}\partial_{y}^{2}-\frac{h^{\prime}(y)}{h(y)}\partial_{y}.$
(38)

The process $\widetilde{X}$ can have only natural, exit-not-entrance, or non-
singular killing boundaries. Now we want to turn it into a killed Brownian
motion by “removing the drift term”. We achieve this through another Doob’s
$h$-transform. First we establish that the function $h$ is $\nu$-excessive for
the process $\widetilde{X}$ with $\nu={\mathfrak{m}}_{h}+\lambda$.

######  Proposition 2.

For $x\in(l,r)$ and $t>0$ we have ${\mathbb{E}}_{x}[h(\widetilde{X}_{t})]\leq
e^{({\mathfrak{m}}_{h}+\lambda)t}h(x)$.

###### Proof.

The process $\widetilde{X}$ is a diffusion process solving the stochastic
differential equation

${\textnormal{d}}\widetilde{X}_{t}=-\frac{h^{\prime}(\widetilde{X}_{t})}{h(\widetilde{X}_{t})}{\textnormal{d}}t+{\textnormal{d}}W_{t},\;\;\;t<\zeta,$

where $\zeta$ is the first time the process reaches the boundary of $(l,r)$.
When both boundaries are natural, we have $\zeta=+\infty$, and if one of the
boundaries is exit-not-entrance or non-singular killing boundary, then the
process is killed when it hits that boundary (in other words, the process at
time $\zeta$ is sent to a cemetery state $\widetilde{X}_{\zeta}=\Delta$). We
apply Ito’s formula to $\ln h(\widetilde{X}_{t})$:

$\ln h(\widetilde{X}_{t})=\ln
h(x)+\int_{0}^{t}\frac{h^{\prime}(\widetilde{X}_{s})}{h(\widetilde{X}_{s})}{\textnormal{d}}\widetilde{X}_{t}+\frac{1}{2}\int_{0}^{t}v(\widetilde{X}_{s}){\textnormal{d}}s,\;\;\;t<\zeta,$
(39)

where

$v(x):=\frac{{\textnormal{d}}^{2}}{{\textnormal{d}}x^{2}}\ln
h(x)=\frac{h^{\prime\prime}(x)}{h(x)}-\Big{(}\frac{h^{\prime}(x)}{h(x)}\Big{)}^{2}=2(c(x)+\lambda)-\Big{(}\frac{h^{\prime}(x)}{h(x)}\Big{)}^{2}.$

In the last step we used the fact that $h$ is a solution to
${\mathcal{L}}_{Y}h=\lambda h$. Let us denote

$U_{t}:=\frac{h^{\prime}(\widetilde{X}_{t})}{h(\widetilde{X}_{t})},\;\;\;Z_{t}:=\int_{0}^{t}U_{s}\;{\textnormal{d}}W_{s},\;\;\;t<\zeta.$

After rearranging the terms in (39) and using the fact that
$h(\widetilde{X}_{t})=h(\Delta)=0$ for $t\geq\zeta$ we obtain

$h(\widetilde{X}_{t})\leq h(x)\times e^{Z_{t}-\frac{1}{2}\langle
Z\rangle_{t}+\frac{1}{2}\int_{0}^{t}(v(\widetilde{X}_{s})-U_{s}^{2}){\textnormal{d}}s},$
(40)

which is now valid for all $t\geq 0$.

Next, we observe the following two facts. The process
$\exp\big{(}Z_{t}-(1/2)\langle Z\rangle_{t}\big{)}$ is a positive local
martingale; thus, it is a supermartingale, and its expected value is less than
one. We also have an upper bound

$\frac{1}{2}\int_{0}^{t}(v(\widetilde{X}_{s})-U_{s}^{2}){\textnormal{d}}s=\int_{0}^{t}\Big{[}c(\widetilde{X}_{s})+\lambda-\Big{(}\frac{h^{\prime}(\widetilde{X}_{s})}{h(\widetilde{X}_{s})}\Big{)}^{2}\Big{]}{\textnormal{d}}s\leq({\mathfrak{m}}_{h}+\lambda)t,$

which follows from (34). Using these two facts and taking expected values of
both sides of (40) gives us the desired result.
$\sqcap\kern-8.0pt\hbox{$\sqcup$}$

Now that we established the fact that $h$ is $\nu$-excessive for the process
$\widetilde{X}_{t}$ (with $\nu={\mathfrak{m}}_{h}+\lambda$), we can define a
new process $\widetilde{Y}$ via Doob’s $h$-transform

${\mathbb{P}}_{y}(\widetilde{Y}_{t}\in
A)=\frac{e^{-({\mathfrak{m}}_{h}+\lambda)t}}{h(y)}{\mathbb{E}}_{y}[{\mathbf{1}}_{\\{\widetilde{X}_{t}\in
A\\}}h(\widetilde{X}_{t})].$ (41)

The infinitesimal generator of $\widetilde{Y}$ is

${\mathcal{L}}_{\widetilde{Y}}=\frac{1}{h}{\mathcal{L}}_{\widetilde{X}}h-({\mathfrak{m}}_{h}+\lambda)=\frac{1}{2}\partial_{y}^{2}-\tilde{c}(y),$
(42)

where

$\tilde{c}(y)={\mathfrak{m}}_{h}+\Big{(}\frac{h^{\prime}(y)}{h(y)}\Big{)}^{2}-c(y).$
(43)

Note that $\tilde{c}(y)$ is a positive continuous function on $(l,r)$. This is
due to the way we defined ${\mathfrak{m}}_{h}$ in (34) and due to our
assumption that ${\mathfrak{m}}_{h}<\infty$. Thus we can identify
$\widetilde{Y}$ as a Brownian motion killed at rate $\tilde{c}(y)$. We call
the process $\widetilde{Y}$ the Darboux transform of the killed Brownian
motion process $Y$. The type of boundary points of $\widetilde{Y}$ can be
determined via the killing term $\tilde{c}(y)$ and integrability conditions
(12). We also know that if $Y$ has a non-singular boundary point with a
reflecting or elastic boundary condition, then this point is also non-singular
for $\widetilde{Y}$ with a killing boundary condition. This follows from the
fact that $X$ has a reflecting boundary condition at this non-singular
boundary point (see Proposition 1) and its Siegmund dual $\widetilde{X}$ has
killing boundary condition, which is preserved under Doob’s $h$-transform
since $h(l+)>0$.

process $Y$ |  | process $\widetilde{Y}$
---|---|---
${\mathcal{L}}_{Y}=\frac{1}{2}\partial_{y}^{2}-c(y)$ |  | ${\mathcal{L}}_{\widetilde{Y}}=\frac{1}{2}\partial_{y}^{2}-\tilde{c}(y)$
$\lambda$-invariant function $h(y)$ |  | $\tilde{c}(y)={\mathfrak{m}}_{h}+\big{(}\frac{h^{\prime}(y)}{h(y)}\big{)}^{2}-c(y)$
$\downarrow$ |  | $\uparrow$
Doob’s $h$-transform of $Y$ |  | Doob’s $h$-transform of $\widetilde{X}$
$\downarrow$ |  | $\uparrow$
process $X$ |  $\longrightarrow\;\;\;\;$ Siegmund dual of $X$ $\;\;\;\;\longrightarrow$ | process $\widetilde{X}$
${\mathcal{L}}_{X}=\frac{1}{2}\partial_{y}^{2}+\frac{h^{\prime}(y)}{h(y)}\partial_{y}$ |  | ${\mathcal{L}}_{\widetilde{X}}=\frac{1}{2}\partial_{y}^{2}-\frac{h^{\prime}(y)}{h(y)}\partial_{y}$
Table 1: The three steps in constructing Darboux transform of killed Brownian
motion process $Y$.

We summarize the steps in our construction of the Darboux transformed process
$\widetilde{Y}$ in Table 1. We want to emphasize that this construction
depends on the following two assumptions:

(i) $h$ is $\lambda$-invariant function for $Y$ for which ${\mathfrak{m}}_{h}$
(defined by (34)) is finite; (44) (ii) $X$ has a reflecting boundary condition
at each non-singular boundary point. (45)

In practical applications, we find a $\lambda$-invariant function $h$ by
solving the equation ${\mathcal{L}}_{Y}h=\lambda h$ with appropriate boundary
conditions at each non-singular boundary. These two conditions alone do not
guarantee that $h$ is $\lambda$-invariant (see [8, 19]); thus, to verify that
this candidate function $h$ is indeed $\lambda$-invariant, we check that (30)
holds by actually computing the left-hand side (an alternative way is to apply
[8][Theorem 2.7] or [19][Corollary 2.2]). The condition
${\mathfrak{m}}_{h}<\infty$ is easy to verify directly (once we have
expressions for $c(x)$ and $h(x)$). The second assumption (45) is covered by
Proposition 1 in some cases, and we conjecture that this result is, in fact,
always true. The intuitive argument for validity of assumption (45) is the
following: the process $X$ is conservative by construction, thus it can not
have elastic or killing boundary conditions at non-singular boundary (as this
would imply ${\mathbb{P}}_{x}(X_{t}\in(l,r))<1$). It can not have sticky
boundaries or traps, as then the distribution of $X_{t}$ would have an atom at
that boundary (which is impossible since the law of $Y_{t}$ has no atoms).
Thus, the only possible boundary conditions are reflecting or non-local
boundary conditions (see Theorem 2 on page 39 in [17]). The latter can be
ruled out if we could show that $X$ has continuous paths, but we did not
pursue this further, as Proposition 1 was sufficient for our examples.

Next, we present our main result in this section, which connects transition
probability densities of the process $Y$ and its Darboux transform
$\widetilde{Y}$.

######  Theorem 2.

Assume that the two conditions (44) and (45) are satisfied. Let
$\widetilde{Y}$ be the Darboux transform of $Y$. Then

$p^{\widetilde{Y}}_{t}(x,y)=e^{-({\mathfrak{m}}_{h}+2\lambda)t}\frac{h(x)}{h(y)}\partial_{x}\Big{[}\frac{1}{h(x)}\int_{y}^{r}p^{Y}_{t}(x,u)h(u){\textnormal{d}}u\Big{]}.$
(46)

###### Proof.

Since $X$ (respectively, $\widetilde{Y}$) is the Doob’s $h$-transform of $Y$
(respectively, $\widetilde{X}$), their transition probability densities
satisfy

$\displaystyle p^{X}_{t}(x,y)=\frac{e^{-\lambda
t}}{h(x)}p_{t}^{Y}(x,y)h(y),\;\;\;p^{\widetilde{Y}}_{t}(x,y)=\frac{e^{-({\mathfrak{m}}_{h}+\lambda)t}}{h(x)}p_{t}^{\widetilde{X}}(x,y)h(y).$

The transition probability densities of $\widetilde{X}$ and $X$ are related
via (29). Combining these identities and using the symmetry of
$p_{t}^{\widetilde{Y}}(x,y)$ with respect to $x$ and $y$, we obtain (46).
$\sqcap\kern-8.0pt\hbox{$\sqcup$}$

Next, we present an example demonstrating the importance of our assumption
(44) of $\lambda$-invariance of $h$. As we discussed in the introduction, to
define the Darboux transform of operator ${\mathcal{L}}$ we need a positive
function $h$ that satisfies ${\mathcal{L}}h=\lambda h$. If $h$ is
$\lambda$-invariant for the process $Y$, that is
${\mathbb{E}}_{y}[h(Y_{t})]=\exp(\lambda t)h(y)$, then $h$ necessarily
satisfies ${\mathcal{L}}_{Y}h=\lambda h$. However, the converse is not true.
For example, consider a process $Y$, which is a Brownian motion on
$(0,\infty)$ reflected at zero. The transition probability of $Y$ (with
respect to Lebesgue measure) is

$p_{t}^{Y}(x,y)=\frac{1}{\sqrt{2\pi
t}}\Big{(}e^{-\frac{1}{2t}(y-x)^{2}}+e^{-\frac{1}{2t}(y+x)^{2}}\Big{)},\;\;\;t,x,y>0.$

The infinitesimal generator is
${\mathcal{L}}_{Y}=\frac{1}{2}\partial_{y}^{2}$. The point $l=0$ is a non-
singular boundary and we have a boundary condition $f^{\prime}(0+)=0$.
Consider a function $h(y)=y$. This function is a solution to
${\mathcal{L}}_{Y}h=0$, but it is not an invariant function for $Y$, since $h$
does not satisfy the reflecting boundary condition at zero. If we take the
above expression for $p^{Y}_{t}(x,y)$ and compute the expression in the right-
hand side of (46) (with $h(y)=y$), we would obtain

$\frac{h(x)}{h(y)}\partial_{x}\Big{[}\frac{1}{h(x)}\int_{y}^{r}p^{Y}_{t}(x,u)h(u){\textnormal{d}}u\Big{]}=\frac{1}{\sqrt{2\pi
t}}\Big{[}e^{-\frac{1}{2t}(y-x)^{2}}\Big{(}1-\frac{t}{xy}\Big{)}-e^{-\frac{1}{2t}(y+x)^{2}}\Big{(}1+\frac{t}{xy}\Big{)}\Big{]}.$

The expression in the right-hand side is negative for small values of $y$,
thus it can not be a transition probability density. This example confirms the
necessity of assumption (44).

Finally, we would like to point out that the second step in our construction
of the Darboux transformed process (see Table 1) is closely related to the
concept of Krein dual strings, see [4, 16]. To demonstrate this connection, we
introduce two increasing functions $u(x)$ and $v(x)$ via
$u^{\prime}(x)=h^{2}(x)$ and $v^{\prime}(x)=2h^{-2}(x)$ and define diffusion
processes $V_{t}=v(X_{t})$ and $U_{t}=u(\widetilde{X}_{t})$. These are
processes in natural scale (the functions $v$ and $u$ are scale functions for
$X$ and $\widetilde{X}$), and their infinitesimal generators can be written as
follows:

${\mathcal{L}}_{U}=\frac{{\textnormal{d}}}{{\textnormal{d}}M(u)}\frac{{\textnormal{d}}}{{\textnormal{d}}u},\;\;\;{\mathcal{L}}_{V}=\frac{{\textnormal{d}}}{{\textnormal{d}}m(v)}\frac{{\textnormal{d}}}{{\textnormal{d}}v},$

where $M$ and $m$ are speed measures of $U$ and $V$ (see (11)):

$\frac{{\textnormal{d}}m}{{\textnormal{d}}v}=\frac{1}{2}h^{4}(x(v)),\;\;\;\frac{{\textnormal{d}}M}{{\textnormal{d}}u}=2h^{-4}(x(u)).$

Here $x(u)$ and $x(v)$ are the inverse functions of $u(x)$ and $v(x)$. One can
check that

$\frac{{\textnormal{d}}M}{{\textnormal{d}}m}\frac{{\textnormal{d}}m}{{\textnormal{d}}v}=1,$

which means that $m(v)$ is the inverse function of $M(u)$. This shows that $m$
is the Krein dual string of $M$, in the terminology of [4, 16]. Thus, the
second step (the Siegmund duality) in Table 1 is closely related to Krein
duality.

##  5 Examples

In this section we present five examples of Darboux transformed processes. In
the first four examples, we take the process $Y$ to be the Brownian motion on
an interval $(l,r)$ with various boundary conditions, and our last example is
related to the Ornstein-Uhlenbeck process. For all killed Brownian motion
processes in this section we take the speed measure
${\textnormal{d}}m(x)={\textnormal{d}}x$, so that the transition probability
density and the Green function for each process are given with respect to the
Lebesgue measure.

Our goal in this section is to compute the transition probability density of
the transformed process $\widetilde{Y}$ and (in most cases) to provide its
spectral expansion. Computing $p_{t}^{\widetilde{Y}}(x,y)$ via (46) is very
time-consuming if done by hand. Instead, we obtained all expressions for
$p_{t}^{\widetilde{Y}}(x,y)$ in this section using symbolic computations, and
then we verified numerically that our formulas were correct. The Matlab
programs for verifying these formulas via symbolic computations can be found
at kuznetsov.mathstats.yorku.ca/code/.

Some computations in this section will require the following simple result
(which is probably well known, but we could not find it in the literature in
this exact form).

######  Lemma 1.

Consider a second order linear differential operator
${\mathcal{L}}=\frac{1}{2}\partial_{y}^{2}-c(y)$, where $c$ is a continuous
function of $y\in(l,r)$. Assume that $h$, $f$ and $g$ are $C^{2}$ functions on
$(l,r)$ such that

  * (i)

$h$ is positive and satisfies ${\mathcal{L}}h=\lambda h$ for $y\in(l,r)$;

  * (ii)

$g$ satisfies ${\mathcal{L}}g=\mu g$ for $y\in(l,r)$.

Denote $\tilde{f}={\mathcal{D}}_{h}f$ and $\tilde{g}=\mathcal{D}_{h}g$. Then

$\int\tilde{f}(y)\tilde{g}(y){\textnormal{d}}y=f(y)\tilde{g}(y)+2(\lambda-\mu)\int
f(y)g(y){\textnormal{d}}y.$ (47)

If $f$ also satisfies ${\mathcal{L}}f=\mu f$ for $y\in(l,r)$, then

${\textnormal{Wr}}[\tilde{f},\tilde{g}]=2(\lambda-\mu){\textnormal{Wr}}[f,g].$
(48)

The proof of (48) is obtained by calculating
${\textnormal{Wr}}[\tilde{f},\tilde{g}]$ and using the fact that
$g^{\prime\prime}=2(c+\mu)g$ and $h^{\prime\prime}=2(c+\lambda)h$. Formula
(47) is derived in a similar way, first by taking derivative in $y$ of both
sides and then simplifying the result using the above identities for
$g^{\prime\prime}$ and $h^{\prime\prime}$.

We also record here the following useful fact (which follows from the
intertwining relation (7) and our construction of Darboux transformed process
$\widetilde{Y}$):

$\;{\mathcal{L}}_{Y}f=(\mu+{\mathfrak{m}}_{h}+2\lambda)f\;\;\;\Longrightarrow\;\;\;{\mathcal{L}}_{\widetilde{Y}}({\mathcal{D}}_{h}f)=\mu({\mathcal{D}}_{h}f).$
(49)

###  5.1 Brownian motion on ${\mathbb{R}}$

We take $Y$ to be a standard Brownian motion on ${\mathbb{R}}$. In this case
we have $l=-\infty$, $r=+\infty$, both boundaries are natural and the
transition probability density of $Y$ is

$p^{Y}_{t}(x,y)=\frac{1}{\sqrt{2\pi
t}}e^{-\frac{1}{2t}(y-x)^{2}},t>0,\;\;x,y\in{\mathbb{R}}.$

The infinitesimal generator is
${\mathcal{L}}_{Y}=\frac{1}{2}\partial_{y}^{2}$, thus the killing term $c(x)$
is zero. We set $\lambda=1/2$ and take $h(x)=\cosh(x)$. Note that $h$ solves
the equation ${\mathcal{L}}_{Y}h=\lambda h$. To check that $h$ is indeed
$\lambda$-invariant for $Y$, we verify that

$\int^{r}_{l}p^{Y}_{t}(x,y)h(y){\textnormal{d}}y=e^{\lambda t}h(x).$ (50)

The integral in the left-hand side can be easily computed explicitly. We also
provide the Matlab code for symbolic verification of this result. Another way
to prove that $h$ is $\lambda$-invariant is via [8][Theorem 2.7] or
[19][Corollary 2.2].

Next, we verify that condition (44) holds with ${\mathfrak{m}}_{h}=0$ and
condition (45) is also satisfied since both boundaries are natural for the
process $X$ with the generator given in (36). Thus, we can construct the
Darboux transformed process $\widetilde{Y}$ and identify it as a Brownian
motion on ${\mathbb{R}}$ killed at rate $\tilde{c}(y)=\tanh(y)^{2}$ (this
latter expression follows from (43)). The transition probability density of
$\widetilde{Y}$ is given by

$p^{\widetilde{Y}}_{t}(x,y)=\frac{1}{\sqrt{2\pi
t}}e^{-t-\frac{1}{2t}(y-x)^{2}}+\frac{e^{-\frac{t}{2}}}{2\cosh(x)\cosh(y)}\Big{[}\Phi\Big{(}\frac{y-x+t}{\sqrt{t}}\Big{)}-\Phi\Big{(}\frac{y-x-t}{\sqrt{t}}\Big{)}\Big{]},$
(51)

where $\Phi(x)$ is the CDF of standard normal distribution. Formula (51) was
obtained from (46) by symbolic computation (and then verified numerically).

Now we turn our attention to the Green function of $\widetilde{Y}$. The
fundamental increasing/decreasing solutions of ${\mathcal{L}}_{Y}f=\mu f$
(with $\mu>0$) are $\psi_{\mu}^{Y}(y)=\exp(\sqrt{2\mu}y)$ and
$\varphi_{\mu}^{Y}(y)=\exp(-\sqrt{2\mu}y)$. Let us denote

$F(z,y)={\mathcal{D}}_{h}e^{zy}=e^{zy}(z-\tanh(y)).$ (52)

The functions

$\psi_{\mu}^{\widetilde{Y}}(y):=F(\sqrt{2(1+\mu)},y),\;\;\;\varphi_{\mu}^{\widetilde{Y}}(y):=F(-\sqrt{2(1+\mu)},y)$

are the increasing/decreasing fundamental solutions to
${\mathcal{L}}_{\widetilde{Y}}f=\mu f$ (this follows from (49)). With the help
of (48) we compute
${\textnormal{Wr}}[\varphi_{\mu}^{\widetilde{Y}};\psi_{\mu}^{\widetilde{Y}}]=-2(1+2\mu)\sqrt{2(1+\mu)}$
and using (14) we find the Green function of $\widetilde{Y}$ in the form

$G^{\widetilde{Y}}_{\mu}(x,y)=-\frac{1}{\sqrt{2(1+\mu)}(1+2\mu)}\times\begin{cases}F(\sqrt{2(1+\mu)},x)F(-\sqrt{2(1+\mu)},y),\;\;&{\textnormal{
if }}\;x<y,\\\ F(-\sqrt{2(1+\mu)},x)F(\sqrt{2(1+\mu)},y),\;\;&{\textnormal{ if
}}\;y<x.\end{cases}$ (53)

We claim that the transition probability density of $\widetilde{Y}$ can be
written in the following spectral representation form:

$p^{\widetilde{Y}}_{t}(x,y)=\frac{e^{-\frac{t}{2}}}{2\cosh(x)\cosh(y)}+\frac{1}{2\pi}\int_{{\mathbb{R}}}e^{-(1+\frac{z^{2}}{2})t}F({\textnormal{i}}z,x)F(-{\textnormal{i}}z,y)\frac{{\textnormal{d}}z}{1+z^{2}}.$
(54)

There are three ways how one can obtain this result. The first method is to
use general Sturm-Liouville theory (see Section 2.2 in [24]). The second is to
write the transition probability density as the inverse Laplace transform of
the Green function in (53)

$p^{\widetilde{Y}}_{t}(x,y)=\frac{1}{2\pi}\int_{c+{\textnormal{i}}{\mathbb{R}}}G^{\widetilde{Y}}_{\mu}(x,y)e^{\mu
t}{\textnormal{d}}\mu,\;\;\;c>0,$

and transform the contour of integration to a Hankel-type contour, which goes
around the interval $(-\infty,-1]$ in the counterclockwise direction (starting
at $-\infty$). While transforming this contour of integration we will collect
a residue at $\mu=-1/2$, which will give us the first term in (54), and the
integral over Hankel’s contour will give the integral term in (54). The third
method is probably the simplest (though the least enlightening): we compute
the integral in the right-hand side of (54) and transform the resulting
expression into the form (51). More details on this last method are provided
in Appendix A.

The formula (54) shows that the effect of Darboux transformation on the
spectral representation of the transition semigroup of the diffusion process
$Y$ is that we have shifted the spectrum by $-1$ and inserted a new eigenvalue
at $-1/2$ (with the corresponding eigenfunction $1/h(x)=1/\cosh(x)$).

###  5.2 Brownian motion on $(0,\infty)$, killed at $0$

Now we take $Y$ to be a Brownian motion on $(0,\infty)$, killed at the first
time it hits $0$. The transition probability density is given by

$p_{t}^{Y}(x,y)=\frac{1}{\sqrt{2\pi
t}}\Big{(}e^{-\frac{1}{2t}(y-x)^{2}}-e^{-\frac{1}{2t}(y+x)^{2}}\Big{)}=\frac{2}{\pi}\int_{0}^{\infty}e^{-\frac{z^{2}}{2}t}\sin(zx)\sin(zy){\textnormal{d}}z,\;t>0,\;\;x,y>0,$

see page 120 in [2]. The infinitesimal generator is
$\widetilde{\mathcal{L}}_{Y}=\frac{1}{2}\partial_{y}^{2}$, the killing term
$c(x)$ is zero. The point $l=0$ is a non-singular boundary for $Y$ and we have
a killing boundary condition $f(0+)=0$. We set $\lambda=1/2$ and find a
function $h(x)=\sinh(x)$ by solving equation ${\mathcal{L}}_{Y}h=\lambda h$
with the boundary condition $h(0+)=0$ and then we verify that (50) holds with
$\lambda=1/2$ (and $l=0$, $r=+\infty$) (or we apply [8][Theorem 2.7] or
[19][Corollary 2.2]). Thus $h$ is $\lambda$-invariant for $Y$.

We check that condition (44) holds with ${\mathfrak{m}}_{h}=0$ and that
condition (45) is also satisfied, since both boundaries are natural for the
process $X$ with the generator given in (36). The process $\widetilde{Y}$ (the
Darboux transform of $Y$) is a killed Brownian motion on $(0,\infty)$ with the
killing rate $\tilde{c}(y)=\coth^{2}(y)$. Both boundaries $0$ and $+\infty$
are natural for $\widetilde{Y}$. The transition probability density of
$\widetilde{Y}$ is

$\displaystyle p_{t}^{\widetilde{Y}}(x,y)$
$\displaystyle=\frac{e^{-t}}{\sqrt{2\pi
t}}\Big{(}e^{-\frac{1}{2t}(y-x)^{2}}+e^{-\frac{1}{2t}(y+x)^{2}}\Big{)}$ (55)
$\displaystyle+\frac{e^{-t/2}}{2\sinh(x)\sinh(y)}\Big{[}\Phi\Big{(}\frac{y-x-t}{\sqrt{t}}\Big{)}+\Phi\Big{(}\frac{y+x+t}{\sqrt{t}}\Big{)}-\Phi\Big{(}\frac{y+x-t}{\sqrt{t}}\Big{)}-\Phi\Big{(}\frac{y-x+t}{\sqrt{t}}\Big{)}\Big{]}.$

The fundamental increasing/decreasing solutions to
${\mathcal{L}}_{\widetilde{Y}}f=\mu f$ are

$\psi^{\widetilde{Y}}_{\mu}(x)=\frac{1}{z^{2}-1}\Big{(}z\cosh(zx)-\sinh(zx)\coth(x)\Big{)},\;\;\;\varphi^{\widetilde{Y}}_{\mu}(x)=e^{-zx}(z+\coth(x)),$

where we denoted $z=\sqrt{2(\mu+1)}$. The Green function for the process
$\tilde{Y}$ is

$G^{\widetilde{Y}}_{\mu}(x,y)=\frac{2}{\sqrt{2(\mu+1)}}\times\begin{cases}\tilde{\psi}^{+}_{\mu}(x)\tilde{\varphi}^{-}_{\mu}(y),\;\;&{\textnormal{
if }}\;x<y,\\\
\tilde{\varphi}^{-}_{\mu}(x)\tilde{\psi}^{+}_{\mu}(y),\;\;&{\textnormal{ if
}}\;y<x.\end{cases}$

The spectral representation for $\tilde{p}^{\widetilde{Y}}_{t}(x,y)$ is given
by

$p^{\widetilde{Y}}_{t}(x,y)=\frac{2}{\pi}\int_{0}^{\infty}e^{-(1+\frac{z^{2}}{2})t}f(z,x)f(z,y)\frac{{\textnormal{d}}z}{1+z^{2}},$
(56)

where we denoted $f(z,y):=z\cos(zx)-\sin(zx)\coth(x)$. In this case the
Darboux transformation has resulted in shifting the spectrum by $-1$, but no
new eigenvalues are created. Intuitively this can be explained by noting that
the function $1/h(x)=1/\sinh(x)$ (which would be the candidate eigenfunction
with the eigenvalue $-1/2$ for the transition semigroup of $\tilde{Y}$) is not
in $L_{2}((0,\infty),{\textnormal{d}}x)$.

###  5.3 Brownian motion on $(0,\infty)$, killed elastically at $0$

Next, we consider the example when $Y$ is a Brownian motion on $(0,\infty)$
killed elastically at zero, see page 125 in [2]. In this case the
infinitesimal generator is ${\mathcal{L}}_{Y}=\frac{1}{2}\partial_{y}^{2}$,
the killing term $c(x)$ is zero. The point $l=0$ is a non-singular boundary
for $Y$ and we have an elastic boundary condition $f^{\prime}(0^{+})=\gamma
f(0+)$, where $\gamma>0$. The transition probability is

$p^{Y}_{t}(x,y)=\frac{1}{\sqrt{2\pi
t}}\Big{[}e^{-\frac{1}{2t}(y-x)^{2}}+e^{-\frac{1}{2t}(y+x)^{2}}\Big{]}-2\gamma
e^{\gamma(x+y)+\frac{\gamma^{2}t}{2}}\Phi\Big{(}-\frac{x+y+\gamma
t}{\sqrt{t}}\Big{)}.$

We set $\lambda=1/2$ and solving the equation ${\mathcal{L}}_{Y}h=\lambda h$
with the boundary condition $h^{\prime}(0^{+})=\gamma f(0+)$ we find a
solution

$h(y)=e^{y}+\beta e^{-y},$ (57)

where $\beta:=(1-\gamma)/(1+\gamma)$ and then we check that $h$ is indeed
$\lambda$-invariant for $Y$ by verifying that (50) holds.

Condition (44) holds with ${\mathfrak{m}}_{h}=0$ and condition (45) is also
satisfied (here we use Proposition 1 for the left boundary). The Darboux
transform of $Y$ is a process $\widetilde{Y}$, which is a Brownian motion on
$(0,\infty)$ killed at rate $\tilde{c}(y)=(h^{\prime}(y)/h(y))^{2}$. The left
boundary $l=0$ is a non-singular boundary for $\widetilde{Y}$, and we have a
killing bondary condition at this point. In order to simplify the expression
for $\tilde{c}$, we introduce a positive parameter $\alpha$ such that
$|\beta|=\exp(-2\alpha)$ and with the help of this parameter we can express
$h$ in the following form: $h(y)=2e^{-\alpha}\cosh(y+\alpha)$ if
$\gamma\in(0,1)$, $h(y)=\exp(y)$ if $\gamma=1$ and
$h(y)=2e^{-\alpha}\sinh(y+\alpha)$ if $\gamma>1$. Thus, the killing rate of
$\widetilde{Y}$ is given by

$\tilde{c}(y)=\begin{cases}\tanh^{2}(y+\alpha),\;\;&{\textnormal{ if
}}\;\gamma\in(0,1),\\\ 1,\;\;&{\textnormal{ if }}\;\gamma=1,\\\
\coth^{2}(y+\alpha),\;\;&{\textnormal{ if }}\;\gamma>1.\end{cases}$

The transition probability density of $\widetilde{Y}$ is

$\displaystyle p^{\widetilde{Y}}_{t}(x,y)$
$\displaystyle=\frac{e^{-t}}{\sqrt{2\pi
t}}\Big{[}e^{-\frac{1}{2t}(y-x)^{2}}-e^{-\frac{1}{2t}(y+x)^{2}}\Big{]}$
$\displaystyle-\frac{8\gamma\beta}{h(x)h(y)}e^{\gamma(x+y)+(\frac{\gamma^{2}}{2}-1)t}\sinh(x)\sinh(y)\Phi\Big{(}-\frac{y+x+\gamma
t}{\sqrt{t}}\Big{)}$ (58) $\displaystyle+\frac{2\beta
e^{-\frac{t}{2}}}{h(x)h(y)}\Big{[}\Phi\Big{(}\frac{y+x-t}{\sqrt{t}}\Big{)}+\Phi\Big{(}\frac{y-x+t}{\sqrt{t}}\Big{)}-\Phi\Big{(}\frac{y+x+t}{\sqrt{t}}\Big{)}-\Phi\Big{(}\frac{y-x-t}{\sqrt{t}}\Big{)}\Big{]},$

where $h$ is given by (57). Note that when $\gamma=1$, we obtain the following
result: Darboux transform of Brownian motion on $(0,\infty)$ killed
elastically at zero with boundary condition $f^{\prime}(0+)=f(0+)$ is the
Brownian motion on $(0,\infty)$ killed at rate one and killed at the first
time it hits zero. The simplicity of this result suggests that it may have a
pathwise explanation, but we were not able to find it.

For the process $\tilde{Y}$ in this section, we did not derive the spectral
representation of the semigroup, as the resulting expressions were rather
complicated.

The results obtained in this section lead to the following

######  Corollary 2.

Let $W$ be a standard Brownian motion, $\gamma\in(0,\infty)\setminus\\{1\\}$,
$\beta=(1-\gamma)/(1+\gamma)$ and $|\beta|=\exp(-2\alpha)$. Then for $t>0$ and
$x,y>\alpha$

${\mathbb{E}}_{x}\Big{[}e^{-\int_{0}^{t}\kappa(W_{s}){\textnormal{d}}s}{\mathbf{1}}_{\\{W_{t}\leq
y,\min\limits_{0\leq s\leq
t}W_{s}>\alpha\\}}]=\int_{0}^{y-\alpha}p_{t}^{\widetilde{Y}}(x-\alpha,u){\textnormal{d}}u,$
(59)

where

$\kappa(x)=\begin{cases}\tanh^{2}(x),\;\;&{\textnormal{ if
}}\;\gamma\in(0,1),\\\ \coth^{2}(x),\;\;&{\textnormal{ if
}}\;\gamma>1,\end{cases}$

and $p_{t}^{\widetilde{Y}}(x,y)$ is given by (5.3).

###### Proof.

Using the representation of $\widetilde{Y}$ as a Brownian motion on
$(0,\infty)$, killed at rate $\tilde{c}(x)=\kappa(x+\alpha)$ and also killed
at the first time it reaches zero, we can write

$\int_{0}^{y-\alpha}p_{t}^{\widetilde{Y}}(x-\alpha,u){\textnormal{d}}u={\mathbb{P}}_{x-\alpha}(\widetilde{Y}_{t}\leq
y-\alpha)={\mathbb{E}}_{x-\alpha}\Big{[}e^{-\int_{0}^{t}\kappa(B_{s}+\alpha){\textnormal{d}}s}{\mathbf{1}}_{\\{B_{t}\leq
y-\alpha,\min\limits_{0\leq s\leq t}B_{s}>0\\}}],$

where $B$ is the standard Brownian motion process. To obtain (59) we denote
$W_{t}=B_{t}+\alpha$ and use translation invariance of Brownian motion.
$\sqcap\kern-8.0pt\hbox{$\sqcup$}$

Taking derivatives with respect to $y$ and $\alpha$ allows us to obtain the
joint probability density function of the Brownian motion killed at rate
$\kappa(x)$ and of its running minimum. This extends the results of Sections
5.1 and 5.2 by providing information on the running minimum of the process
$\widetilde{Y}$ constructed in those sections.

Remark 2. As discussed in the introduction, the Darboux transformation of
second order linear differential operators is invertible. Starting with an
operator ${\mathcal{L}}$ and a positive formal eigenfunction $h$ such that
${\mathcal{L}}h=\lambda h$, we obtain Darboux transformed operator
$\widetilde{\mathcal{L}}$, for which $1/h$ is also a formal eigenfunction.
Applying Darboux transformation to $\tilde{\mathcal{L}}$ with the seed
function $1/h$, we recover the original operator ${\mathcal{L}}$. However,
this inversion does not work on the level of transition semigroups of
diffusion processes. Consider the process $Y$ from this section and the
function $h$ given by (57). The function $f=1/h$ is positive and satisfies the
equation ${\mathcal{L}}_{\widetilde{Y}}f=\lambda f$ (with $\lambda=1/2$), but
it is not a $\lambda$-invariant function for $\widetilde{Y}$, since it does
not satisfy the boundary condition at $0$ (recall that the process
$\widetilde{Y}$ has a killing boundary condition at zero).

###  5.4 Brownian motion on $(0,1)$, killed at $0$ or $1$

Let $Y$ be a Brownian motion on $(0,1)$ killed at both boundaries, see page
122 in [2]. The transition probability density is

$p^{Y}_{t}(x,y)=2\sum\limits_{n\geq 1}e^{-\frac{1}{2}\pi^{2}n^{2}t}\sin(\pi
nx)\sin(\pi ny).$

The above formula can be seen as a spectral representation of the transition
density of $Y$ in $L_{2}((0,1),{\textnormal{d}}x)$. Using orthogonality of
$\\{\sin(n\pi y)\\}_{n\geq 1}$ on the interval $(0,1)$, it is easy to check
that $h(y)=\sin(\pi y)$ is $\lambda$-invariant function for $Y$ with
$\lambda=-\pi^{2}/2$. We check that condition (44) holds with
${\mathfrak{m}}_{h}=0$ and condition (45) is also satisfied, since both
boundaries are natural for the process $X$ with the generator given in (36).
The Darboux transformed process $\widetilde{Y}$ is a Brownian motion on
$(0,1)$, killed at rate $\tilde{c}(x)=\pi^{2}\cot^{2}(\pi x)$. Both boundaries
are natural for $\widetilde{Y}$. The transition probability density is
obtained using (46) and has the form

$p^{\widetilde{Y}}_{t}(x,y)=2\sum\limits_{n\geq
2}\frac{e^{-\frac{1}{2}\pi^{2}(n^{2}-2)t}}{n^{2}-1}f_{n}(x)f_{n}(y),$ (60)

where we denoted

$f_{n}(x):=\frac{1}{\pi}{\mathcal{D}}_{h}\sin(n\pi x)=n\cos(\pi nx)-\sin(\pi
nx)\cot(\pi x).$

The functions $\\{f_{n}\\}_{n\geq 0}$ are orthogonal in
$L_{2}((0,1),{\textnormal{d}}x)$:

$\int_{0}^{1}f_{n}(x)f_{m}(x){\textnormal{d}}x=\frac{1}{2}(n^{2}-1)\delta_{n,m},\;\;\;n,m\geq
2.$

This result follows by applying (47) with $f=\sin(n\pi x)$ and $g=\sin(m\pi
x)$. Thus formula (60) gives a spectral representation of the transition
probability density of $\widetilde{Y}$ in $L_{2}((0,1),{\textnormal{d}}x)$. In
this case, the Darboux transformation removes the first eigenvalue from the
spectrum of ${\mathcal{L}}_{Y}$ and shifts all eigenvalues by $-2\lambda$.

Now we can repeat this construction. We note that $f_{2}(x)=-2\sin^{2}(\pi
x)$, thus $h_{1}(y)=-f_{2}(y)=2\sin^{2}(\pi y)$ is a positive
$\lambda$-invariant function for $\widetilde{Y}$, with $\lambda=-\pi^{2}$.
Applying Darboux transformation to $\widetilde{Y}$ with this
$\lambda$-invariant function $\tilde{h}$, we obtain a new diffusion process,
which we denote by $Y^{(2)}$. This process is a Brownian motion on $(0,1)$
killed at rate $c^{(2)}(x)=3\pi^{2}\cot^{2}(\pi x)$ and having transition
probability density

$p^{(2)}_{t}(x,y)=2\sum\limits_{n\geq
3}\frac{e^{-\frac{1}{2}\pi^{2}(n^{2}-6)t}}{(n^{2}-1)(n^{2}-4)}f^{(2)}_{n}(x)f^{(2)}_{n}(y),$
(61)

where

$f^{(2)}_{n}(x):=\frac{1}{\pi}{\mathcal{D}}_{h_{1}}f_{n}(x)=-(2+n^{2})\sin(\pi
nx)+3\sin(\pi nx)\csc^{2}(\pi x)-3n\cos(\pi nx)\cot(\pi x).$

The formula (61) again gives us a spectral representation of the transition
probability density of the process $Y^{(2)}$ in
$L_{2}((0,1),{\textnormal{d}}x)$.

This process can be repeated indefinitely. The first eigenfunction of the
transition semigroup must have a constant sign on $(0,1)$, thus we can take
$\lambda$-invariant function $h$ to equal to this first eigenfunction or its
negative. After performing Darboux transform $m$ times, we will obtain a
process $Y^{(m)}$, which is a Brownian motion on $(0,1)$, killed at rate
$c^{(2)}(x)=\frac{1}{2}m(m+1)\pi^{2}\cot^{2}(\pi x)$. The transition
probability can be written in the following form

$p^{(m)}_{t}(x,y)=2\sum\limits_{n\geq
m+1}\frac{e^{-\frac{1}{2}\pi^{2}(n^{2}-m(m+1))t}}{\prod\limits_{j=1}^{m}(n^{2}-j^{2})}f^{(m)}_{n}(x)f^{(m)}_{n}(y),$
(62)

where the functions $f^{(m)}_{n}$ can be given explicitly in terms of
Gegenbauer polynomials and also in terms of certain Wronskian determinants.
The proof of (62) will appear in the forthcoming paper [26].

###  5.5 Brownian motion on ${\mathbb{R}}$ killed at rate $y^{2}/2$

Let $Y$ be a Brownian motion on ${\mathbb{R}}$ killed at rate $c(y)=y^{2}/2$.
It is known that its transition probability density function is given by
Mehler’s kernel

$\displaystyle p^{Y}_{t}(x,y)$
$\displaystyle=\frac{1}{\sqrt{2\pi\sinh(t)}}\exp\Big{(}-\frac{1}{2}\coth(t)(x^{2}+y^{2})+\frac{xy}{\sinh(t)}\Big{)}$
$\displaystyle=\frac{1}{\sqrt{\pi}}e^{-\frac{1}{2}(x^{2}+y^{2})}\sum\limits_{n\geq
0}\frac{e^{-(n+\frac{1}{2})t}}{2^{n}n!}H_{n}(x)H_{n}(y).$

Here $H_{n}(x)$ denote Hermite polynomials. The process $Y$ is closely related
to the Ornstein-Uhlenbeck process $Z$, which has infinitesimal generator
${\mathcal{L}}_{Z}=\frac{1}{2}\partial_{z}^{2}-z\partial_{z}$. It is known
(and easy to check) that $Y$ is Doob’s $h$-transform of $Z$:

${\mathbb{P}}_{z}(Y_{t}\in
A)=e^{-\frac{t}{2}-\frac{1}{2}z^{2}}{\mathbb{E}}_{z}\Big{[}e^{\frac{1}{2}Z_{t}^{2}}{\mathbf{1}}_{\\{Z_{t}\in
A\\}}\Big{]}.$

Next, we will discuss how to find a $\lambda$-invariant function for the
process $Y$. It is known (see Example 2.3 in [12]) that the functions
$f_{n}(y)=e^{-y^{2}/2}H_{n}(y)$ and
$g_{n}(y)={\textnormal{i}}^{n}e^{y^{2}/2}H_{n}({\textnormal{i}}y)$ are formal
eigenfunctions of the operator
${\mathcal{L}}=\frac{1}{2}(\partial_{y}^{2}-y^{2})$. Hermite polynomials
$H_{n}(y)$ have real zeros if $n\geq 1$ and the functions
$H_{n}({\textnormal{i}}y)$ have a zero at $y=0$ when $n$ is odd, so these are
not good candidates for a positive $\lambda$-invariant function and we should
take $h$ as a multiple of $f_{0}$ or $g_{2n}$. Taking $h(y)=\exp(\pm y^{2}/2)$
does not lead to a new process, as $\tilde{c}(y)$ given by (43) is equal to
$c(y)$ up to a constant. Thus, the first non-trivial example can be obtained
when we take $h$ as a multiple of $g_{2}(y)$. We set
$h(y)=e^{y^{2}/2}(2y^{2}+1)$ and check (by computing the integral in (50)
symbolically or by applying [8][Theorem 2.7] or [19][Corollary 2.2]) that $h$
is $\lambda$-invariant for $Y$ with $\lambda=5/2$. Now we can construct the
Darboux transformed process $\widetilde{Y}$, which is a Brownian motion on
${\mathbb{R}}$ killed at rate

$\tilde{c}(y)=\frac{y^{2}}{2}+\frac{8y^{2}(2y^{2}+3)}{(2y^{2}+1)^{2}}.$

The transition probability density of $\widetilde{Y}$ (computed via (46)) has
a surprisingly simple form:

$p^{\widetilde{Y}}_{t}(x,y)=p^{Y}_{t}(x,y)\times
e^{-4t}\Big{[}1+\frac{4\sinh(t)(e^{t}-2xy)}{(2x^{2}+1)(2y^{2}+1)}\Big{]}.$
(63)

Next, we present the spectral expansion of the transition probability density.
The operator ${\mathcal{L}}_{Y}=\frac{1}{2}(\partial_{y}^{2}-y^{2})$ has a
complete set of orthogonal eigenfunctions $f_{n}(y)=e^{-y^{2}/2}H_{n}(y)$ in
$L_{2}({\mathbb{R}},{\textnormal{d}}x)$. We know that $1/h(y)$ is a formal
eigenfunction of $\widetilde{Y}$. We can expect that
$\tilde{f}_{n}(y)={\mathcal{D}}_{h}f_{n}(y)$ to be the eigenfunctions of
${\mathcal{L}}_{\widetilde{Y}}$ in $L_{2}({\mathbb{R}},{\textnormal{d}}x)$. We
compute

$\tilde{f}_{n}(y)={\mathcal{D}}_{h}f_{n}(y)=f_{n}^{\prime}(y)-\frac{h^{\prime}(y)}{h(y)}f_{n}(y)=\frac{Q_{n}(y)}{h(y)}$

where $Q_{n}(y)$ is the following polynomial

$Q_{n}(y)={\textnormal{Wr}}[h,f_{n}]=2n(2y^{2}+1)H_{n-1}(y)-(4y^{3}+6y)H_{n}(y),\;\;\;n\geq
0.$

Note that $Q_{n}$ are polynomials of degree $n+3$.

######  Proposition 3.

The set of polynomials $\\{1\\}\cup\\{Q_{n}\\}_{n\geq 0}$ is a complete
orthogonal set in $L_{2}({\mathbb{R}},\nu({\textnormal{d}}y))$, where

$\nu({\textnormal{d}}y):=h^{-2}(y){\textnormal{d}}y=\frac{e^{-y^{2}}}{(2y^{2}+1)^{2}}{\textnormal{d}}y.$

The transition probability density of $\widetilde{Y}$ has spectral
representation

$p^{\widetilde{Y}}_{t}(x,y)=\frac{1}{\sqrt{\pi}}\frac{e^{-\frac{1}{2}(x^{2}+y^{2})}}{(2x^{2}+1)(2y^{2}+1)}\times\Big{[}2e^{-\frac{5}{2}t}+\sum\limits_{n\geq
0}e^{-(n+11/2)t}\frac{Q_{n}(x)Q_{n}(y)}{2^{n+1}n!(n+3)}\Big{]}.$ (64)

###### Proof.

Let $S=\\{0,3,4,5,6,\dots\\}$ be the set of non-negative integers that is
missing $\\{1,2\\}$ and let us define polynomials $H^{(1)}_{n}(y)$ via

$H^{(1)}_{n}(y)={\textnormal{Wr}}[H_{1},H_{2},H_{n}](y)=16H_{n}(y)-16xH_{n}^{\prime}(y)+4(2x^{2}+1)H_{n}^{\prime\prime}(y).$

Note that $H^{(1)}_{n}\neq 0$ if and only if $n\in S$. The polynomials
$\\{H^{(1)}_{n}\\}_{n\in S}$ are the exceptional Hermite polynomials
corresponding to the partition $2=1+1$, see [10][Definition 5.1]. It is known
(see Propositions 5.7 and 5.8 in [10]) that $\\{H^{(1)}_{n}\\}_{n\in S}$ form
a complete orthogonal set in $L_{2}(R,\nu({\textnormal{d}}x))$ with the
squared norm

$\int_{{\mathbb{R}}}H^{(1)}_{n}(y)^{2}\nu({\textnormal{d}}y)=\sqrt{\pi}2^{n+6}n!(n-1)(n-2).$
(65)

Using [11][Theorem 3.1], we can write

$H^{(1)}_{n}(y)=-16(n+1)(n+2)Q_{n}(y),\;\;\;n\geq 0.$ (66)

and clearly $H^{(1)}_{n}(y)=16$. Thus the polynomials $\\{H^{(1)}_{n}\\}_{n\in
S}$ and $\\{1\\}\cup\\{Q_{n}\\}_{n\geq 0}$ are identical, up to a
multiplicative constant. This implies that the functions
$\\{1/h\\}\cup\\{Q_{n}/h\\}_{n\geq 0}$ form a complete orthogonal set of
eigenfunctions of ${\mathcal{L}}_{\widetilde{Y}}$ in
${\mathcal{L}}_{2}({\mathbb{R}},{\textnormal{d}}x)$ with eigenvalues $-5/2$
and $-(n+11/2)$, $n\geq 0$. The $L_{2}$-norm of these functions can be found
via (65) and (66) (though one could find the norms of $Q_{n}/h$ more directly
via Lemma 1:

$\int_{{\mathbb{R}}}\frac{1}{h(y)^{2}}{\textnormal{d}}y=\nu({\mathbb{R}})=\frac{\sqrt{\pi}}{2},\;\;\;\int_{{\mathbb{R}}}\frac{Q_{n}(y)^{2}}{h(y)^{2}}{\textnormal{d}}y=\int_{{\mathbb{R}}}Q_{n}(y)^{2}\nu({\textnormal{d}}y)=\sqrt{\pi}2^{n+1}n!(n+3).$

One can now write down the transition probability density of $\widetilde{Y}$
in the form (64).

$\sqcap\kern-8.0pt\hbox{$\sqcup$}$

In conclusion, we would like to discuss the connection of our results with the
existing literature on Darboux transformation and propagators of for the one-
dimensional Schrödinger equation. Our formula (46) has an analogue in
[22][Theorem 2]. Example 1 in Section 5.1 should be compared with the results
of [21] and example 4 in Section 5.4 with example 5.1 in [22]. Our example 5
is closely related to the propagators of rational extensions of harmonic
oscillator considered in [20]. Moreover, the results in [20] imply that if the
potential

$u(y)=y^{2}+\frac{a(y)}{b(y)}$

is a rational extension of harmonic oscillator (as defined in [10]) then the
process $\widetilde{Y}$, which is the Brownian motion on ${\mathbb{R}}$ killed
at rate $\tilde{c}(y)=u(y)/2$, will have transition probability density of the
form similar to (63):

$p_{t}^{\widetilde{Y}}(x,y)=p_{t}^{Y}(x,y)\times
e^{-ct}\frac{Q(x,y,e^{t})}{R(x,y)},$

where $c$ is a constant and $Q(x,y,z)$ and $R(x,y)$ are certain polynomials
that can be computed explicitly. Thus, there exists a hierarchy of diffusion
processes on ${\mathbb{R}}$, which could be called rational extensions of
Ornstein-Uhlenbeck process, for which the transition probability density can
be computed in a fairly simple form. We plan to investigate this family of
diffusion processes in future work.

## Acknowledgements

The research was supported by the Natural Sciences and Engineering Research
Council of Canada. The authors would like to thank Mateusz Kwaśnicki for
stimulating discussions and for helping with the proof of Theorem 1.

## References

  * [1] V. G. Bagrov and B. F. Samsonov.  Darboux transformation and elementary exact solutions of the Schrödinger equation.  Pramana - J. Phys., 49(6):563 – 580, 1997.  https://doi.org/10.1007/BF02848330.
  * [2] A. N. Borodin and P. Salminen.  Handbook of Brownian motion.  Srpinger, 2nd edition, 2002.
  * [3] P. Clifford and A. Sudbury.  A sample path proof of the duality for stochastically monotone Markov processes.  The Annals of Probability, 13(2):558–565, 1985.  https://doi.org/10.1214/aop/1176993008.
  * [4] A. Comtet and Y. Tourigny.  Excursions of diffusion processes and continued fractions.  Annales de l’Institut Henri Poincaré, Probabilités et Statistiques, 47(3):850 – 874, 2011.  https://doi.org/10.1214/10-AIHP390.
  * [5] J. T. Cox and U. Rösler.  A duality relation for entrance and exit laws for Markov processes.  Stochastic Processes and their Applications, 16(2):141–156, 1984\.  https://doi.org/10.1016/0304-4149(84)90015-2.
  * [6] M. M. Crum.  Associated Sturm-Liouville systems.  The Quarterly Journal of Mathematics, 6(1):121–127, 01 1955.  https://doi.org/10.1093/qmath/6.1.121.
  * [7] G. Darboux.  Sur une proposition relative aux équations linéaires.  Comptes Rendus Acad. Sci., 94:1456–1459, 1882.
  * [8] S. Desmettre, G. Leobacher, and L. C. G. Rogers.  Change of drift in one-dimensional diffusions.  Finance and Stochastics, 25:359–381, 2021.  https://doi.org/10.1007/s00780-021-00451-w.
  * [9] M. Garcia-Ferrero, D. Gómez-Ullate, and R. Milson.  A Bochner type characterization theorem for exceptional orthogonal polynomials.  Journal of Mathematical Analysis and Applications, 472(1):584–626, 2019.  https://doi.org/10.1016/j.jmaa.2018.11.042.
  * [10] D. Gómez-Ullate, Y. Grandati, and R. Milson.  Rational extensions of the quantum harmonic oscillator and exceptional polynomials.  J. Phys. A: Math Theor., 47(1):015203, 2014.  http://dx.doi.org/10.1088/1751-8113/47/1/015203.
  * [11] D. Gómez-Ullate, Y. Grandati, and R. Milson.  Durfee rectangles and pseudo-Wronskian equivalences for Hermite polynomials.  Studies in Applied Mathematics, 141(4):596–625, 2018.  https://doi.org/10.1111/sapm.12225.
  * [12] D. Gómez-Ullate and R. Milson.  Exceptional orthogonal polynomials and rational solutions to Painlevé equations.  In Orthogonal Polynomials, pages 335–386. Springer International Publishing, 2020, editors M. Foupouagnigni and W. Koepf.  https://doi.org/10.1007/978-3-030-36744-2_15.
  * [13] C. Gu, H. Hu, and Z. Zhou.  Darboux Transformations in Integrable Systems, volume 26 of Mathematical Physics Studies.  Springer, 2005.  https://doi.org/10.1007/1-4020-3088-6.
  * [14] K. Itô and H. P. McKean, Jr.  Diffusion processes and their sample paths.  Springer, 1974.
  * [15] V. N. Kolokol’tsov.  Stochastic monotonicity and duality for one-dimensional Markov processes.  Mathematical Notes, 89(5):652–660, 2011.  https://doi.org/10.1134/S0001434611050063.
  * [16] S. Kotani and S. Watanabe.  Krein’s spectral theory of strings and generalized diffusion processes.  In M. Fukushima, editor, Functional Analysis in Markov Processes, pages 235–259. Springer Berlin Heidelberg, 1982.  https://doi.org/10.1007/BFb0093046.
  * [17] P. Mandl.  Analytical treatment of one-dimensional Markov processes.  Springer, 1968.
  * [18] V. B. Matveev and M. A. Salle.  Darboux Transformations and Solitons.  Nonlinear Dynamics. Springer, 1991.
  * [19] A. Mijatović and M. Urusov.  On the martingale property of certain local martingales.  Probability Theory and Related Fields, 152:1–30, 2012.  https://doi.org/10.1007/s00440-010-0314-7.
  * [20] A. M. Pupasov-Maksimov.  Propagators of isochronous an-harmonic oscillators and Mehler formula for the exceptional Hermite polynomials.  Annals of Physics, 363:122–135, 2015.  https://doi.org/10.1016/j.aop.2015.09.021.
  * [21] A. M. Pupasov-Maksimov and B. F. Samsonov.  Exact Propagators for Soliton Potentials.  SIGMA, 1:20, 2005.  https://doi.org/10.3842/SIGMA.2005.020.
  * [22] A. M. Pupasov-Maksimov, B. F. Samsonov, and U. Günther.  Exact propagators for SUSY partners.  Journal of Physics A: Mathematical and Theoretical, 40(34):10557, 2007.  https://dx.doi.org/10.1088/1751-8113/40/34/013.
  * [23] D. Siegmund.  The equivalence of absorbing and reflecting barrier problems for stochastically monotone Markov processes.  The Annals of Probability, 4(6):914 – 924, 1976.  https://doi.org/10.1214/aop/1176995936.
  * [24] R. Sousa, M. Guerra, and S. Yakubovich.  Convolution-like Structures, Differential Operators and Diffusion Processes, volume 2315 of Lecture Notes in Mathematics.  Springer, 2022.
  * [25] B. Toth.  Generalized Ray-Knight theory and limit theorems for self-interacting random walks on ${\mathbb{Z}}^{1}$.  The Annals of Probability, 24(3):1324–1367, 1996.  http://www.jstor.org/stable/2244976.
  * [26] M. Yuan.  On Gegenbauer polynomials and Wronskian determinants of trigonometric functions.  manuscript in preparation, 2024.

##  Appendix A Proof of identities (54) and (56)

We begin with the following result

$\frac{e^{w}}{2\pi}\int_{{\mathbb{R}}}\frac{e^{-t\frac{z^{2}}{2}+{\textnormal{i}}zw}}{{\textnormal{i}}z+1}{\textnormal{d}}z=e^{\frac{t}{2}}\Phi\Big{(}\frac{w-t}{\sqrt{t}}\Big{)},\;\;\;t>0,\;\;w\in{\mathbb{R}},$
(A.1)

which can be established by taking the derivative in $w$ of both sides and
computing the integral in the left-hand side. From (A.1) we obtain the
following pair of Fourier transform identities:

$\displaystyle\frac{1}{2\pi}\int_{{\mathbb{R}}}\frac{e^{-t\frac{z^{2}}{2}+{\textnormal{i}}zw}}{1+z^{2}}{\textnormal{d}}z=\frac{1}{2}\Big{[}e^{\frac{t}{2}+w}\Phi\Big{(}-\frac{w+t}{\sqrt{t}}\Big{)}+e^{\frac{t}{2}-w}\Phi\Big{(}\frac{w-t}{\sqrt{t}}\Big{)}\Big{]}$
(A.2)
$\displaystyle\frac{{\textnormal{i}}}{2\pi}\int_{{\mathbb{R}}}\frac{e^{-t\frac{z^{2}}{2}+{\textnormal{i}}zw}}{1+z^{2}}z{\textnormal{d}}z=\frac{1}{2}\Big{[}e^{\frac{t}{2}+w}\Phi\Big{(}-\frac{w+t}{\sqrt{t}}\Big{)}-e^{\frac{t}{2}-w}\Phi\Big{(}\frac{w-t}{\sqrt{t}}\Big{)}\Big{]}.$
(A.3)

We write

$F({\textnormal{i}}z,x)F(-{\textnormal{i}}z,y)=e^{{\textnormal{i}}z(x-y)}\Big{[}(1+z^{2})+{\textnormal{i}}z(\tanh(x)-\tanh(y))+(\tanh(x)\tanh(y)-1)]$

and now we can compute the integral in (54):

$\displaystyle\frac{1}{2\pi}\int_{{\mathbb{R}}}e^{-(1+\frac{z^{2}}{2})t}F({\textnormal{i}}z,x)F(-{\textnormal{i}}z,y)\frac{{\textnormal{d}}z}{1+z^{2}}=\frac{1}{\sqrt{2\pi
t}}e^{-t-\frac{1}{2t}(x-y)^{2}}$
$\displaystyle+\frac{1}{2}(\tanh(x)-\tanh(y))\Big{[}e^{-\frac{t}{2}+(x-y)}\Phi\Big{(}\frac{y-x-t}{\sqrt{t}}\Big{)}-e^{-\frac{t}{2}-(x-y)}\Phi\Big{(}\frac{x-y-t}{\sqrt{t}}\Big{)}\Big{]}$
$\displaystyle+\frac{1}{2}\big{(}\tanh(x)\tanh(y)-1)\Big{[}e^{-\frac{t}{2}+(x-y)}\Phi\Big{(}\frac{y-x-t}{\sqrt{t}}\Big{)}+e^{-\frac{t}{2}-(x-y)}\Phi\Big{(}\frac{x-y-t}{\sqrt{t}}\Big{)}\Big{]}.$

Plugging the above expression into the right-hand side of (54) and simplifying
the result will give us $p^{\widetilde{Y}}_{t}(x,y)$.

The proof of (56) is very similar. We write the integral in the right-hand
side of (56) as a sum of five integrals

$\displaystyle\frac{2}{\pi}\int_{0}^{\infty}e^{-(1+\frac{z^{2}}{2})t}f(z,x)f(z,y)\frac{{\textnormal{d}}z}{1+z^{2}}=\frac{2}{\pi}\int_{0}^{\infty}e^{-(1+\frac{z^{2}}{2})t}\cos(zx)\cos(zy){\textnormal{d}}z$
$\displaystyle-\frac{2}{\pi}\coth(x)\int_{0}^{\infty}e^{-(1+\frac{z^{2}}{2})t}\cos(zy)\sin(zx)\frac{z}{1+z^{2}}{\textnormal{d}}z-\frac{2}{\pi}\coth(y)\int_{0}^{\infty}e^{-(1+\frac{z^{2}}{2})t}\cos(zx)\sin(zy)\frac{z}{1+z^{2}}{\textnormal{d}}z$
$\displaystyle+\frac{2}{\pi}\coth(x)\coth(y)\int_{0}^{\infty}e^{-(1+\frac{z^{2}}{2})t}\frac{\sin(zx)\sin(zy)}{1+z^{2}}{\textnormal{d}}z-\frac{2}{\pi}\int_{0}^{\infty}e^{-(1+\frac{z^{2}}{2})t}\frac{\cos(zx)\cos(zy)}{1+z^{2}}{\textnormal{d}}z.$

Each of these integrals can be computed by using product-to-sum identities for
the trigonometric functions and applying Fourier transform identities (A.2)
and (A.3). The resulting expression, when plugged into the right-hand side of
(56), after some long and tedious simplifications, will give us
$p^{\widetilde{Y}}_{t}(x,y)$.

##  Appendix B Proof of (32)

We present the proof in the case when $l$ is a non-singular boundary for the
process $Y$ from Section 4. The proof for the case when $r$ is non-singular is
identical. We recall that ${\mathcal{L}}_{Y}=\frac{1}{2}\partial_{y}^{2}-c(y)$
(where $c$ is a positive continuous function on $(l,r)$) and that the boundary
$l$ is non-singular for $Y$ if and only if $l$ is finite and

$\int_{l}^{z}c(y){\textnormal{d}}y<\infty,$ (B.1)

for some $z\in(l,r)$. Let $\lambda$ be a real number and $h$ be a
$\lambda$-invariant function for $Y$.

First, let us establish the following result: any positive solution to
${\mathcal{L}}_{Y}f=\lambda f$ must be bounded on $(l,z)$. This result is
undoubtedly known in Sturm-Liouville theory, but we were unable to find a
reference, so we will include a short proof here. We denote
$g(y)=f^{\prime}(y)/f(y)$ and check that $g$ satisfies the Ricatti equation

$g^{\prime}(y)=2(c(y)+\lambda)-g(y)^{2},\;\;\;y\in(l,r),$

which implies (by integrating over $y\in(w,z)$)

$g(z)-g(w)\leq
2\lambda(z-w)+2\int_{w}^{z}c(y){\textnormal{d}}y,\;\;\;l<w<z<r.$

Therefore $-g(w)$ is bounded from above on the interval $(l,z)$ (due to
(B.1)), and since

$f(y)=f(z)\exp\Big{(}-\int_{y}^{z}g(w){\textnormal{d}}w\Big{)},\;\;\;l<y<z,$

we conclude that $f$ is also bounded on $(l,z)$.

Now we are ready to prove that $h$ satisfies the appropriate boundary
condition at $l$. We denote by $\psi_{\mu}$ and $\varphi_{\mu}$ the
fundamental increasing/decreasing solutions to ${\mathcal{L}}_{Y}f=\mu f$. We
will show $h$ satisfies the same boundary condition at $l$ as $\psi_{\mu}$.
The starting point is the identity (33), which, coupled with the
representation (14) for the Green function of the process $Y$, gives us

$h(x)=(\mu-\lambda)\Big{[}\varphi_{\mu}(x)\int_{l}^{x}\psi_{\mu}(y)h(y){\textnormal{d}}y+\psi_{\mu}(x)\int_{x}^{r}\varphi_{\mu}(y)h(y){\textnormal{d}}y\Big{]}.$

The above holds for $\mu>\max(0,\lambda)$ and $x\in(l,r)$. As we established
above, both functions $\varphi_{\mu}$ and $h$ are bounded (thus, integrable)
on $(l,z)$, thus we can rewrite the above identity in the form

$\displaystyle h(x)$
$\displaystyle=(\mu-\lambda)\psi_{\mu}(x)\int_{l}^{r}\varphi_{\mu}(y)h(y){\textnormal{d}}y$
$\displaystyle+(\mu-\lambda)\varphi_{\mu}(x)\int_{l}^{x}\big{(}\psi_{\mu}(y)-\psi_{\mu}(x))h(y){\textnormal{d}}x-(\mu-\lambda)\psi_{\mu}(x)\int_{l}^{x}\big{(}\varphi_{\mu}(y)-\varphi_{\mu}(x))h(y){\textnormal{d}}x.$

Since $\psi_{\mu}$, $\varphi_{\mu}$ and $h$ are bounded on $(l,z)$, the above
identity implies that as $x\to l+$ we have

$\displaystyle h(x)$ $\displaystyle=C\psi_{\mu}(x)+O(x-l),$ (B.2)
$\displaystyle h^{\prime}(x)$
$\displaystyle=C\psi^{\prime}_{\mu}(x)+\varphi^{\prime}_{\mu}(x)O(x-l)+\psi^{\prime}_{\mu}(x)O(x-l),$
(B.3)

where we denoted

$C=C(\mu,\lambda):=(\mu-\lambda)\int_{l}^{r}\varphi_{\mu}(y)h(y){\textnormal{d}}y.$

If we have a killing boundary condition for $Y$ at point $l$, that
is$f(l+)=0$, then (B.2) implies $h(l+)=\psi_{\mu}(l+)=0$. If we have a
reflecting boundary condition at $l$, then $\psi_{\mu}^{\prime}(l+)=0$ and
(B.3) implies $h^{\prime}(l+)=0$. Here we also need the fact that
$\psi_{\mu}^{\prime}(x)$ and $\phi_{\mu}^{\prime}(x)$ are bounded near point
$l$ (the first one is bounded near $l$ since $\psi_{\mu}(x)$ is convex and
increasing and the second one is bounded because the Wronskian
${\textnormal{Wr}}[\varphi_{\lambda};\psi_{\lambda}]$ is a non-zero constant).
In the case when we have an elastic boundary condition for $Y$ at point $l$,
that is $f(l+)=\gamma f^{\prime}(l+)$, we can obtain $h(l+)=\gamma
h^{\prime}(l+)$ using the same argument.