$\displaystyle\times\left[\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\mathtt{Im}\left(\frac{\left(\phi^{\dagger}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right)^{k_{3}+1}}{\rho^{k_{3}+1-(k_{2}+k_{3}+l)}e^{i(k_{3}-(l+k_{2}+k_{3}))\left(\pi-\frac{\theta}{2}\right)}}e^{t\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}}\right)\,d\rho+\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon,\theta}}\frac{(\phi^{\dagger}(z))^{k_{3}+1}}{z^{k_{3}+1-(k_{2}+k_{3}+l)}}e^{tz}dz\right]$
$\displaystyle=\sum_{j=0}^{+\infty}\sum_{k_{1}+k_{2}+k_{3}=k+j}\frac{(k+j)!}{k_{1}!k_{2}!k_{3}!}(-1)^{k+j}\frac{x^{j}}{j!}q^{k_{1}}\mathfrak{b}^{k_{2}}\frac{d^{k_{2}+k_{3}+l}}{dt^{k_{2}+k_{3}+l}}\mu^{\ast(k_{3}+1)}(t),$

where we used (5.23) in the last equality. Now we only have to prove that we
can exchange the series with the integrals. This is clear for the integral
over $\gamma_{\varepsilon,\theta}$, thus let us only consider the one over
$(\varepsilon,+\infty)$. Indeed, we have

$\displaystyle\int_{\varepsilon}^{+\infty}\sum_{j=0}^{+\infty}\left|(-1)^{k+j}\frac{\phi^{\dagger}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\left(\phi\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right)^{{k+j}}}{j!\;\rho^{1-l}e^{-il\left(\pi-\frac{\theta}{2}\right)}}\;x^{j}\;e^{t\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}}\right|d\rho$ $\displaystyle=\,$
$\displaystyle\int_{\varepsilon}^{\infty}\frac{\left|\phi^{\dagger}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right|\left|\phi\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right|^{k}}{\rho^{k+1}}\rho^{l+k}e^{x\left|\phi\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right|-t\rho\cos\frac{\theta}{2}}d\rho$
$\displaystyle\leq\,$ $\displaystyle
e^{xq}\int_{\varepsilon}^{+\infty}\frac{\left|\phi^{\dagger}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right|\left|\phi\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right|^{k}}{\rho^{k+1}}\rho^{l+k}e^{x\mathfrak{b}\rho\cos\left(\frac{\theta}{2}\right)-t\rho\cos\frac{\theta}{2}+\phi^{\dagger}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)}d\rho.$ (5.28)

It is clear that we have to check integrability in the right-hand side of
(5.28) only in a neighbourhood of infinity. To do this, set
$\delta=t-\mathfrak{b}x$ and $p,p^{\prime}>1$ as in the proof of Proposition
5.1. By ($\mathbb{B}_{2}$) we know that there exists $K$ big enough such that
$\frac{\left|\phi^{\dagger}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right|\left|\phi\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right|^{k}}{\rho^{k+1}}$ is
bounded and $\frac{\left|\phi^{\dagger}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right|}{\rho\cos(\theta/2)}<\frac{\delta}{4}$,
whenever $\rho>K$. Hence we get

$\displaystyle\begin{split}&\int_{K}^{+\infty}\frac{\left|\phi^{\dagger}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right|\left|\phi\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right|^{k}}{\rho^{k+1}}e^{-\frac{t}{p^{\prime}}\rho\cos\left(\frac{\theta}{2}\right)}\exp\left(-\rho\cos\left(\frac{\theta}{2}\right)\left(\frac{t}{p}-x\mathfrak{b}-\frac{\left|\phi^{\dagger}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right|}{\rho\cos\left(\frac{\theta}{2}\right)}\right)\right)d\rho\\\
\leq\,&C\int_{K}^{+\infty}\rho^{k+l}e^{-\frac{t}{p^{\prime}}\rho\cos\left(\frac{\theta}{2}\right)}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d\rho}<+\infty.\end{split}$

This concludes the proof. ∎

###  5.3. Proof of Theorem 3.18

The behaviour at zero provided in Theorem 3.18 can be shown as a direct
consequence of the series representation given in Theorem 3.15.

###### Proof of Theorem 3.18.

Let $[t_{1},t_{2}]\subset(0,+\infty)$ and observe that for $t\in[t_{1},t_{2}]$
we have by (5.23)

$\begin{split}|\mathcal{I}_{j,k,l}(t)|&\leq\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\frac{\left|\phi^{\dagger}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right|\left|\phi\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right|^{j+k}}{\rho^{1-l}}e^{-t_{1}\rho\cos\left(\frac{\theta}{2}\right)}d\rho\\\
&+\frac{\varepsilon^{l}}{2\pi}\int_{\frac{\theta}{2}-\pi}^{\pi-\frac{\theta}{2}}|\phi^{\dagger}(\varepsilon
e^{i\varphi})||\phi(\varepsilon e^{i\varphi})|^{j+k}e^{\varepsilon
t_{2}|\cos(\varphi)|}d\varphi.\end{split}$ (5.29)

Hence, for $0<x<\frac{t_{1}}{\mathfrak{b}}$ and $t\in[t_{1},t_{2}]$, recalling
the definition of $\mathcal{P}_{n,k,l}(x,t)$ given in (3.20), we have

$\displaystyle\left|\frac{\partial^{k}\partial^{l}}{\partial x^{k}\partial
t^{l}}f_{\phi}(x,t)-\mathcal{P}_{n,k,l}(x,t)\right|\leq$
$\displaystyle\sum_{j=n+1}^{+\infty}\frac{x^{j}}{j!}|\mathcal{I}_{j,k,l}(t)|$
$\displaystyle\leq\,$
$\displaystyle\frac{x^{n+1}}{(n+1)!}\sum_{j=0}^{+\infty}\frac{x^{j}}{j!}\left[\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\frac{\left|\phi^{\dagger}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right|\left|\phi\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right|^{j+k+n+1}}{\rho^{1-l}}e^{-t_{1}\rho\cos\left(\frac{\theta}{2}\right)}d\rho\right.$
$\displaystyle+\left.\frac{\varepsilon^{l}}{2\pi}\int_{\frac{\theta}{2}-\pi}^{\pi-\frac{\theta}{2}}|\phi^{\dagger}(\varepsilon
e^{i\varphi})||\phi(\varepsilon e^{i\varphi})|^{j+k+n+1}e^{\varepsilon
t_{2}|\cos(\varphi)|}d\varphi\right]$ (5.30)

where, in the last step, we used (5.29). The convergence of the series in
(5.30) can be ascertained as in (5.28). Taking the supremum in $[t_{1},t_{2}]$
in (5.30), we get (3.19). ∎

###  5.4. Proof of Proposition 3.19

A slight modification of the arguments used in the proofs of Proposition 5.3
and Theorem 3.14 leads to the desired result.

###### Proof of Proposition 3.19.

Since $\phi$ is a complete Bernstein function that can be extended by
continuity over $\overline{\mathbb{C}(0,\pi)}$ with extension $\phi_{+}$,
then, by using the relation $\overline{\phi(z)}=\phi(\overline{z})$, it is
clear that it can be also extended by continuity over
$\overline{\mathbb{C}(-\pi,0)}$. If we denote such an extension $\phi_{-}$,
for any $z\in\overline{\mathbb{C}(0,\pi)}$ we get

$\overline{\phi_{+}(z)}=\phi_{-}(\overline{z}).$

The proof is then carried on exactly as in Propositions 5.1 and 5.3, by
setting $\theta=0$ (see Figure 2), where we use $\phi^{\dagger}$ when we
integrate over $\ell_{1}$, $\Gamma_{R}^{+}$, $\Gamma_{R}^{-}$ and
$-\gamma_{\varepsilon}$, $\phi^{\dagger}_{+}$ over $\ell_{2}$ and
$\phi^{\dagger}_{-}$ over $\ell_{1}$. ∎

###  5.5. Proof of Theorem 3.20

In order to prove Theorem 3.20 we first need to provide an integral
representation for $G_{\phi}$, $g_{\phi}$ and its derivatives, analogously to
what we did for Theorems 3.14 and 3.15. This is done in the following
proposition, whose proof is almost identical to the one of Proposition 5.1 and
thus is omitted.

######  Proposition 5.4.

Let $\phi$ be the Laplace exponent of a potentially killed subordinator
satisfying assumptions $\mathbb{B}_{1}$ and ($\mathbb{B}_{2}$) for some
$\theta\in(0,\pi)$. Fix any $\varepsilon>0$ and let
$\gamma_{\varepsilon,\theta}$ be defined as in Proposition 5.1. Then, on
$\mathbb{D}$,

$\begin{split}G_{\phi}(x,t)\,=\,&\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\mathtt{Im}\left(\frac{e^{-x\phi\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)+t\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}}}{\rho}\right)\,d\rho+\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon,\theta}}\frac{e^{-x\phi(z)+tz}}{z}dz.\end{split}$
(5.31)

In particular, $G_{\phi}\in C^{\infty}(\mathbb{D})$, $g_{\phi}\in
C^{\infty}(\mathbb{D})$ is well defined and for any $k,l\geq 0$ we have

$\begin{split}\frac{\partial^{k}}{\partial x^{k}}\frac{\partial^{l}}{\partial
t^{l}}g_{\phi}(x,t)\,=\,&\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\mathtt{Im}\left(\left(\phi\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right)^{k}\rho^{l}e^{il\left(\pi-\frac{\theta}{2}\right)}e^{-x\phi\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)+t\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}}\right)\,d\rho\\\ &+\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon,\theta}}(\phi(z))^{k}z^{l}e^{-x\phi(z)+tz}dz.\end{split}$
(5.32)

We employ the latter together with Proposition 5.3 in the following proof.

###### Proof of Theorem 3.20.

Let us first consider $G_{\phi}$. Starting from (5.31), we have, assuming that
we can exchange the order of the series and the integral,

$\displaystyle G_{\phi}(x,t)$
$\displaystyle=\sum_{j=0}^{+\infty}(-1)^{j}\frac{x^{j}}{j!}\left[\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\mathtt{Im}\left(\frac{\left(\phi\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right)^{j}}{\rho}e^{t\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}}\right)d\rho+\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon,\theta}}\frac{(\phi(z))^{j}}{z}e^{tz}dz\right]$
$\displaystyle=\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\mathtt{Im}\left(\frac{1}{\rho}{e^{t\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}}}\right)\,d\rho+\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon,\theta}}\frac{1}{z}{e^{tz}}\,dz$
$\displaystyle+\sum_{j=1}^{+\infty}(-1)^{j}\frac{x^{j}}{j!}\left[\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\mathtt{Im}\left(\frac{\left(\phi\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right)^{j}}{\rho}e^{t\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}}\right)d\rho+\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon,\theta}}\frac{(\phi(z))^{j}}{z}e^{tz}dz\right]$
$\displaystyle={\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\mathtt{Im}\left(\frac{1}{\rho}e^{t\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}}\right)\,d\rho+\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon,\theta}}\frac{1}{z}e^{tz}\,dz}$
$\displaystyle+\sum_{j=1}^{+\infty}\sum_{k_{1}+k_{2}+k_{3}=j-1}(-1)^{j}\frac{x^{j}}{k_{1}!k_{2}!(k_{3}+1)!}q^{k_{1}}\mathfrak{b}^{k_{2}}$
$\displaystyle\times\left[\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\mathtt{Im}\left(\frac{\left(\phi^{\dagger}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right)^{k_{3}+1}}{\rho^{k_{3}+1-(k_{2}+k_{3})}e^{i(k_{3}-(k_{2}+k_{3}))}}e^{t\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}}\right)d\rho+\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon,\theta}}\frac{(\phi^{\dagger}(z))^{k_{3}+1}}{z^{k_{3}+1-(k_{2}+k_{3})}}e^{tz}dz\right]$
$\displaystyle{+\sum_{j=1}^{+\infty}(-1)^{j}\frac{x^{j}}{j!}\left[\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\mathtt{Im}\left(\frac{\left(q+\mathfrak{b}\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)^{j}}{\rho}e^{t\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}}\right)d\rho+\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon,\theta}}\frac{(q+\mathfrak{b}z)^{j}}{z}e^{tz}dz\right]}$
(5.33)

where in the last step we used the multinomial theorem, for $j\geq 1$, to
expand $\left(\phi\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right)^{j}=\left(q+\mathfrak{b}\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}+\phi^{\dagger}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right)^{j}$. Incorporating the
first summand of (5.33) into the last summation, we achieve

$\displaystyle G_{\phi}(x,t)$
$\displaystyle=\sum_{j=0}^{+\infty}(-1)^{j}\frac{x^{j}}{j!}\left[\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\mathtt{Im}\left(\frac{\left(q+\mathfrak{b}\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)^{j}}{\rho}e^{t\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}}\right)d\rho+\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon,\theta}}\frac{(q+\mathfrak{b}z)^{j}}{z}e^{tz}dz\right]$
$\displaystyle+\sum_{j=1}^{+\infty}\sum_{k_{1}+k_{2}+k_{3}=j-1}(-1)^{j}\frac{x^{j}}{k_{1}!k_{2}!(k_{3}+1)!}q^{k_{1}}\mathfrak{b}^{k_{2}}$
$\displaystyle\times\left[\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\mathtt{Im}\left(\frac{\left(\phi^{\dagger}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)\right)^{k_{3}+1}}{\rho^{k_{3}+1-(k_{2}+k_{3})}e^{i(k_{3}-(k_{2}+k_{3}))}}e^{t\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}}\right)d\rho+\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon,\theta}}\frac{(\phi^{\dagger}(z))^{k_{3}+1}}{z^{k_{3}+1-(k_{2}+k_{3})}}e^{tz}dz\right]$
$\displaystyle=\sum_{j=0}^{+\infty}(-1)^{j}\frac{x^{j}}{j!}\left[\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\mathtt{Im}\left(\frac{\left(q+\mathfrak{b}\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)^{j}}{\rho}e^{t\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}}\right)d\rho+\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon,\theta}}\frac{(q+\mathfrak{b}z)^{j}}{z}e^{tz}dz\right]$
$\displaystyle+\sum_{j=1}^{+\infty}\sum_{k_{1}+k_{2}+k_{3}=j-1}(-1)^{j}\frac{x^{j}}{k_{1}!k_{2}!(k_{3}+1)!}q^{k_{1}}\mathfrak{b}^{k_{2}}\frac{d^{k_{2}+k_{3}}}{dt}\mu^{\ast(k_{3}+1)}(t)$
$\displaystyle{=:S_{\varepsilon}(x,t;q,\mathfrak{b})+\sum_{j=1}^{+\infty}\sum_{k_{1}+k_{2}+k_{3}=j-1}(-1)^{j}\frac{x^{j}}{k_{1}!k_{2}!(k_{3}+1)!}q^{k_{1}}\mathfrak{b}^{k_{2}}\frac{d^{k_{2}+k_{3}}}{dt}\mu^{\ast(k_{3}+1)}(t)}.$

We now evaluate $S_{\varepsilon}$. Consider the function
$F_{j}(z;t)=\frac{(q+\mathfrak{b}z)^{j}}{z}e^{tz}$, that is holomorphic on
$\mathbb{C}\setminus\\{0\\}$, Hermitian and admits a simple pole in $0$ with
residue $q^{j}$. We have that

$\displaystyle\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\mathtt{Im}\left(\frac{\left(q+\mathfrak{b}\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)^{j}}{\rho}e^{t\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}}\right)d\rho+\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon,\theta}}\frac{(q+\mathfrak{b}z)^{j}}{z}e^{tz}dz$
$\displaystyle=\frac{1}{2\pi
i}\left[\int_{\varepsilon}^{+\infty}F_{j}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)};t\right)e^{i\left(\pi-\frac{\theta}{2}\right)}d\rho-\int_{\varepsilon}^{+\infty}F_{j}\left(\rho
e^{i\left(\pi+\frac{\theta}{2}\right)};t\right)e^{i\left(\pi+\frac{\theta}{2}\right)}d\rho+\int_{\gamma_{\varepsilon,\theta}}F_{j}(z;t)dz\right].$

Taking into account the notation in Figure 1, denote by $-\Gamma_{R,\theta}$
the clockwise oriented circular arc joining $F(R)$ and $C(R)$, $-\ell_{2}$ the
oriented segment joining $D(\varepsilon)$ to $C(R)$ and $-\ell_{3}$ the
oriented segment joining $F(R)$ to $E(\varepsilon)$. Denote by
$\mathfrak{D}_{R,\theta}$ the domain whose negatively oriented contour is
given by $\gamma_{\varepsilon,\theta}$, $-\ell_{2}$, $-\Gamma_{R,\theta}$ and
$-\ell_{3}$. Then, since $F_{j}$ is holomorphic on $\mathfrak{D}_{R,\theta}$,
that is a simply connected domain, we have, by Cauchy’s theorem

$\int_{-\partial\mathfrak{D}_{R,\theta}}F_{j}(z;t)dz=0.$ (5.34)

On the other hand, if we denote by $\Gamma_{R}$ and $\gamma_{\varepsilon}$ the
counterclockwise oriented circles of radius respectively $R$ and
$\varepsilon$, we have by Cauchy’s residue theorem

$\frac{1}{2\pi i}\int_{\Gamma_{R}}F_{j}(z;t)dz=\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon}}F_{j}(z;t)dz=q^{j}.$

Hence, denoting by $A_{R,\varepsilon}$ the annulus delimited by the images of
$\Gamma_{R}$ and $\gamma_{\varepsilon}$, we have

$\int_{-\partial
A_{R,\varepsilon}}F(z;t)dz=0=\int_{-\partial\mathfrak{D}_{R,\theta}}F_{j}(z;t)dz,$
(5.35)

where $-\partial A_{R,\varepsilon}$ is the negatively oriented contour of
$A_{R,\varepsilon}$. If we use the parametrization of $-\ell_{2}$ and
$\ell_{3}$ as $z=\rho e^{i\left(\pi\mp\frac{\theta}{2}\right)}$ for
$\rho\in(\varepsilon,R)$, (5.35) implies

$\int_{-\Gamma_{R}}F(z;t)dz+\int_{\gamma_{\varepsilon}}F(z;t)dz=\int_{\varepsilon}^{R}F_{j}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)};t\right)e^{i\left(\pi-\frac{\theta}{2}\right)}d\rho\\\
-\int_{\varepsilon}^{R}F_{j}\left(\rho
e^{i\left(\pi+\frac{\theta}{2}\right)};t\right)e^{i\left(\pi+\frac{\theta}{2}\right)}d\rho+\int_{\gamma_{\varepsilon,\theta}}F_{j}(z;t)dz+\int_{-\Gamma_{R,\theta}}F(z;t)dz.$
(5.36)

Furthermore, if we denote by $-\Gamma^{\dagger}_{R,\theta}$ the clockwise
oriented circular arc joining $C(R)$ to $F(R)$ it is clear that

$\int_{-\Gamma_{R}}F(z;t)dz-\int_{-\Gamma_{R,\theta}}F(z;t)dz=\int_{-\Gamma^{\dagger}_{R,\theta}}F(z;t)dz$

hence (5.36) leads to

$\int_{-\Gamma^{\dagger}_{R,\theta}}F(z;t)dz+\int_{\gamma_{\varepsilon}}F(z;t)dz\\\
=\int_{\varepsilon}^{R}F_{j}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)};t\right)e^{i\left(\pi-\frac{\theta}{2}\right)}d\rho-\int_{\varepsilon}^{R}F_{j}\left(\rho
e^{i\left(\pi+\frac{\theta}{2}\right)};t\right)e^{i\left(\pi+\frac{\theta}{2}\right)}d\rho+\int_{\gamma_{\varepsilon,\theta}}F_{j}(z;t)dz.$
(5.37)

By the estimation lemma we have

$\left|\int_{-\Gamma^{\dagger}_{R,\theta}}F(z;t)dz\right|\leq\theta
R\max_{z\in\Gamma^{\dagger}_{R,\theta}}|F(z;t)|\leq\theta(q+\mathfrak{b}R)^{j}e^{-t\cos\left(\frac{\theta}{2}\right)R}$

hence, taking the limit as $R\to+\infty$ in (5.37) we get

$2\pi iq^{j}=\int_{\gamma_{\varepsilon}}F(z;t)dz\\\
=\int_{\varepsilon}^{\infty}F_{j}\left(\rho
e^{i\left(\pi-\frac{\theta}{2}\right)};t\right)e^{i\left(\pi-\frac{\theta}{2}\right)}d\rho-\int_{\varepsilon}^{\infty}F_{j}\left(\rho
e^{i\left(\pi+\frac{\theta}{2}\right)};t\right)e^{i\left(\pi+\frac{\theta}{2}\right)}d\rho+\int_{\gamma_{\varepsilon,\theta}}F_{j}(z;t)dz.$
(5.38)

Substituting this value into $S_{\varepsilon}(x,t;q,\mathfrak{b})$ we get
(3.23).

The other two series (3.24) and (3.25) are obtained analogously, once one
notices that for any $j\geq 0$

$\frac{1}{\pi}\int_{\varepsilon}^{+\infty}\mathtt{Im}\left(\left(q+\mathfrak{b}\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}\right)^{j}e^{t\rho
e^{i\left(\pi-\frac{\theta}{2}\right)}}\right)d\rho+\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon,\theta}}(q+\mathfrak{b}z)^{j}e^{tz}dz=\frac{1}{2\pi
i}\int_{\gamma_{\varepsilon}}(q+\mathfrak{b}z)^{j}e^{tz}dz=0,$

since $(q+\mathfrak{b}z)^{j}e^{tz}$ is holomorphic in the disc
$\\{z\in\mathbb{C}:\ |z|<\varepsilon\\}$. Finally, the fact that one can
actually exchange the series with the integral is proven exactly as in the
proof of Theorem 3.15. ∎

##  6\. Examples and Further Discussion

Here we provide computational examples and further discussion of our results.

###  6.1. Examples for the results in Section 3.1

We start with examples concerning Section 3.1.

Take a Bernstein function with $q=\mathfrak{b}=0$ and a Lévy measure with a
compactly supported density that is regularly varying at $0$ with index
$-1-\alpha$ for some $\alpha\in(0,1)$, as for instance

$\mu_{\phi}(dy)=m(y)\mathbb{I}_{\\{y<\frac{1}{2}\\}}dy,\quad\,m(y)\sim
y^{-\alpha-1},\mbox{ as }x\to 0.$ (6.1)

We can opt for general regular variation and $\mathfrak{b},q>0$ but we leave
it to the interested reader. Then,

$\phi^{\prime\prime\prime}(x)\sim C_{3}x^{\alpha-3}\text{ and
}-\phi^{\prime\prime}(x)\sim C_{2}x^{\alpha-2},\mbox{ as }x\to\infty$

for positive constants $C_{2},C_{3}$. Clearly, ($\mathbb{A}_{2}$) is
satisfied. Also, due to the finite support of $\mu_{\phi}$ we have that
$\left|\phi^{\prime\prime}(0^{+})\right|<\infty,\phi^{\prime\prime\prime}(0^{+})<\infty$
and ($\mathbb{A}^{\prime}_{2}$) is valid too. Furthermore, for some $C>0$,

$x^{2}\Delta(x)=x^{2}\int_{0}^{\frac{1}{x}}y^{2}m(y)dy\sim Cx^{\alpha},\
\mbox{ as }\ x\to\infty,$

and ($\mathbb{A}_{1}$) holds with $L=\infty$. It is clear that, actually,
assumptions ($\mathbb{A}_{1}$) and ($\mathbb{A}_{2}$) hold for any Bernstein
function $\phi$ that is regularly varying at $\infty$ with index
$\alpha\in(0,1)$. Going back to the example, let $t(x)/x\downarrow 0$. Using
$-\phi^{\prime\prime}(x)\sim C_{2}x^{\alpha-2}$, as $x\to\infty$, we
establish, for any $k,l\geq 0$, as $x\to\infty$,

$\begin{split}&\sup_{(t,x)\in\mathbb{D}^{\prime}}\left|(-1)^{k}\frac{\sqrt{2\pi
C_{2}}}{C^{k+1}_{0}}\frac{\sqrt{x}}{c^{k+l+1-\frac{\alpha}{2}}}e^{-ct+x\phi(c)}\frac{\partial^{k}\partial^{l}}{\partial
x^{k}\partial
t^{l}}f_{\phi}\left(x,t\right)-1\right|=\mathrm{O}\left(\sqrt{\frac{\ln\left(x\right)}{a^{\alpha}_{*}x}}\right),\end{split}$
(6.2)

where we used Theorem 3.2 with $\phi(y)=\phi^{\dagger}(y)\sim
C_{0}y^{\alpha}$, as $y\to\infty$, and the monotonicity in $c$ of the first
asymptotic term in (3.5). Finally, employing that $\phi^{\prime}(y)\sim
C_{1}y^{\alpha-1}$, as $y\to\infty$, we get from (3.6)

$C_{1}a^{\alpha-1}_{*}\sim\frac{t(x)}{x}\quad\mbox{ as }x\to\infty\iff
a_{*}\sim\left(\frac{1}{C_{1}}\frac{t(x)}{x}\right)^{\frac{1}{\alpha-1}},\quad\mbox{
as }x\to\infty,$

and we can plug in this expression in (6.2). The speed of convergence of the
first asymptotic term in (3.5) is hence
$(t(x)/x)^{\alpha/(2(1-\alpha))}x^{-\frac{1}{2}}$ which offers faster decay
when $\alpha$ approaches $1$ and we note that it is faster than the second
term of the asymptotic in (3.5) which in this case is of order
$\mathrm{O}\left(x^{-1/2}\right)$.

We continue the example above with an illustration of Theorem 3.6, for which
only ($\mathbb{A}_{1}$) is required. Take $t(x)=x$ and then
$a_{*}=a_{*}(x)=(\phi^{\prime})^{-1}(1)$. Hence, from (3.8) we deduce that the
explicit in $x$

$\begin{split}&\frac{\partial^{k}\partial^{l}}{\partial x^{k}\partial
t^{l}}f_{\phi}\left(x,t\right)=\frac{(-1)^{k}\phi^{\dagger}(a_{*})\phi^{k}(a_{*})}{\sqrt{2\pi}a^{1-l}_{*}\sqrt{-\phi^{\prime\prime}(a_{*})}}\frac{e^{x\left(a_{*}-\phi(a_{*})\right)}}{\sqrt{x}}\left(1+\mathrm{O}\left(\sqrt{\frac{\ln(x)}{x}}\right)\right),\quad\mbox{
as }x\to\infty.\end{split}$

Since $x\phi^{\prime}(x)\leq\phi(x)$, see (4) of Lemma 2.1, we double check
that $a_{*}-\phi\left(a_{*}\right)\leq 0$.

Finally, for this example, consider the case when
$\lim\limits_{x\to\infty}t(x)/x=\infty=\phi^{\prime}(0^{+})$. Then as above
$a_{*}\sim C(t(x)/x)^{-\frac{1}{1-\alpha}}$ and
$x=\mathrm{o}\left(t(x)\right)$, as $x\to\infty$, thus

$-a^{2}_{*}\phi^{\prime\prime}(a_{*})x\sim Ca^{\alpha}_{*}x\sim
C^{\prime}t^{-\frac{\alpha}{1-\alpha}}(x)x^{\frac{1}{1-\alpha}},\quad\mbox{ as
}x\to\infty,$

which goes to infinity if $t(x)=\mathrm{o}\left(x^{\frac{1}{\alpha}}\right)$
and the first requirement of (3.9) holds. Also, for any $\delta>0$,

$\lim\limits_{x\to\infty}e^{-\delta
x}t^{-\frac{\alpha}{1-\alpha}}(x)x^{\frac{1}{1-\alpha}}=0,$

for $t(x)=\mathrm{o}\left(x^{\frac{1}{\alpha}}\right)$ and the third condition
in (3.9) holds. Under the same restrictions $\ln(1/a_{*})$ is at most of
logarithmic growth and the second imposition of (3.9) is valid and Theorem 3.8
is applicable.

Let us discuss the route to the derivation of new fine local estimates in the
region of the lower envelope. We use the example above with $\alpha=1/2$. In
this case one has to consider

$\mathbb{P}\left(\sigma(x)\leq c\frac{x^{2}}{\log_{2}x}\right),$

where $\log_{2}x=\log\log x,c>0$, see [10, Chapter III]. Set
$t(x)=c\frac{x^{2}}{\log_{2}x}$. Clearly, $t(x)/x\to\infty$ and from above
($\mathbb{A}_{1}$),($\mathbb{A}_{2}$) and ($\mathbb{A}^{\prime}_{2}$) hold.
Also,

$\phi^{\prime}(a_{*})=\frac{t(x)}{x}\Rightarrow a_{*}\sim
C\frac{\log^{2}_{2}x}{x^{2}},\quad\mbox{ as }x\to\infty.$

Furthermore, the first requirement of (3.9) is fulfilled since as
$x\to\infty$,

$-xa_{*}\phi^{\prime\prime}(a_{*})\sim Cxa^{-\frac{1}{2}}_{*}\sim
C\log_{2}(x)\to\infty.$

The second and third impositions of (3.9) are then obvious and Theorem 3.10
holds true and by virtue of its claims it yields local estimates for the
densities and all derivatives of the probabilities above.

Assume that $\left|\phi^{\prime\prime}(0^{+})\right|<\infty$ and hence
$\phi^{\prime}(0^{+})<\infty$. Let us determine the speed by which $t(x)/x$
may approach $\phi^{\prime}(0^{+})$ so that our results hold. First, we note
that since $a_{*}\to 0$ then, as $x\to\infty$,

$-a_{*}\phi^{\prime\prime}(0^{+})\sim\phi^{\prime}(0)-\phi^{\prime}(a_{*})=\phi^{\prime}(0)-\frac{t(x)}{x}.$

From Remark 3.9 we have to ensure that
$\lim\limits_{x\to\infty}a_{*}\sqrt{x}=\infty$. Hence, from the last relation
we must have

$\lim\limits_{x\to\infty}\frac{\sqrt{x}}{x\phi^{\prime}(0)-t(x)}=0.$

Thus, $\sqrt{x}=\mathrm{o}\left(x\phi^{\prime}(0)-t\right)$, as $x\to\infty$,
or alternatively $t<x\phi^{\prime}(0)-K\sqrt{x}$, for all $K>0$, and all $x$
large enough. This captures the region below that of the central limit theorem
for the density $g_{\phi}$, see Theorem 3.10, and therefore we can approximate
with high precision, as $x\to\infty,$ quantities of the type

$\mathbb{P}\left(\sigma(x)\in\left[a,b\right]\right)=\int_{a}^{b}g_{\phi}(x,t)dt.$

###  6.2. Examples for the results in Section 3.2

First we develop some explicit examples for Section 3.2. It is clear that the
representation (3.17) provides some completely explicit result if we evaluate
the derivatives of the convolution powers $\bar{\mu}_{\phi}^{\ast n}$. For
instance, this is doable for stable subordinators in which case the tail of
the Lévy measure of $\phi(z)=z^{\alpha}$ is given by

$\bar{\mu}_{\alpha}(t)=\frac{t^{-\alpha}}{\Gamma(1-\alpha)},$

where we use the subscript $\alpha$ instead of $\phi$ to underline the
dependence on the single parameter. For these specific tails, it is well-known
that

$\bar{\mu}_{\alpha}^{\ast(j+1)}(t)=\frac{t^{j-(j+1)\alpha}}{\Gamma(j+1-(j+1)\alpha)}$

and then

$\frac{d^{j}}{dt^{j}}\bar{\mu}_{\alpha}^{\ast
j}(t)=\frac{t^{-j\alpha-1}}{\Gamma(-j\alpha)}=\frac{t^{-j\alpha-1}}{\pi}\sin(\pi
j\alpha)\Gamma(1+j\alpha),$

where we have used Euler’s reflection formula, provided that
$\alpha\not\in\mathbb{Q}$, while for $\alpha\in\mathbb{Q}$ we have to pay
attention to the case in which $j\alpha$ is an integer, for which it can be
simply proven that $\bar{\mu}_{\alpha}^{\ast j}(t)$ is a monomial of degree
less than $j$ and thus $\frac{d^{j}}{dt^{j}}\bar{\mu}_{\alpha}^{\ast j}(t)=0$
as expected. For $l=0$ substituting this into (3.26) and then into (3.24), we
get that

$g_{\alpha}(x,t)=\sum_{j=1}^{+\infty}(-1)^{j+1}\frac{x^{j}}{j!}\frac{t^{-j\alpha-1}}{\pi}\sin(\pi
j\alpha)\Gamma(1+j\alpha).$ (6.3)

This series expansion is well-known in literature, see e.g. [61, Equation
$(7)$]. The same argument can be also adopted to obtain the series
representation of $f_{\alpha}$. Indeed, with the same arguments as before

$\frac{d^{j}}{dt^{j}}\bar{\mu}_{\alpha}^{\ast(j+1)}(t)=\frac{t^{-(j+1)\alpha}}{\pi(j+1)\alpha}\sin(\pi\alpha(j+1))\Gamma(\alpha(j+1)+1)$

for any $j\geq 0$ and $\alpha\in(0,1)$. Thus, by using (3.17), we get

$f_{\alpha}(x,t)=\sum_{j=0}^{\infty}(-1)^{j}\frac{\Gamma(1+(j+1)\alpha)}{\alpha(j+1)!}\frac{\sin(\pi\alpha(j+1))}{\pi}x^{j}t^{-\alpha(j+1)}.$
(6.4)

Such a series can also be deducted by combining (6.3) with the relation

$f_{\alpha}(x,t)=\frac{t}{\alpha}x^{-1-\frac{1}{\alpha}}g_{\alpha}(1,tx^{-\frac{1}{\alpha}}),$

see [50] for more details. The series representation (6.4) has also been
obtained by means of a limit argument in [44, Remark 2.3].

A similar approach can be used to deduce some information on the relativistic
(or tempered) stable subordinator, i.e. when
$\phi(z)=(\lambda+z)^{\alpha}-\lambda^{\alpha}$. In [44] the authors provide
both an integral and a series representation for the density
$f_{\alpha,\lambda}$ of the inverse tempered stable subordinator. In the
proof, they exploit the possibility to extend $\phi$ to the whole complex
half-plane $\overline{\mathbb{C}(0,\pi)}$ and then they integrate on a keyhole
contour centered in $-\lambda$. This seems to be slightly different from our
contour, which is not really a keyhole contour and it is always centered in
$0$. However, Proposition 3.19 lets us extend the approach to a full keyhole
contour. Furthermore, if we choose $\varepsilon<\lambda$, then for any $j\geq
0$ we get

$\int_{\gamma_{\varepsilon}}e^{tz}\frac{(\phi(z))^{j+1}}{z}dz=0,$

since the integrand is holomorphic on the disc $\\{z\in\mathbb{C}:\
|z|<\varepsilon\\}$. Since $\phi_{+}(\rho)$ is real for $\rho>-\lambda$, we
can rewrite

$\displaystyle\frac{d^{j+1}}{dt^{j+1}}\bar{\mu}_{\phi}^{\ast j}(t)$
$\displaystyle=\frac{1}{\pi}\int_{\lambda}^{+\infty}\mathtt{Im}\left[\frac{(\phi_{+}(-\rho))^{j+1}}{\rho}e^{-t\rho}\right]d\rho=\frac{e^{-t\lambda}}{\pi}\int_{0}^{+\infty}\mathtt{Im}\left[\frac{(\phi_{+}(\lambda-\rho))^{j+1}}{\rho+\lambda}e^{-t\rho}\right]d\rho.$

This leads to

$\displaystyle f_{\phi}(x,t)$
$\displaystyle=\frac{e^{-t\lambda}}{\pi}\sum_{j=0}^{+\infty}\frac{x^{j}}{j!}(-1)^{j}\int_{0}^{+\infty}\mathtt{Im}\left[\frac{(\phi_{+}(\lambda-\rho))^{j+1}}{\rho+\lambda}e^{-t\rho}\right]d\rho$
$\displaystyle=\frac{e^{x\lambda^{\alpha}-t\lambda}}{\pi}\sum_{j=0}^{+\infty}\frac{x^{j}}{j!}(-1)^{j}\int_{0}^{+\infty}\mathtt{Im}\left[\frac{(\phi_{+}(\lambda-\rho))^{j+1}}{\rho+\lambda}e^{-x(-\rho)^{\alpha}-t\rho}\right]d\rho$
$\displaystyle=\frac{e^{x\lambda^{\alpha}-t\lambda}}{\pi}\int_{0}^{+\infty}\frac{e^{x\rho^{\alpha}\cos(\alpha\pi)-t\rho}}{\rho+\lambda}[\rho^{\alpha}\sin(\alpha\pi-x\rho^{\alpha}\sin(i\alpha\pi))+\lambda^{\alpha}\cos(\alpha\pi)\sin(x\rho^{\alpha}\sin(i\alpha\pi))]d\rho,$

where the last equality follows by simple algebraic manipulations. It is the
integral representation of [44, Theorem 2.1] and thus, arguing as in [44,
Proposition 2.1], we get the series representation

$\begin{split}f_{\phi}(x,t)\,=\,\frac{e^{\lambda^{\alpha}x}}{\pi}\sum_{j=0}^{+\infty}\frac{(-1)^{j}x^{j}}{j!}\lambda^{\alpha(j+1)}&\left[\Gamma(1+\alpha(j+1))\Gamma(-\alpha(k+1),\lambda
t)\sin((j+1)\alpha\pi)\right.\\\ &\left.-\Gamma(1+\alpha j)\Gamma(-\alpha
j,\lambda t)\sin(j\alpha\pi)\right],\end{split}$ (6.5)

where $\Gamma(x,y)$ is the upper-incomplete Gamma function.

With the same arguments, we can use Proposition 3.19 to obtain the integral
representation of the density of the inverse Gamma subordinator (i.e. the case
$\phi(z)=\log(1+z)$), as in [45, Proposition 1].

In [14], the authors studied a special class of Thorin subordinators. Let us
consider now an example taken from this paper. Fix $\alpha\in(0,1)$, let us
consider the subordinator $\sigma$ whose Laplace exponent is given by
$\phi(z)=\varphi(z)-\varphi(0)-z$, where $\varphi(z)$ is the unique solution
of $\varphi(z)-(\varphi(z))^{\alpha}=z$ for $z\geq 0$. Clearly $\phi(0)=0$ and
then $q=0$. Furthermore, by [14, Theorem 1], we know that $\mathfrak{b}=0$
and, by [14, Proposition $5$], we get that $g_{\phi}$ is well-defined. We can
rewrite [14, Equation (21)] in the following form, for $x,t>0$,

$g_{\phi}(x,t)=\sum_{j=1}^{+\infty}\frac{x^{j}}{j!}(A_{j}(t)+B_{j}(t)+C_{j}(t)),$
(6.6)

where

$\displaystyle A_{j}(t)$
$\displaystyle=j\sum_{n=1}^{+\infty}(-1)^{n+1}\frac{\Gamma(1+n\alpha)}{n!}\sin(\pi
n\alpha)t^{-n\alpha+n}$ $\displaystyle B_{j}(t)$
$\displaystyle=j(j-1)\sum_{n=1}^{+\infty}(-1)^{n+1}\frac{(n+1)\Gamma(1+n\alpha)}{n!}\sin(\pi
n\alpha)t^{-n\alpha+n-1}$ $\displaystyle C_{j}(t)$
$\displaystyle=\begin{cases}0,&j=1,2,\\\
\displaystyle\sum_{k=0}^{j}\frac{j!}{(j-k)!}\sum_{n=k+1}^{+\infty}(-1)^{n+1}\binom{n+1}{k+2}\frac{\Gamma(1+n\alpha)}{n!}\sin(\pi
n\alpha)t^{-n\alpha+n-2-k},&j\geq 3.\end{cases}$

In particular, comparing (3.24) and (6.6), we get, by means of (3.26),

$\frac{d^{j}}{dt^{j}}\bar{\mu}_{\phi}^{\ast
j}(t)=A_{j}(t)+B_{j}(t)+C_{j}(t),\quad t>0,\ j\geq 1.$

For $j=1$, we get again [14, Equation (22)]. Moreover, setting

$\displaystyle\overline{A}_{j}(t)$
$\displaystyle=j\sum_{n=1}^{+\infty}(-1)^{n}\frac{\Gamma(1+n\alpha)}{n!(n-n\alpha+1)}\sin(\pi
n\alpha)t^{-n\alpha+n+1}$ $\displaystyle\overline{B}_{j}(t)$
$\displaystyle=j(j-1)\sum_{n=1}^{+\infty}(-1)^{n}\frac{(n+1)\Gamma(1+n\alpha)}{n!(n-n\alpha)}\sin(\pi
n\alpha)t^{-n\alpha+n}$ $\displaystyle\overline{C}_{j}(t)$
$\displaystyle=\begin{cases}0,&j=1,2,\\\
\displaystyle\sum_{k=0}^{j}\frac{j!}{(j-k)!}\sum_{n=k+1}^{+\infty}(-1)^{n}\binom{n+1}{k+2}\frac{\Gamma(1+n\alpha)}{n!(n-n\alpha-1-k)}\sin(\pi
n\alpha)t^{-n\alpha+n-1-k},&j\geq 3,\end{cases}$

we have

$\frac{d^{j}}{dt^{j}}\bar{\mu}_{\phi}^{\ast
j+1}(t)=\bar{A}_{j+1}(t)+\bar{B}_{j+1}(t)+\bar{C}_{j+1}(t),\quad t>0,\ j\geq
0,$

and then we achieve the density of the inverse subordinator $L_{\phi}$ by
means of (3.17).

Finally, we show a specific application of Theorem 3.18 to the context of
time-nonlocal equations. For example, in [2] the following assumption is used
to prove the main result:

  * •

For any $0<a\leq b$ there exists $\delta_{a,b}>0$ and a function
$F_{a,b}:(0,\delta_{a,b})\to(0,+\infty)$ such that

$\int_{0}^{\delta_{a,b}}x^{-\frac{1}{2}}F_{a,b}(x)dx<\infty$

and, for any $x\in(0,\delta_{a,b})$ and $t\in[a,b]$,

$\left|\frac{\partial f_{\phi}}{\partial x}(x,t)\right|\leq F_{a,b}(x).$

It is clear that we can Theorem 3.18 for $n=0$ and since the remainder is
locally uniform in $t$, the latter condition is verified with $F_{a,b}$ being
independent of $x$ if $\phi$ satisfies ($\mathbb{B}_{1}$) and
($\mathbb{B}_{2}$). In particular, the latter holds whenever (but not
exclusively if) $\phi$ is a complete Bernstein function, as discussed in
Section 3.2.1.

## Acknowledgements

The authors are gratefull the anonymous referees whose remarks and suggestions
have considerably improved a previous draft of the paper.

The author Giacomo Ascione has been partially supported by the MIUR-PRIN 2022
project “Anomalous Phenomena on Regular and Irregular Domains: Approximating
Complexity for the Applied Sciences”, no. 2022XZSAFN. Mladen Savov
acknowledges - “This study is financed by the European Union-NextGenerationEU,
through the National Recovery and Resilience Plan of the Republic of Bulgaria,
project No BG-RRP-2.004-0008”. The author Bruno Toaldo has been partially
supported by the MIUR-PRIN 2022 project “Non-Markovian dynamics and non-local
equations”, no. 202277N5H9. The authors Giacomo Ascione and Bruno Toaldo are
partially supported by Gruppo Nazionale per l’Analisi Matematica, la
Probabilità e le loro Applicazioni (GNAMPA-INdAM).

## References

  * [1] W. Arendt, C.J.K. Batty, M. Hieber and F. Neubrander.  Vector valued Laplace transform and Cauchy problem. Second Edition.  _Birkhäuser_ , Berlin, 2010.
  * [2] G. Ascione, P. Patie and B. Toaldo.  Non-local heat equation with moving boundary and curve-crossing of delayed Brownian motion.  _Submitted_ , available at arXiv:2203.09850, 2022.
  * [3] G. Ascione, E. Pirozzi and B. Toaldo.  On the exit time from open sets of some semi-Markov processes.  _Ann. Appl. Probab._ , 30(3): 1130 – 1163, 2020.
  * [4] G. Ascione and B. Toaldo.  A semi-Markov leaky integrate-and-fire model.  Mathematics, 7 (11): 1022, 2019.
  * [5] B. Baeumer and P. Straka. Fokker-Planck and Kolmogorov backward equations for continuous time random walk limits. _Proceedings of the American Mathematical Society_ , 145: 399 – 412, 2017.
  * [6] A. Barker and M. Savov.  Bivariate Bernstein–gamma functions and moments of exponential functionals of subordinators.  _Stochastic Process. Appl._ , 131: 454-497, 2021.
  * [7] B. Baeumer and M.M. Meerschaert. Stochastic solutions for fractional Cauchy problems. Fract. Calc. Appl. Anal. 4: 481–500, 2001.
  * [8] P. Becker Kern, M.M. Meerschaert and H. P. Scheffler. Limit theorems for coupled continuous time random walks, Ann. Probab., 32: 730 – 756, 2004.
  * [9] L. Beghin, C. Macci and C. Ricciuti. Random time-change with inverses of multivariate subordinators: Governing equations and fractional dynamics. _Stoch. Proc. Appl._ , 130(10): 6364 – 6387, 2020.
  * Bertoin [1996] J. Bertoin.  Lévy processes.  _Cambridge University Press_ , Cambridge, 1996.
  * Bertoin [1997] J. Bertoin.  Subordinators: examples and applications.  _Lectures on Probability Theory and Statistics (Saint-Flour, 1997)_ , 1 – 91. _Lectures Notes in Math._ , 1717, Springer, Berlin, 1999.
  * [12] N.H. Bingham, C. M. Goldie and J. L. Teugels.  Regular variation  Encyclopedia of Mathematics and its Applications  Cambridge University Press, Cambridge, 1987
  * [13] J. Brown, R. Churchill  Complex Variables and Applications  _McGraw Hill_ , 2014
  * [14] J. Burridge, A. Kuznetsov, M. Kwaśnicki, A. E. Kyprianou  New families of subordinators with explicit transition probability semigroup.  _Stoch. Proc. Appl._ 124(10): 3480–3495, 2014.
  * [15] Z.-Q. Chen. Time fractional equations and probabilistic representation. _Chaos, Solitons and Fractals_ , 102: 168 – 174, 2017.
  * [16] Z.-Q. Chen, P. Kim, T. Kumagai and J. Wang,  Heat kernel estimates for time fractional equations.  _Forum Mathematicum_ , 30.5: 1163–1192, 2018.
  * [17] Z.-Q. Chen, P. Kim, T. Kumagai and J. Wang  Time fractional Poisson equations: Representations and estimates.  _Journal of Functional Analysis_ , 278.2: 108311, 2020.
  * [18] S. Cho and P. Kim.  Estimates on transition densities of subordinators with jumping density decaying in mixed polynomial orders.  _Stochastic Process. Appl._ , 139: 229–279, 2021.
  * [19] E. Cinlar.  Markov additive processes and semi-regeneration.  _Discussion Paper No. 118_ , Northwestern University, 1974.
  * [20] F. Colantoni and M. D’Ovidio. On the inverse gamma subordinator. _Stoch. Anal. Appl._ , 41(5): 999 – 1024, 2023.
  * [21] R. Doney and V. Rivero.  Asymptotic behaviour of first passage time distributions for Lévy processes.  _Probab. Theor. Rel._ , 157.1-2: 1–45, 2013.
  * [22] R. Doney and V. Rivero.  Asymptotic behaviour of first passage time distributions for subordinators.  _Electron. J. Probab._ , 20: 1-28, 2015.
  * [23] L. Döring and M. Savov.  (Non)Differentiability and Asymptotics for Potential Densities of Subordinators.  _Electron. J. Probab._ , 16: 470-503, 2011.
  * [24] M. Fahrenwaldt.  Heat kernel asymptotics of the subordinator and subordinate Brownian motion.  _J. Evol. Equ._ , 19: 33–70, 2019.
  * [25] W. Feller. An Introduction to Probability Theory and Its Applications. Vol. II, 2nd Ed., Wiley, New York, 1971.
  * [26] J. E. Figueroa-López and C. Houdré.  Small-time expansions for the transition distributions of Lévy processes  _Stoc. Proc. Appl._ , 119.11: 3862 – 3889, 2009.
  * [27] J. E. Figueroa-López, R. Gong and C. Houdré.  Small-time expansions of the distributions, densities, and option prices of stochastic volatility models with Lévy jumps  _Stoc. Proc. Appl._ 122.4: 1808 – 1839, 2012.
  * [28] J. E. Figueroa-López and S. Olafsson.  Short-term asymptotics for the implied volatility skew under a stochastic volatility model with Lévy jumps  Finance and Stochastics 20.4: 973 – 1020, 2016.
  * [29] B. Fridstedt and W.E. Pruitt.  Lower functions for increasing random walks and subordinators  _Z. Wahrscheinlichkeitstheorie und Verw. Gebiete_ , 18: 167–182, 1971.
  * [30] R. Gorenflo, F. Mainardi, D. Moretti, G. Pagnini and P. Paradisi.  Discrete random walk models for space–time fractional diffusion.  _Chemical physics_ , 284(1-2): 521-541, 2002.
  * [31] R. Gorenflo and S. Vessella.  Abel Integral Equations. Analysis and Applications.  _Springer-Verlag_ , Berlin - Heidelberg, 1991.
  * [32] T. Grzywny, Ł. Leżaj and B. Trojan.  Transition densities of subordinators of positive order.  _J. Inst. Math. Jussieu_ , 1-61, 2021.
  * [33] M. Hairer, G. Iyer, L. Koralov, A. Novikov, and Z. Pajor-Gyulai. A fractional kinetic process describing the intermediate time behaviour of cellular flows. _Ann. Probab._ , 46(2): 897 – 955, 2018.
  * [34] J. M. Howie,  Complex analysis.  _Springer Science & Business Media_, 2003.
  * [35] A. Jacquier and L. Torricelli, Anomalous Diffusions in Option Prices: Connecting Trade Duration and the Volatility Term Structure,  _SIAM J. Financial Math.,_ 11(4), 1137–1167, 2020.
  * [36] N. C. Jain and W. E. Pruitt.  Lower tail probability estimates for subordinators and non-decreasing random walks  _Ann. Probab.,_ 15: 75 – 101, 1987.
  * [37] H. Kaspi and B. Maisonneuve.  Regenerative systems on the real line.  _Ann. Probab._ , 16: 1306 – 1332, 1988.
  * [38] V. Knopova and A. Kulik.  Intrinsic small time estimates for distribution densities of Lévy processes.  _Random Operators and Stochastic Equations_ , 21.4: 321–344, 2013.
  * [39] K. Kobayashi and H. Park  Large-time and small-time behaviors of the spectral heat content for time-changed stable processes.  _Electronic Communications in Probability_ , 27: 1–11, 2022.
  * [40] K. Kobayashi and H. Park.  Spectral heat content for time-changed killed Brownian motions.  _Journal of Theoretical Probability_ , 36.2: 1148–1180, 2023.
  * [41] A.N. Kochubei, Y.G. Kondratiev and J.L. da Silva.  From Random Times to Fractional Kinetics.  _Interdisciplinary Studies of Complex Systems,_ 16: 5 – 32, 2020.
  * Kochubei et al. [2020] A.N. Kochubei, Y.G. Kondratiev and J.L. da Silva.  Random time change and related evolution equations. Time asymptotic behaviour.  _Stochastics and Dynamics_ , 20(5): 2050034, 2020.
  * [43] M. Kovacs and M.M. Meerschaert. Ultrafast subordinators and their hitting times. Publications de L’Institut Mathematique, 94(71): 193 – 206, 2006.
  * [44] A. Kumar and P. Vellaisamy. Inverse tempered stable subordinators. _Statistics and Probability Letters_ , 103: 134 – 141, 2015.
  * [45] A. Kumar, A. Wyłomańska, A., R. Połoczański, and S. Sundar  Fractional Brownian motion time-changed by gamma and inverse gamma process.  _Physica A: Statistical Mechanics and Its Applications_ , 468: 648–667, 2017.
  * [46] A.E. Kyprianou.  Fluctuations of Lévy Processes with Applications.  _Springer-Verlag Berlin Heidelberg_ , 2014.
  * Meerschaert and Scheffler [2004] M.M. Meerschaert and H.P. Scheffler.  Limit theorems for continuous-time random walks with infinite mean waiting times.  _J. Appl. Probab._ , 41.3: 623–638, 2004.
  * Meerschaert and Scheffler [2008] M.M. Meerschaert and H.P. Scheffler.  Triangular array limits for continuous time random walks.  _Stochastic Process. Appl._ , 118(9): 1606 – 1633, 2008.
  * Meerschaert and Sikorskii [2012] M.M. Meerschaert and A. Sikorskii. Stochastic Models for Fractional Calculus. De Gruyter, Berlin, 2012.
  * MeerschaertStraka2013 [2013] M. M. Meerschaert, and P. Straka.  Inverse stable subordinators.  _Mathematical modelling of natural phenomena_ , 8.2: 1 – 16, 2013.
  * Meerschaert and Straka [2014] M.M. Meerschaert and P. Straka.  Semi-Markov approach to continuous time random walk limit processes.  _Ann. Probab._ , 42(4) : 1699 – 1723, 2014.
  * [52] R. Metzler and J. Klafter. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. _Phys. Rep._ , 339: 1 – 77, 2000.
  * [53] R. Metzler, J.-H. Jeon, A.G. Chertsvy and E. Barkai. Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking. _Phys. Chem. Chem. Phys._ , 16: 24128, 2014.
  * [54] M. Minchev and M. Savov.  Asymptotics for densities of exponential functionals of subordinators.  _Bernoulli_ , 29(4): 3307–3333, DOI: 10.3150/23-BEJ1584, 2023
  * [55] A. Mura and G. Pagnini. Characterizations and simulations of a class of stochastic processes to model anomalous diffusion. _J. Phys. A.: Math. Theor._ , 41: 285003, 2008.
  * [56] F.Olver. Error bounds for the Laplace approximation for definite integrals. _SIAM Review_ , 1: 293–313, 1968.
  * [57] F.Olver. Why steepest descents? _SIAM Review_ , 12(2): 228–248, 1970.
  * [58] P. Patie and M. Savov.  Bernstein-gamma functions and exponential functionals of Lévy processes.  _Electron. J. Probab._ , 23(75): 1 – 101, 2018.
  * [59] P. Patie and M. Savov.  Spectral expansions for generalized Laguerre semigroups.  _Mem. Amer. Math. Soc._ , 272(1336): 1 – 183, 2021.
  * [60] P. Patie and A. Vaidyanathan.  Non-classical Tauberian and Abelian type criteria for the moment problem.  _Math. Nachr._ , 295(5): 970 - 990, 2022.
  * [61] K. A. Penson and K. Górska.  Exact and explicit probability densities for one-sided Lévy stable distributions.  Phys. Rev. Lett. 105.21: 210604, 2010.
  * [62] C. Ricciuti and B. Toaldo. From semi-Markov random evolutions to scattering transport and superdiffusion. Commun. Math. Phys., 401(3): 2999 – 3042, 2023.
  * [63] M. Savov and B. Toaldo. Semi-Markov processes, integro-differential equations and anomalous diffusion-aggregation. Ann. Inst. Henri Poincare Probab. Stat., 56(4): 2640 – 2671, 2020.
  * [64] E. Scalas. Five years of continuous-time random walks in econophysics. The Complex Networks of Economic Interactions, _Lecture Notes in Econom. and Math. Systems_ , 567: 3–16, Springer, Berlin, 2006.
  * [65] E. Scalas and B. Toaldo. Limit theorems for prices of options written on semi-Markov processes. _Theory Probab. Math. Statist._ , 105: 3 – 33, 2021.
  * [66] R.L. Schilling, R. Song and Z. Vondraček.  Bernstein functions: theory and applications. 2nd edition.  _Walter de Gruyter GmbH & Company KG_, Vol 37 of De Gruyter Studies in Mathematics Series, 2010.
  * [67] J.L. da Silva and Y.G. Kondratiev.  Random Time Dynamical Systems I: General Structures.  _https://doi.org/10.48550/arXiv.2012.15201_ , 2020.
  * [68] E. M. Stein and R. Shakarchi.  Complex analysis. Vol. 2  Princeton University Press, 2010.
  * [69] L. Torricelli. Volatility targeting using delayed diffusions. _Appl. Math. Finance_ , 25(3): 213 – 246, 2018.
  * [70] S. Umarov, M. Hahn and K. Kobayashi. Beyond the Triangle: Brownian Motion, Itô Calculus, and Fokker-Planck Equation-Fractional Generalizations. _Wold Scientific_ , 2018.
  * [71] M. Veillette and M.S. Taqqu.  Using differential equations to obtain joint moments of first-passage times of increasing Lévy processes.  _Statistics and Probability Letters_ , 80: 697 – 705, 2010.


Inverse extremal problem for an anti-tumor
therapy model

Andrey Kovtanyuk1,$\star$, Christina Kuttler1, Kristina Koshel2,

and Alexander Chebotarev2

1 Technical University of Munich, Mathematical Department, Boltzmannstr. 3,
85748, Garching, Germany

2 Far Eastern Federal University, AjaxBay 10, Russky Island,

690922, Vladivostok, Russia

email: $\star$<EMAIL_ADDRESS>
###### Abstract.

An optimal control problem for a model of tumor growth is studied. In a given
subdomain, it is required to minimize the density of tumor cells, while the
drug concentration in tissue is limited by given minimal and maximal values.
Based on derived estimates of the solution of the controlled system, the
solvability of the control problem is proved. The problem is reduced to an
optimal control problem with a penalty. An algorithm for solving the optimal
control problem with a penalty is constructed and implemented. The efficiency
of the algorithm is illustrated by a numerical example.

##  1\. Introduction

Currently, there are many mathematical models of tumor evolution based on ODEs
or PDEs [1, 2, 3, 4]. Often such models are quite complex and include several
state variables. For example, the article [2] considers a model that contains
five state variables, two of which are volume fractions of proliferative and
necrotic tumor cells. A 3D–1D model studied in [3] includes seven state
variables of which three relate to volume fractions of different species of
tumor cells (proliferative, hypoxic, and necrotic cells). A model considered
in [4] contains eleven state variables including three species of tumor cells.
Note that the models from [2, 3, 4] also take into account the influence of
anti-tumor drug therapy.

Some studies of tumor evolution are based on reaction-diffusion models in
which tumor cells are not divided into species (see, e.g., [5, 6, 7, 8], where
tumor cell growth simulation results were presented). Many diffusion models on
tumor cell growth and death contain nonlinear reaction terms. Note that the
equation for glioma cell growth in [6], in addition to the nonlinear reaction
term, also contains a nonlinear diffusion term.

Optimal control of the volume fraction of tumor cells under the influence of
drug therapy is important for assessing the effectiveness of the treatment
plan and its correction. There are a series of works on the optimal control of
tumor treatment that do not take into account its spatial distribution (see,
e.g., [9, 10]). Optimal control problems for diffusion models of tumor
treatment are studied in [7, 11].

In the current paper, a model of tumor growth under the influence of drug
therapy is under consideration. The model describes the behavior of two state
variables: the normalized density of tumor cells and normalized drug
concentration. An optimal control problem for the considered model is studied.
Following this, in a given subdomain it is necessary to minimize the density
of tumor cells, while the drug concentration is limited to the specified
minimum and maximum values. The solvability of the optimal control problem is
proved. An algorithm for solving the optimal control problem by minimizing the
objective functional with a penalty is proposed. The efficiency of the
algorithm is illustrated by a numerical example.

##  2\. Formalization of the optimal control problem

We consider the following initial-boundary value problem describing the growth
and death of tumor cells in the domain $\Omega$ with boundary $\Gamma$ during
time interval $(0,T)$:

(1) $y^{\prime}-\nabla\cdot(k\nabla
y)=d(s)y,\;\;\;y|_{\Gamma}=0,\;\;\;y(x,0)=y_{0};$ (2)
$s^{\prime}(t)+M_{0}s(t)=u(t),\;\;\;s(0)=0,$

where $y$ is the normalized density of tumor cells, $s$ is the normalized drug
concentration in tissue, $u$ describes the drug supply, $k=k(x)$ is the
diffusion coefficient, $M_{0}$ is positive constants.

Let $\Omega$ be a Lipschitz bounded domain, $Q=\Omega\times(0,T)$,
$\Sigma=\Gamma\times(0,T)$. We denote by $L^{p}$, $1\leq p\leq\infty$, the
Lebesgue space, by $H^{1}$ the Sobolev space $W^{1}_{2}$, and by
$L^{p}(0,T;X)$ the Lebesgue space of functions from $L^{p}$, defined on
$(0,T)$, with values in Banach space $X$. Also let
$H=L^{2}(\Omega),\;V=H^{1}_{0}(\Omega)=\\{v\in
H^{1}(\Omega):\;v|_{\Gamma}=0\\}$, and the space $V^{\prime}$ be the dual to
$V$. Then we identify $H$ with its dual space $H^{\prime}$ such that $V\subset
H^{1}(\Omega)\subset H=H^{\prime}\subset(H^{1}(\Omega))^{\prime}\subset
V^{\prime}$, and denote by $\|\cdot\|$ the norm in $H$, and by $(h,v)$ the
value of functional $h\in V^{\prime}$ on the element $v\in V$ coinciding with
the inner product in $H$ if $h\in H$. We also define the control space
$U=L^{2}(0,T)$ and the state space

$Y=\\{y\in L^{2}(0,T;V)\colon y^{\prime}=dy/dt\in L^{2}(0,T,V^{\prime})\\}.$

We assume that the input data satisfies the following conditions:
$(i)\;\;\;0<k_{0}\leq k(x)\leq k_{1},\;\;k_{j}=Const$, $y_{0}\in H$;
$(ii)\;\;d(s)(s-s_{m})<0,\;\;s\neq s_{m}$, $|d(s_{1})-d(s_{2})|\leq
L_{a}|s_{1}-s_{2}|\;\;\forall\;|s_{1,2}|\leq a.$

Let us define the operator $A\colon V\to V^{\prime}$ using the equality:

$\forall\;y,v\in V\quad(Ay,v)=(k\nabla y,\nabla v).$

Therefore, for $y\in Y$, the problem (1) can be rewritten in the following
weak form:

(3) $y^{\prime}+Ay=d(s)y,\;\;\;y(0)=y_{0}.$

Here and further, when writing differential equations, we assume that they
hold over the entire time interval $(0,T)$.

Problem ($P$) For $\lambda>0$, it is required to minimize the objective
functional

(4)
$J(y,u)=\frac{1}{2}\|y\|^{2}_{L^{2}(Q)}+\frac{\lambda}{2}\|u\|^{2}_{L^{2}(0,T)}\rightarrow\inf$

on functions $y\in Y$, $u\in U$ satisfying the conditions (2) and (3) such
that

$s(t)\leq s_{+},\;\;t\in(0,T);\;\;s(t)\geq s_{-},\;\;t\in(t_{0},T),$

where $s_{-}$, $t_{0}$ are positive constants.

##  3\. Solvability of the optimal control problem

First, consider the properties of problems (2) and (3). Define the operator
$B:L^{2}(0,T)\to H^{1}(0,T)$ so that $s=B(u)$ if $s$ is a solution of problem
(2).

######  Lemma 1.

For $s=B(u)$ the following estimates hold:

(5) $|s(t)|\leq\sqrt{t}\|u\|_{U},\;\;\;\|s\|_{U}\leq
T\|u\|_{U},\;\;\;\|s^{\prime}\|_{U}\leq(1+M_{0}T)\|u\|_{U}.$

###### Proof.

The first two estimates (5) immediately follow from the representation of the
solution of problem (2),

$s(t)=\int_{0}^{t}e^{-M_{0}(t-\tau)}u(\tau)d\tau.$

The inequality for $\|s^{\prime}\|_{U}$ follows from the equation (2) and the
estimate for $\|s\|_{U}$. ∎

######  Lemma 2.

Let conditions (i), (ii) hold. Then for $u\in U$ there is a unique solution
$y\in Y$ of problem (3), where $s=B(u)$.

###### Proof.

Let $\chi=d(s)$. By Lemma 1, $|s|\leq a$, where $a=\sqrt{T}\|u\|_{U}$, and
therefore $|\chi|\leq|d(0)|+|d(s)-d(0)|\leq|d(0)|+L_{a}a=a_{1}$. Problem (3)
can be written as

$y^{\prime}+A_{1}y=0,\;\;\;y(0)=y_{0},$

where $A_{1}y=Ay-\chi y$. Note that

$(A_{1}z,z)+a_{1}\|z\|^{2}\geq(k\nabla z,\nabla z)\geq k_{0}\|z\|^{2}_{V}.$

Therefore [12, Ch.3, Th.1.2] there is a unique solution $y\in Y$ of problem
(3). ∎

######  Theorem 1.

Let conditions (i), (ii) hold and

(6) $\frac{1-e^{-M_{0}T}}{1-e^{-M_{0}t_{0}}}s_{-}\leq s_{+}.$

Then there is a solution of problem (P).

###### Proof.

Let us show that the set of controls ensuring the fulfillment of the
inequalities $s(t)\leq s_{+},\;t\in[0,T]$, $s(t)\geq s_{-},\;t\in[t_{0},T]$,
is non-empty. Indeed, if we take $u:=M_{0}(1-e^{-M_{0}t_{0}})^{-1}s_{-}$ as a
control, then the solution to problem (2) has the form
$s(t)=u(1-e^{-M_{0}t})/M_{0}$ and then $s(t)\geq s_{-},\;t\in[t_{0},T]$.
Condition (6) guarantees that $s(t)\leq s_{+},\;t\in[0,T]$.

Let $J(y_{m},u_{m})\to j=\inf J$,

(7)
$y_{m}^{\prime}+Ay_{m}=d(s_{m})y_{m},\;\;\;s_{m}=B(u_{m}),\;\;\;y_{m}(0)=y_{0};$
$s_{m}(t)\leq s_{+},\;\;\;t\in[0,T],\;\;\;s_{m}(t)\geq
s_{-},\;\;\;t\in[t_{0},T].$

The structure of $J$ implies the estimate $\|u_{m}\|_{U}\leq C$ and therefore,
by Lemma 1,

$|s_{m}(t)|\leq C,\;\;\;\|s_{m}\|_{H^{1}(0,T)}\leq C,\;\;\;|d(s_{m}(t))|\leq
C.$

Here and below, by $C>0$, we denote constants independent of $m$. Multiplying
in the sense of the inner product the equation for $y_{m}$ in (7) by $y_{m}$
and taking into account the estimates of $s_{m}$, we derive the inequality

$\frac{1}{2}\frac{d}{dt}\|y_{m}\|^{2}+k_{0}\|\nabla y_{m}\|^{2}\leq
C\|y_{m}\|^{2}.$

Integrating this inequality over $t$ and applying Gronwall’s lemma, we obtain
the estimate $\|y_{m}\|_{Y}\leq C.$ Based on the estimates obtained, we
conclude, passing to subsequences if necessary, that there exist functions
$u\in U,\;s\in H^{1}(0,T),\;y\in Y$ such that

$u_{m}\to u\,\text{ weakly in }U,\;\;\;s_{m}\to s\,\text{ weakly in
}H^{1}(0,T),\text{ strongly in }U,$ $y_{m}\to y\,\text{ weakly in
}L^{2}(0,T;V),\;\;\text{ strongly in }L^{2}(0,T;H).$

Moreover, the Lipschitz property of $d$ guarantees that $d(s_{m})\to d(s)$ in
$U=L^{2}(0,T)$. The obtained results on convergences allow us to pass to the
limit in (7). As a result,

$y^{\prime}+Ay=d(s)y,\;\;\;s=B(u),\;\;\;y(0)=y_{0};$ $s(t)\leq
s_{+},\;\;t\in[0,T];\;\;\;s(t)\geq s_{-},\;\;t\in[t_{0},T].$

Since the objective functional is weakly lower semicontinuous, then $j\leq
J(y,u)\leq\liminf J(y_{m},u_{m})=j$ and therefore the pair $\\{y,u\\}$ is a
solution of problem ($P$). ∎

##  4\. Problem with penalty

Let us define the constraint operator

$F:W\times H^{1}(0,T)\times U\to L^{2}(0,T;V^{\prime})\times U\times
H\times\mathbb{R}$

such that

$F(y,s,u)=\\{y^{\prime}+Ay-d(s)y,s^{\prime}+M_{0}s-u,y(0)-y_{0},s(0)\\}.$

Further, we consider the following extremal problem with a penalty.

Problem ($P_{\varepsilon}$)

$J_{\varepsilon}(y,s,u)=J(y,u)+\frac{1}{\varepsilon}\int\limits_{0}^{T}f_{1}(s(t))dt+\frac{1}{\varepsilon}\int\limits_{t_{0}}^{T}f_{2}(s(t))dt\rightarrow\inf,\\\
F(y,s,u)=0,\;\;u\in U.$

Here, $\varepsilon>0$,

$f_{1}(s)=\begin{cases}0,&\text{if $s\leq s_{+}$,}\\\ (s-s_{+})^{2},&\text{if
$s>s_{+}$,}\end{cases}$ $f_{2}(s)=\begin{cases}0,&\text{if $s\geq s_{-}$,}\\\
(s-s_{-})^{2},&\text{if $s<s_{-}$.}\end{cases}$

Estimates for the solution of the controllable system make it possible to
prove the solvability of the problem with a penalty ($P_{\varepsilon}$)
similar to the proof of Theorem 1.

######  Theorem 2.

Let conditions (i), (ii) hold. Then a solution of the problem
($P_{\varepsilon}$) exists.

Consider the approximation properties of solutions of the problem with a
penalty.

######  Theorem 3.

Let conditions (i), (ii) hold. If
$\\{y_{\varepsilon},s_{\varepsilon},u_{\varepsilon}\\}$ is a solution of the
problem ($P_{\varepsilon}$) for $\varepsilon>0$, then there is a sequence
$\varepsilon\to+0$ such that

$u_{\varepsilon}\rightarrow\widehat{u}\;\text{ weakly in
}U,\;\;\;y_{\varepsilon}\rightarrow\widehat{y}\;\text{ strongly in
}L^{2}(0,T;H),\\\ s_{\varepsilon}\rightarrow\widehat{s}\;\text{ strongly in
}U,$

where $\\{\widehat{y},\,\widehat{u}\\}$ is a solution of Problem ($P$),
$\widehat{s}=B(\widehat{u})$.

###### Proof.

Consider the pair $\\{y,u\\}$, a solution to problem ($P$), the existence of
which is guaranteed by Theorem 1. Since $s=B(u)$ satisfies the inequalities
$s(t)\leq s_{+},\;t\in[0,T]$; $s(t)\geq s_{-},\;t\in[t_{0},T]$ that means
$f_{1}(s(t))=0,\;t\in[0,T]$, $f_{2}(s(t))=0,\;t\in[t_{0},T]$, then

$J(y_{\varepsilon},u_{\varepsilon})+\frac{1}{\varepsilon}\int\limits_{0}^{T}f_{1}(s_{\varepsilon}(t))dt+\frac{1}{\varepsilon}\int\limits_{t_{0}}^{T}f_{2}(s_{\varepsilon}(t))dt\\\
\leq J(y,u).$

Here, $s_{\varepsilon}=B(u_{\varepsilon})$. Therefore

(8) $J(y_{\varepsilon},u_{\varepsilon})\leq
J(y,u),\;\;\;\int\limits_{0}^{T}f_{1}(s_{\varepsilon}(t))dt+\int\limits_{t_{0}}^{T}f_{2}(s_{\varepsilon}(t))dt\leq\varepsilon
J(y,u).$

From the first inequality (8) it follows that $\|u_{\varepsilon}\|_{U}\leq C$
and similarly to the proof of Theorem 1 we derive the estimates

$|s_{\varepsilon}(t)|\leq C,\;\;\;\|s_{\varepsilon}\|_{H^{1}(0,T)}\leq
C,\;\;\;|d(s_{\varepsilon}(t))|\leq C,\;\;\;\|y_{\varepsilon}\|_{Y}\leq C.$

Here, by $C>0$, we denote constants independent of $\varepsilon$. Based on the
estimates obtained, we conclude, passing to subsequences if necessary, that
there are functions $\widehat{u}\in U,\;\widehat{s}\in
H^{1}(0,T),\;\widehat{y}\in Y$ such that

$u_{\varepsilon}\to\widehat{u}\;\,\text{ weakly in
}U,\;\;s_{\varepsilon}\to\widehat{s}\;\,\text{ weakly in }H^{1}(0,T),\text{
strongly in }U,$ $y_{\varepsilon}\to\widehat{y}\;\,\text{ weakly in
}L^{2}(0,T;V),\;\,\text{ strongly in }L^{2}(0,T;H).$

At the same time $d(s_{\varepsilon})\to d(\widehat{s})$ in $U=L^{2}(0,T)$. The
obtained results on convergence allow us, by passing to the limit in the
equality $F(y_{\varepsilon},s_{\varepsilon},u_{\varepsilon})=0$, to arrive at
the relation $F(\widehat{y},\widehat{s},\widehat{u})=0$. Moreover,

$\int\limits_{0}^{T}f_{1}(s_{\varepsilon}(t))dt\to\int\limits_{0}^{T}f_{1}(\widehat{s}(t))dt,\;\;\;\int\limits_{t_{0}}^{T}f_{2}(s_{\varepsilon}(t))dt\to\int\limits_{t_{0}}^{T}f_{2}(\widehat{s}(t))dt,$

and since by virtue of (8) these limits equal to $0$, then $\widehat{s}(t)\leq
s_{+},\;t\in[0,T]$; $\widehat{s}(t)\geq s_{-},\;t\in[t_{0},T].$ Therefore, the
pair $\\{\widehat{y},\widehat{u}\\}$ is admissible for problem ($P$). Notice
that

$j=\inf J\leq J(\widehat{y},\widehat{u})\leq\liminf
J(y_{\varepsilon},u_{\varepsilon})\leq J(y,u)=j,$

and therefore $\\{\widehat{y},\widehat{u}\\}$ is a solution of problem ($P$).
∎

To obtain the optimality system, we will use the Lagrange principle for
smooth-convex extremal problems [13]. Let us establish the non-degeneracy of
the optimality conditions by showing that the image of the derivative of the
constraint operator with respect to $s$ coincides with the space
$W=L^{2}(0,T;V^{\prime})\times U\times H\times\mathbb{R}$.

######  Lemma 3.

Let conditions (i), (ii) hold. For any $\widehat{y}\in Y$, $\widehat{s}\in
H^{1}(0,T)$, $\widehat{u}\in U$ the following equality is true:

$\texttt{Im}\,F^{\prime}_{\nu}(\widehat{y},\widehat{s},\widehat{u})=W.$

Here, $\nu=\\{y,s\\}$.

###### Proof.

Consider an arbitrary element $\eta=\\{\zeta,q,z_{0},h_{0}\\}\in W$ and check
that there is a solution $z\in Y$, $h\in H^{1}(0,T)$ of the system

$z^{\prime}+Az-d(\widehat{s})z=\zeta+d^{\prime}(\widehat{s})h\widehat{y},\;\;z(0)=z_{0},$
$h^{\prime}+M_{0}h=q,\;\;h(0)=h_{0}.$

The existence of a solution $h\in H^{1}(0,T)$ of the Cauchy problem for the
linear ordinary differential equation is obvious. The solvability of the
parabolic problem for $z$ with the right-hand side
$\zeta+d^{\prime}(\widehat{s})h\widehat{y}\in L^{2}(0,T;V^{\prime})$ follows
from [12, Ch.3, Th. 1.2]. ∎

######  Theorem 4.

Let conditions (i), (ii) hold. If
$\\{y_{\varepsilon},s_{\varepsilon},u_{\varepsilon}\\}$ is a solution of
problem ($P_{\varepsilon}$) for $\varepsilon>0$, then there is a unique
solution $\\{p_{1},p_{2}\\}\in W\times H^{1}(0,T)$ of the dual system

(9)
$-p_{1}^{\prime}+Ap_{1}-d(s_{\varepsilon})p_{1}=-y_{\varepsilon},\;\;\;p_{1}(T)=0;$
(10)
$-p_{2}^{\prime}+M_{0}p_{2}+\frac{1}{\varepsilon}f_{1}^{\prime}(s_{\varepsilon})+\frac{1}{\varepsilon}\chi(s_{\varepsilon},t)\\\
=d^{\prime}(s_{\varepsilon})(y_{\varepsilon},p_{1}),\;\;\;t\in(0,T),\;\;\;p_{2}(T)=0,$

where $\chi(s,t)=f_{2}^{\prime}(s)$ if $t\in[t_{0},T]$ and $\chi(s,t)=0$ if
$t\in[0,t_{0})$. Moreover, the parameter $u_{\varepsilon}$ is defined by the
following equality $u_{\varepsilon}=\lambda^{-1}p_{2}.$

###### Proof.

The statement of the theorem follows from [13, Chapter 2, Theorem 1.5]. By
Lemma 3, the Lagrange function of problem ($P$) is defined by the equality

$L(y,s,u,p_{1},p_{2},q_{1},q_{2})=J_{\varepsilon}(y,s,u)+\int_{0}^{T}(y^{\prime}+Ay-d(s)y,p_{1})dt+\\\
\int_{0}^{T}(s^{\prime}+M_{0}s-u)p_{2}dt+(y(0)-y_{0},q_{1})+s(0)q_{2},$

where $\\{p_{1},p_{2}\\}\in Y\times H^{1}(0,T),\;\,q_{1}\in H,\;\,\text{and}$
$q_{2}\in\mathbb{R}$. In accordance with Lagrange’s principle, equalities
$L_{y}^{\prime}(y_{\varepsilon},s_{\varepsilon},u_{\varepsilon},p_{1},p_{2},q_{1},q_{2})=0$
and
$L_{s}^{\prime}(y_{\varepsilon},s_{\varepsilon},u_{\varepsilon},p_{1},p_{2},q_{1},q_{2})=0$
give the conjugate system (9), (10), and the condition
$L_{u}^{\prime}(y_{\varepsilon},s_{\varepsilon},u_{\varepsilon},p_{1},p_{2},q_{1},q_{2})=0$
implies the equality $u_{\varepsilon}=\lambda^{-1}p_{2}.$ ∎

##  5\. Numerical algorithm

1\. Choose the value of the gradient step $\delta$, the number of iterations
$N$, the value of the penalty coefficient $\varepsilon$, and the initial
approximation for the control $u_{0}\in U$.

2\. Loop for $k=0,1,2,\dots,N$.
For a given $u_{k}$, calculate the state $y_{k},\,s_{k}$ that is the solution
of the problem (2), (3).

3\. Calculate the value of the objective functional
$J_{\varepsilon}(y_{k},s_{k},u_{k})$.

4\. Calculate the conjugate state $\\{p_{1k},\,p_{2k}\\}$ from the equations
(9), (10) setting $y_{\varepsilon}=y_{k},\;s_{\varepsilon}=s_{k}$.

5\. Recalculate the control $u_{k+1}=u_{k}-\delta(\lambda u_{k}-p_{2k})$.

The value of the parameter $\delta$ is chosen empirically in such a way that
the value $\delta(\lambda u_{k}-p_{2k})$ is a significant correction for
$u_{k+1}$. The number of iterations $N$ is chosen sufficient to satisfy the
condition
$J_{\varepsilon}(y_{k},s_{k},u_{k})-J_{\varepsilon}(y_{k+1},s_{k+1}u_{k+1})<\epsilon$,
where $\epsilon>0$ determines the accuracy of calculations.

##  6\. Numerical example

As the computational domain, we consider a square with an edge of 3 cm. The
tumor occupies the circle located in the center of the computational domain
with a diameter of 1 cm. We set the following values of the problem
parameters: $s_{c}=0.2$, $s_{-}=0.4$, $s_{+}=0.8$, $T=28$ (days). The value of
the diffusion coefficient $k$ is $2.5\cdot 10^{-9}$ (cm2/s) [14]. The initial
distribution of the control is set as follows: it equals 0.00014 (s-1) during
the first hour of each day. The value of the penalty coefficient $\varepsilon$
is set equal to 0.2. The number of iterations is 10.

The implementation of the computational algorithm was fulfilled by the finite
element method package FreeFEM ++ [15]. In Figure 1, the behavior of the
objective functional $J_{\varepsilon}$ in dependence on the number of
iterations is shown.

In Figure 2, the behavior of the control $u_{\varepsilon}$ during 28 days of
observation is shown. The dots indicate the amount of medication taken every 8
hours. The red plot is their continuous approximation. In Figure 3, the
behavior of the drug concentration in tissue $s_{\varepsilon}$ (blue plot) in
comparison with the minimal and maximal constraints ($s_{-}$ and $s_{+}$, red
lines) and the critical level of drug concentration ($s_{c}$, violet line) is
shown. In Figure 4, the distributions of tumor cell density $y_{\varepsilon}$
in the central cross-section of the computational domain for different time
moments (initial time moment, 1 week, 2 weeks, 3 weeks, and 4 weeks) are
shown. As can be seen, during the whole time interval, a monotonic decrease in
tumor cell density is observed.

##  7\. Conclusion

The optimal control problem for the model of tumor evolution under the
influence of drug therapy has been studied. The solvability of the problem has
been proved. The numerical algorithm for finding a solution has been
constructed, its convergence has been established.

Subsequent efforts of the authors will be aimed at studying a model of tumor
evolution that includes a nonlinear reaction term and describes the behavior
of two or more species of tumor cells.

Figure 1: Objective functional.

Figure 2: Control (dots) and its continuous approximation (red plot).

Figure 3: Drug concentration (blue plot), constraints on drug concentration
(red lines), and the critical value of drug concentration (violet line).

Figure 4: Density of tumor cells in different time moments: initial
distribution (blue), 1 week (red), 2 weeks (black), 3 weeks (green), and 4
weeks (yellow).

## References

  * [1] Schönfeld, S., Ozkan, A., Scarabosio, L., Rylander, M.N., Kuttler, C., 2023, Environmental stress level to model tumor cell growth and survival, _Math. Biosci. Eng._ , Vol. 19(6), pp. 5509–5545.
  * [2] Wagner, A., Schlicke, P., Fritz, M., Kuttler, C., Oden, J.T., Schumann, C., Wohlmuth, B., 2023, A phase-field model for non-small cell lung cancer under the effects of immunotherapy, _Math. Biosci. Eng._ , Vol. 20(10), pp. 18670–18694.
  * [3] Fritz, M., Jha, P.K., Köppl, T., Oden, J.T., Wohlmuth, B., 2021, Analysis of a new multispecies tumor growth model coupling 3D phase-fields with a 1D vascular network, _Nonlinear Anal. Real World Appl._ , Vol. 61, n.a. 103331.
  * [4] Kuznetsov, M.B., Kolobov, A.V., 2018, Transient alleviation of tumor hypoxia during first days of antiangiogenic therapy as a result of therapy-induced alterations in nutrient supply and tumor metabolism – Analysis by mathematical modeling, _J. Theor. Biol._ , Vol. 451, pp. 86–100.
  * [5] Harpold, H.L.P., Alvord, E.C., Swanson, K.R., 2007, The evolution of mathematical modeling of glioma, _J. Neuropathol. Exp. Neurol._ , Vol. 66(1), pp. 1–9.
  * [6] Alfonso, J., Köhn-Luque, A., Stylianopoulos, T., Feuerhake, F., Deutsch, A., Hatzikirou, H., 2016, Why one-size-fits-all vaso-modulatory interventions fail to control glioma invasion: in silico insights, _Sci. Rep._ , Vol. 6, n.a. 37283.
  * [7] Colli, P., Gomez, H., Lorenzo, G., Marinoschi, G., Reali, A., Rocca, E., 2021, Optimal control of cytotoxic and antiangiogenic therapies on prostate cancer growth, _Math. Mod. Meth. Appl. Sci._ , Vol. 31(7), pp. 1419–1468.
  * [8] Viguerie, A., Grave, M., Barros, G.F., Lorenzo, G., Reali, A., Coutinho, A.L.G.A., 2022, Data-driven simulation of Fisher–Kolmogorov tumor growth models using dynamic mode decomposition, _J. Biomech. Eng._ , Vol. 144(12), n.a. 121001.
  * [9] Schättler, H., Ledzewicz, U., 2015, Optimal control for mathematical models of cancer therapies, _An application of geometric methods_ , Springer.
  * [10] Badri, H., Leder, K., 2016, Optimal treatment and stochastic modeling of heterogeneous tumors, _Biol. Direct_ , Vol. 11(1), n.a. 40.
  * [11] Sowndarrajan, P.T., Manimaran, J., Debbouche, A., Shangerganesh, L., 2019, Distributed optimal control of a tumor growth treatment model with cross-diffusion effect, _Eur. Phys. J. Plus_ , Vol. 134, n.a. 463
  * [12] Lions, J.-L., 1971, _Optimal Control of Systems Governed by Partial Differential Equations_. Springer: Berlin, Heidelberg.
  * [13] Fursikov, A.V., 2000, _Optimal Control of Distributed Systems: Theory and Applications_. Am. Math. Soc.: Providence.
  * [14] Kolobov, A.V., Gubernov, V.V., Polezhaev, A.A., 2007, Fisher type waves in a model of invasive tumor growth, _Math. Model._ , Vol. 19(6), pp. 31–42.
  * [15] Hecht, F., 2012, New development in FreeFem++, _J. Numer. Math._ , Vol. 20(3-4), pp. 251–265.