Non-Local Pearson diffusions

In this paper we focus on strong solutions of some heat-like problems with a non-local derivative in time induced by a Bernstein function and an elliptic operator given by the generator or the Fokker-Planck operator of a Pearson diffusion. Such kind of non-local equations naturally arise in the treatment of particle motion in heterogeneous media. In particular, we use spectral decomposition results for the usual Pearson diffusion to exploit explicit solutions of the aforementioned equations. Moreover, we provide stochastic representation of such solutions in terms of time-changed Pearson diffusions. Finally, we exploit some further properties of these processes, such as limit distributions and long/short-range dependence.

1. Introduction. Pearson distributions [55] constitute a family of probability distributions whose density functions m(x) satisfy the so called Pearson equation as x ∈ E ⊆ R. As one defines the polynomial µ(x) = a 0 + a 1 x and D(x) = d 0 + d 1 x + d 2 x 2 and sets σ(x) = 2D(x), Pearson equation is satisfied by the stationary measure of the solutions of the Stochastic Differential Equation (SDE) dX(t) = µ(X(t))dt + σ(X(t))dW (t), where W (t) is a standard Brownian motion. For this reason, such kind of diffusions, called Pearson diffusions [26], are statistically tractable with many different tools. For instance, they naturally arise while studying nonlinear time series linked to dynamical models with random perturbations (see, for instance, [52]). An important tool to deal with Pearson diffusions is the spectral decomposition, proposed in [44]. Indeed Pearson diffusions can be subdivided in three spectral categories depending on the spectrum of the generator. In the first spectral category, we have Pearson diffusions whose generator admit purely discrete spectrum and then a spectral decomposition of the transition density of the process follows easily. In the second spectral category, we have Pearson diffusions whose generator admits a discrete spectrum and an absolutely continuous spectrum (both simple) separated by a cut-off value. Such Pearson diffusions are more difficult to study due to the fact that the eigenfunctions in the continuous spectrum are expressed in terms of hypergeometric functions, as shown in [42] and [11]. Finally, in the third spectral category we only have Student diffusions, whose generator admits a simple discrete spectrum and an absolutely continuous spectrum of multiplicity two separated by a cut-off value. This case has been studied in [43,5]. Second and third spectral category go under the name of heavy tailed Pearson diffusions and their spectral properties have been analysed in [10]. In the modern theory of diffusion processes, anomalous diffusions have shown to be quite useful (see [50,29]) and stochastic models for such diffusions have been widely studied (see [60,28,47]). Anomalous diffusions models naturally arise when studying the motion of particles in heterogeneous media. For such kind of models, different analytic techniques are known (see, for instance, [18]). Such models fall into the class of fractional motions, which are described via a Langevin-type equation (see [25]). A class of anomalous diffusions can be seen as limit of continuous time random walks whose inter-jump times have infinite mean. Such limits lead to the appearance of the fractional Caputo derivative in the Fokker-Planck equation (see [50]). Despite Laplace transforms methods let us obtain the solution up to inverse Laplace transform, fractional Fokker-Planck equations are not usually solved explicitly. In this context, the spectral decomposition of Pearson diffusions has been revealed to be a quite powerful tool to explicitly express strong solutions of fractional diffusion equations. This lead to the introduction of fractional Pearson diffusions: in [40] the first spectral category was covered, in [41] the authors focused on the second spectral category, while, up to our knowledge, there were no known methods to extend the spectral decomposition to the fractional case in the third spectral category. In the second case, to express strong solutions of fractional Kolmogorov equations, the semigroup approach presented in [13] has been used. In both cases, the stochastic representation of such solutions is given by means of time-changed Pearson diffusions (with inverse stable subordinators). Similar strategies have been shown to work for lattice approximation of fractional Pearson diffusions: in [7] the spectral decomposition of a fractional immigration-death process (that is the lattice approximation of the Ornstein-Uhlenbeck process) is presented. On the other hand, subordinated Pearson diffusions (in particular the subordinated Jacobi process) have been shown to be useful tools to obtain large deviation principles (see [23]). Let us also remark that spectral properties of time-changed Markov processes can be used to determine the behaviour of their correlation structure, as shown in [39,53]. However, the classical fractional model is not the unique way one can achieve an anomalous diffusion. Indeed, anomalous diffusions can be introduced also by considering more general relaxation patterns than the exponential or the Mittag-Leffler ones (see [51]). Different inter-jump times in continuous time random walks and different relaxation patterns lead to a wider class of anomalous diffusions obtained via subordination of a random process (see [65]). In this case, one has to consider a more general non-local derivative in place of the Caputo one in Fokker-Planck equations. In [34] and [67] other non-local derivatives have been constructed, with the aid of Bernstein functions [61], that are the Laplace exponents of subordinators. A non-local derivative that can be defined in this way is, for instance, the tempered fractional derivative, which is an alternative to the classical fractional derivative that preserve the finiteness of the mean of the associated subordinator. The tempered model has been widely used to describe anomalous diffusions in heterogeneous media (in particular in geophysics, as in [45,71]) and its properties are well-known (see, for instance, [2]). A wide discussion on relaxation patterns can be found in [49]. Generalized fractional calculus is strictly linked with the definition of time-changed Markov processes. Indeed, for instance, in [67,20] a link between abstract generalized fractional differential equations and time-changed semigroups is established, while in [22] properties of the Green measures of time-changed Markov properties are exploited. Moreover, such non-local derivatives have been used for instance to define a class of non-local birth-death processes (see [8]) whose stationary distributions fall into the Katz family, which are discrete analogous of Pearson diffusions of the first spectral category. On the other hand, a first generalization to the general non-local setting of a Pearson diffusion (in particular the Ornstein-Uhlenbeck process) has been achieved in [27]. Let us also recall that in [23] the authors study a subordinated Jacobi process (with an inverse Gaussian subordinator), thus leading to a time-changed version of the aforementioned process.
The aim of the paper is to give a unified framework for working with non-local generalizations of Pearson diffusions and non-local equations that are linked to them, exhibiting exact strong solutions for this particular class of non-local advection-diffusion equations. In particular, we focus on non-local Cauchy problems of the form x ∈ E where G and F are respectively the generator and the Fokker-Planck operator of a Pearson diffusion and ∂ Φ t is a Caputo-type non-local derivative linked to a Bernstein function Φ. We use the notion of time-changed process to introduce the family of non-local Pearson diffusions and then provide the spectral decomposition of the transition densities of such diffusions. Moreover, we use both the spectral decomposition and a semigroup approach to exploit strong solutions of the aforementioned non-local Cauchy problems. Let us remark that in this paper we cover also the case of the third spectral category, for which spectral results for non-local generalizations were (up to our knowledge) unknown even in the standard fractional case. Section 2 is dedicated to a complete recap on spectral properties of classical Pearson diffusions, while Section 3 presents some basic properties of Bernstein functions, inverse subordinators and non-local derivatives. We also specify the semigroup theory results given in [13,67,20] to the case of complete Bernstein functions, in which better regularity can be proven. The proof of such result is quite technical, thus it is left in Appendix A. The definition and some preliminary properties of non-local Pearson diffusions are presented in Section 4, with particular attention to the existence of a transition probability density. In Section 5 we focus on non-local Pearson diffusions of the first spectral category, extending the results given in [40] to general Bernstein functions.
In Section 6 we consider non-local Pearson diffusions of the second (adapting the approach given in [41]) and the third spectral categories. Section 7 is devoted to stochastic representation results of strong solutions of (1.1) and (1.2). Finally, in Section 8 we explore some further properties of the non-local Pearson diffusions, as limit distributions, first-order stationarity and long/short-range dependence.
2. Pearson diffusions and their spectral classification. From now on let us fix a filtered probability space (Ω, Σ, F t , P). Let us give the definition of Pearson diffusion as presented in [26].
is a diffusion process satisfying the following stochastic differential equation where W (t) is a standard Brownian motion and µ(x) and σ 2 (x) are polynomials respectively of at most first and second degree.
In particular let us set In particular the diffusion space of X(t) is given by an interval E = (l, L) in which the polynomial D(x) is positive (hence the function σ(x) is real and non-zero). Moreover, let us recall that the family of operators (T (t)) t≥0 acting on C 0 (E), i.e. the space of continuous functions on the closureĒ of E that are 0 at infinity (if E is bounded then C 0 (E) = C(Ē)), given by where P x (·) = P(·|X 0 = x), is a uniformly bounded strongly continuous C 0 -semigroup on L 2 (E) (this is proved, for instance, in [40,41]). In particular its generator G is defined as with operator core C 2 b (E). Concerning the process X(t), the transition probability density p(t, x; y) is well defined for x, y ∈ E. In particular, it is solution of the following Cauchy problem (the backward problem) where δ x is a Dirac delta centred in x ∈ E.
On the other hand, we can also define the Fokker-Planck operator, or just forward operator, as with operator core C 2 (E). In particular p(t, x; y) is also solution of the Cauchy problem (the forward problem): In particular, if the process X(t) admits a stationary measure m, then its density m(x) must satisfy a stationary version of the forward equation, that becomes: Such equation is called Pearson equation, by the fact that it was introduced in [55] to classify some important classes of distributions. From now on we will only work with Pearson diffusions that admit a stationary measure that is also the limit measure of the process. In particular this means that, up to a re-parametrization, After such observation, we can recognize six different Pearson diffusions depending on the coefficients of D: • If D ≡ d 0 , then X(t) is a Ornstein-Uhlenbeck (OU) process and the stationary distribution is a Gaussian distribution; is a Cox-Ingersoll-Ross (CIR) process and the stationary distribution is a Gamma distribution; is a Jacobi process and the stationary distribution is a Beta distribution; (FS) process and the stationary distribution is a Fisher-Snedecor distribution; (RG) process and the stationary distribution is a reciprocal Gamma distribution; is a Student process and the stationary distribution is a Student distribution. Depending on the property of the spectrum of the generator G, in [41] these distributions were subdivided in three categories depending on the spectral category of Linetski classification (see [44,Theorem 3.2]): • The first spectral category contains the OU, the CIR and the Jacobi processes: their generator G admits purely discrete spectrum with infinitely many single non-positive eigenvalues (λ n ) n∈N ; • The second spectral category contains the FS and the RG processes: their generator G admits a discrete part and an absolutely continuous part that are disjoint and both of multiplicity one; • The third spectral category contains only the Student processes: its generator G admits a discrete part of multiplicity one and a disjoint absolutely continuous part of multiplicity two. Such classification is based on the oscillatory/non-oscillatory behaviour of the endpoints of the Sturm-Liouville equations G f = −λf (see [24,69,3]). Let us give some details on each Pearson diffusion. In particular we will re-parametrize again the polynomials D(x) and µ(x) in a form that will make the writing of the parameters of the stationary distributions easier. Moreover, we define the normalized polynomials, where the normalization constant is chosen with respect to the stationary density m(x), which is also the orthogonality density of the polynomials. In the following we will denote by Q n the normalized polynomials and with K n the normalization constants.
as θ > 0 and µ, σ ∈ R. The diffusion space is given by E = R and its stationary density is a Gaussian one, given by Concerning the eigenvalue equation G f = −λf , it admits solutions for λ n = θn as n ≥ 0 and its solutions are the Hermite polynomials (see [62]), defined by means of the Rodrigues formula (see [31]) with normalization constants K n = σ n √ n! . 2.1.2. The CIR process. The CIR process is solution of the SDE where θ, a, b > 0. The diffusion space is given by E = (0, +∞) and its stationary density is the Gamma one, given by n (x) are defined by the Rodrigues formula (see [31]) with normalization constants K n = Γ(b)n! Γ(b+n) . 2.1.3. The Jacobi process. The Jacobi process is solution of the SDE The diffusion space is E = (−1, 1) and its stationary density is a Beta one The eigenvalue equation G f = −λf admits solutions for λ n = nθ(n+a+b+1) a+b+2 and the eigenfunctions are Jacobi polynomials defined by the Rodrigues formula (see [31]) with normalization constants K n = 2 a+b+1 Γ(n+a+1)Γ(n+b+1) (2n+a+b+1)Γ(n+a+b+1)n! . Let us recall that some properties of statistical interest for the Jacobi process have been highlighted in [23].
as α, θ > 0 and β > 2. The diffusion space is E = (0, +∞) and its stationary density is a Fisher-Snedecor one: A spectral analysis of the FS process has been carried in [11]. In particular it is not difficult to see that the Fisher-Snedecor density admits finite even moments up to 2N 1 as N 1 = ⌊ β 4 ⌋. Thus we have a finite number of simple eigenvalues in the discrete spectrum of − G. The eigenvalues are given by The absolutely continuous spectrum σ ac (− G) = (Λ, +∞) where the cut-off Λ is given by .
For λ ∈ (Λ 1 , +∞), the fundamental solutions of the Sturm-Liouville equation G f = −λf are given in terms of hypergeometric functions. A detailed study is made in [11]. In any case, the solution that appears in the absolutely continuous part of the spectral decomposition is Let us recall that the solution f 1 is found by making use of the theory of Kummer's solutions for Hypergeometric Equations (see [64,56]).

2.2.2.
The RG process. The RG process is solution of the SDE as α, θ > 0 and β > 1. The diffusion space is E = (0, +∞) and its stationary density is a reciprocal Gamma one: As for the Fisher-Snedecor distribtuon, let us recall that such density admits finite even moments up to 2N 2 as N 2 = β 2 . A complete spectral analysis of the RG process has been made in [42]. We have a finite number of simple eigenvalues in the discrete spectrum of − G, given by The eigenfunctions are Bessel polynomials B (α,β) n defined by the Rodrigues formula with normalizing constant The absolutely continuous spectrum σ ac (− G) = (Λ, +∞) where the cut-off Λ is given by .

Spectral decomposition theorem for spectral category II.
Let us recall the spectral decomposition theorem for Pearson diffusions of spectral category II as given in [11,42].
Theorem 2.2. Let X(t) be a Pearson diffusion of spectral category II and, if X(t) is a FS process, let α > 2 with α = 2(m + 1) for any m ∈ N. The density p(t, x; x 0 ) admits the following spectral decomposition: where j = 1, 2 and

Pearson diffusions of spectral category III: Student diffusions. The Student process is solution of the SDE
where θ, δ > 0, µ, µ ′ ∈ R and ν > 1. If µ = µ ′ we refer to it as the symmetric Student process, while for µ = µ ′ it is called skew Student process. In [10] it is shown that such process admits diffusion space E = R and stationary density The Student density admits finite even moments up to 2N 3 as N 3 = ν 2 . Thus the discrete spectrum of (− G) is made of a finite number of simple eigenvalues given by while the eigenfunctions are the generalized Routh-Romanovski polynomials R n (x) defined by the Rodrigues formula where c(ν, δ, µ, µ ′ ) is a suitable constant. If µ = µ ′ we obtain the classical Routh-Romanovski polynomials R n (x) and their normalizing constant (up to a multiplicative constant). Concerning the absolutely continuous spectrum σ ac (− G) = (Λ 3 , +∞), where the cut-off Λ 3 is given by we have to recall that its elements are of multiplicity 2.
Concerning the eigenfunctions, in place of the monotonic solutions considered in [10], we can consider two independent linear combination of them. Indeed, by using the same strategy adopted in [43], referring to Kummer solutions (see [56,64]), we can consider one of the solutions of G f = −λf and its complex conjugate.
In [43] the eigenfunctions are explicitly expressed in the case µ = µ ′ , while for µ = µ ′ , the eigenfunctions are studied in [5], by reconsidering the seminal paper [70]. For the sake of shortness, we omit the explicit formula for such eigenfunctions, but we refer to it as f 3 (x; −λ) and f 3 (x; −λ). Using Linetski approach (see [44]) as done in [10], we obtain the following spectral decomposition theorem.
Theorem 2.3. The Student process is ergodic for ν = 2k − 1 for any k ∈ N. Under such choice of parameters, the density p(t, x; x 0 ) admits the following spectral decomposition: Let us also underline the expression of the speed density that will be useful in what follows stressing that +∞ −∞ sp(x)dx = M < +∞. 3. Inverse subordinators and non-local convolution derivatives. Now let us introduce our main object of study. Let us denote by BF the convex cone of Bernstein functions, that is to say Φ ∈ BF if and only if Φ ∈ C ∞ (R + ), Φ(λ) ≥ 0 and for any n ∈ N In particular it is known that for Φ ∈ BF the following Lévy-Khintchine representation ( [61]) is given is called the Lévy triple of Φ. Also the vice versa can be shown, i.e. for any Lévy triple (a Φ , b Φ , ν Φ ) such that ν Φ is a Lévy measure satisfying the integral condition (3.2) there exists a unique Bernstein function Φ such that Equation (3.1) holds. In the following we will consider a Φ = b Φ = 0 and ν Φ (0, +∞) = +∞. Moreover, let us denote byν Φ (t) = ν Φ (t, +∞). It is also known (see [61]) that for each Bernstein function Φ ∈ BF there exists a unique subordinator σ Φ = {σ Φ (y), y ≥ 0} (i. e. an increasing Lévy process) such that For general notions on subordinators we refer to [15,Chapter 3] and [16]. Our hypothesis on a Φ , b Φ and ν Φ imply that σ Φ is a pure jump strictly increasing process. For any Bernstein function Φ we can define the inverse subordinator L Φ as, for any t > 0 L Φ (t) is absolutely continuous for any t > 0 and we can denote by f Φ (s; t) its density. Let us recall (see [46] Moreover, let us observe that, under our assumptions, the sample paths of L Φ (t) are almost surely increasing and continuous. Now let us recall the definition of non-local convolution derivative, defined in [34] and [67].
Definition 3.1. Let u : R + → R be an absolutely continuous function. Then we define the non-local convolution derivative induced by Φ of u as Let us observe that one can define also the regularized version of the non-local convolution derivative as observing that it coincides with the previous definition on absolutely continuous functions. It can be shown, by Laplace transform arguments (see, for instance [35,6]) or by Green functions arguments (see [36]), that the Cauchy problem admits a unique solution for any λ ∈ R and it is given by e Φ (t; λ) := E[e λLΦ(t) ]. In [8] the following proposition is proved.
A Bernstein function Φ is said to be complete if its Lévy measure ν Φ (dt) admits a density ν Φ (t) that is completely monotone. Following the approach given in [13] we are able to prove the following theorem (the proof is given in Appendix).
Theorem 3.3. Let (X, · ) be a Banach space and (T (t)) t≥0 be a uniformly bounded and strongly continuous C 0 -semigroup on X. Define the family of linear operators on X (T Φ (t)) t≥0 as where Φ is a driftless complete Bernstein function, f Φ (s; t) is the density of the inverse subordinator L Φ (t) associated to Φ and the integral has to be intended in Bochner sense. Then (T Φ (t)) t≥0 is a uniformly bounded and strongly continuous family of linear operators. Moreover, it is also strongly analytic in a suitable sector is the generator of the semigroup (T (t)) t≥0 and u ∈ D(A), then also T Φ (t)u ∈ D(A) and it solves the Cauchy problem where the equality holds in X (and not necessarily pointwise).
Concerning, strong continuity, we refer to the definition given in [54] for general families of operators of one parameter. Moreover, let us observe that a similar result has been shown in [20, Theorem 2.1]. However, following the lines of [13], we are also able to show better regularity of the family of operators (T Φ (t)) t≥0 . Let us give some examples of complete Bernstein functions and associated subordinators.

Definition of the non-local Pearson diffusions.
Let us now define the non-local Pearson diffusions.
Definition 4.1. Let X(t) be a Pearson diffusion and Φ ∈ BF be a driftless Bernstein function. Let L Φ (t) be an inverse subordinator associated to Φ and independent of X(t). Then we define the non-local Pearson diffusion induced by X and L Φ as X Φ (t) := X(L Φ (t)).
The first property one has to recall is that X Φ (t) is not a Markov process, but it is still a semi-Markov one (see [21]). The fact that L Φ (t) is a.s. continuous implies that the Brownian motion W (t) in the stochastic differential equation (2.1) is in synchronization (in the sense of [33]) with L Φ (t) for any driftless Bernstein function Φ. Thus the first change-of-variable formula for synchronized processes (see [33,Lemma 2.3]) and the duality theorem [33,Theorem 4.2] lead us to the following characterization of non-local Pearson diffusions (following the lines of [41]).
Thus we can characterize the non-local Peason diffusions as the unique strong solutions of Equation (4.1), in analogy of what is done with the classical Pearson diffusion. In particular, from this result, we still have that non-local Pearson diffusions are semimartingales, but with respect to a time-changed filtration F LΦ(t) . Now let us show that non-local Pearson diffusions admit in some sense a transition probability density.
Moreover, the following integral representation holds: where p is the transition probability density of X(t).
Proof. Let us first observe that since L Φ (0) = 0 almost surely, then X Φ (0) = X(0) almost surely. Now fix B ∈ B(E) and observe, by the independence of L Φ and X and by definition of p(t, x; x 0 ), where we used Fubini's theorem since all the integrands are non-negative. Thus, by uniqueness of the Radon-Nikodym derivative of measures, we have Equation (4.2).
We call p Φ (t, x; x 0 ) the transition probability density of X Φ (t).
Let us stress out that, for any Pearson diffusion, if we consider the semigroup T (t)g(x) = E x [g(X(t))] acting on C 0 (E), as a direct consequence of the previous lemma, the family of operators (T Φ (t)) t≥0 acting on C 0 (E) and defined as T Φ (t)g(x) = E x [g(X Φ (t))] is obtained by time-changing T (t). Moreover, we can characterize the family of adjoint operator (T * Φ (t)) t≥0 of T Φ (t). Lemma 4.4. Let X Φ (t) be a non-local Pearson diffusion and T Φ (t) be the family of operators defined on C 0 (E) as T Φ (t)g(x) = E x [g(X Φ (t))]. Then the adjoint operators (T * Φ (t)) t≥0 are defined as x; y)f (y)dy and represent the time-changed semigroup of the strongly continuous semigroup (T * (t)) t≥0 defined as We omit the proof since it easily follows from the definition of T * Φ and the integral representation given in Lemma 4.3. In the following sections we will investigate the spectral decomposition of the transition probability density of non-local Pearson diffusions and the existence of strong solutions for the associated backward and forward Kolmogorov problems.
5. Spectral decomposition of non-local Pearson diffusions of spectral category I. 5.1. Spectral decomposition of the transition probability density. Let us consider a Pearson diffusion X(t) of spectral category I with diffusion space E, generator G, stationary density m(x) and orthonormal polynomials Q n (x). Then it is well known that the spectral decomposition of p(t, x; x 0 ) is given by where λ n are the eigenvalues of − G for each Q n . We want to show a similar spectral decomposition for p Φ (t, x; x 0 ) where, in place of e −λnt , we have e Φ (t; −λ n ). To do this, we first need to show the following technical Lemma.
Lemma 5.1. Let Φ ∈ BF be a driftless Bernstein function, X(t) be a Pearson diffusion of spectral category I with diffusion space E, associated family of classical orthonormal polynomials Q n (x) and λ n the respective eigenvalues. Then the series +∞ n=0 e Φ (t; −λ n )Q n (x)Q n (y) absolutely converges for fixed t > 0 and x, y ∈ E and totally converges for x, y belonging to compact sets in E and t ≥ t 0 .
Proof. Let us recall that, by Equation (3.5), for fixed t > 0 there exists K > 0 such that e Φ (t; −λ n ) ≤ K λn . Let us first work with the OU process. We can suppose without loss of generality that µ = 0 and σ = 1. Then there exists a constant C > 0 such that where the estimate follows from the definition of the normalizing constant and the estimate on Hermite polynomials given in [59,Page 369]. Moreover, λ n = θn, thus obtaining the absolute convergence. Total convergence follows from the same estimates and the fact that e Φ (t; −λ n ) ≤ e Φ (t 0 ; −λ n ).
Concerning the CIR process, we can suppose without loss of generality a = 1. We have, considering the normalizing constant and the estimate on Laguerre polynomials given in [59,Page 348], that there exists a constant C independent of x and n such that, for n big enough, it holds thus we have again (5.2) and both absolute and total converge follow. Now let us consider the Jacobi process case. By [66,Theorem 8.21.8] and the definition of the normalizing constant we have where the remainder term and the function C(x, a, b) are uniform in x as x ∈ K x ⊂ (−1, 1) with K x compact set and cos(θ) = x. Thus, in this case, we have the bound being uniform as t ≥ t 0 > 0, x ∈ K x and y ∈ K y where K x and K y are compact sets in E = (−1, 1). Now we are ready to give the spectral decomposition of the transition probability density p Φ (t, x; x 0 ).
Theorem 5.2. Let Φ ∈ BF be a driftless Bernstein function, X(t) be a Pearson diffusion of spectral category I with diffusion space E, associated family of classical orthonormal polynomials Q n (x) with eigenvalues λ n . Let X Φ (t) be the respective non-local Pearson diffusion. Then for any t > 0 and x, x 0 ∈ E.
Proof. By using Equation (4.2) we know that Now we have to show that we can exchange the order of the integral with the summation. To do this, let us observe that +∞ 0 e −λns |Q n (x)Q n (x 0 )|f Φ (s; t)ds = e Φ (t; −λ n )|Q n (x)Q n (x 0 )| < +∞ by the previous Lemma and then we can use Fubini's theorem to conclude the proof. Now that we have the spectral decomposition of the transition density, we can focus on showing the existence of strong solutions to the non-local backward Kolmogorov equation under suitable assumptions on the initial datum.

The non-local backward Kolmogorov equation.
Here we want to focus on the following non-local Cauchy problem y ∈ E for suitable initial data g. • t ≥ 0 → u(t, ·) ∈ L 2 (m(x)dx) is strongly continuous; • The function ∂ Φ t u(t, y) is well-defined for any t > 0 and y ∈ E; • The convolution integral involved in ∂ Φ t u(t, y) is a Bochner integral and the derivative is a strong derivative in L 2 (m(x)dx); • The function t > 0 → ∂ Φ t u(t, ·) ∈ L 2 (m(x)dx) is strongly continuous; • For fixed t > 0, u(t, ·) ∈ C 2 (E); • The equations of (5.5) hold pointwise for any t > 0 and almost all y ∈ E.
Let us first show the next result, following the lines of [40].
Theorem 5.4. Let X(t) be a Pearson diffusion of spectral category I with diffusion space E, associated family of orthonormal polynomials Q n with respective eigenvalues λ n , generator G and stationary density m(x). Let Φ ∈ BF be a driftless Bernstein function. Let g ∈ L 2 (m(x)dx) be decomposed as g(y) = +∞ n=0 g n Q n (y), where the series converges in L 2 (m(x)dx), absolutely for fixed y ∈ E and uniformly on compact intervals [y 1 , y 2 ] ⊂ E. Then the unique strong solution in L 2 (m(x)dx) of the problem (5.5) is given by Moreover, p Φ (t, x; y) is the fundamental solution of the problem (5.5), in the sense that, for any initial datum g ∈ L 2 (m(x)dx) such that g = +∞ n=0 g n Q n converges in L 2 (m(x)dx), absolutely for fixed y ∈ E and uniformly on compact interval [y 1 , y 2 ] ⊂ E, it holds u(t, y) = E p Φ (t, x; y)g(x)dx.
Proof. First of all, let us observe that the function in Equation (5.6) is well defined, by showing that the involved series converges in L 2 (m(x)dx). To do this, let us first define for N ∈ N u N (t, y) = N n=0 e Φ (t; −λ n )Q n (y)g n .
By Parseval's identity, setting N > M , we have that However, still using Parseval's identity, we have +∞ n=0 g 2 n = g 2 L 2 (m(x)dx) < +∞, thus concluding the convergence in L 2 (m(x)dx) of the sequence u N (t, ·) by Cauchy's criterion. Now let us show that t → u(t, ·) is strongly continuous in L 2 (m(x)dx). Concerning strong continuity at 0, let us fix ε > 0 and consider N > 0 such that +∞ n=N +1 g 2 n < ε. Thus we have, by Parseval's identity, Taking the limit superior as t → 0 + we have Since ε > 0 is arbitrary, we have continuity at 0 + . Strong continuity at any t > 0 is proven in an analogous way. Now let us observe that for the single summand it holds (5.7) ∂ Φ t e Φ (t; −λ n )Q n (y)g n = −λ n e Φ (t; −λ n )Q n (y)g n = G e Φ (t; −λ n )Q n (y)g n for any n ∈ N. Thus, we actually have to show that we can exchange the operators G and ∂ Φ t with the summation. Let us first consider ∂ Φ t +∞ n=0 e Φ (t; −λ n )Q n (y)g n . First of all, let us denote by I Φ (t) = t 0ν Φ (s)ds the integrated tail of the Lévy measure. I Φ (t) is an increasing non-negative function with derivativeν Φ (t), hence we can rewrite t 0ν Φ (t − τ )(u(τ, y) − g(y))dτ = where the right-hand side converges for fixed y ∈ E. Thus the series in the right-hand side of (5.8) totally converges and we have, by using [58, Theorem 7.16], Now we only need to show that we can exchange the derivative d dt with the summation sign on the right-hand side of Equation (5.9). To do this, we have to show that the series of the derivatives converges uniformly in any compact set containing t (see [58,Theorem 7.17]). However, the series of the derivatives is given by +∞ n=0 −λ n e Φ (t; −λ n )g n Q n (y). Now fix t 0 > 0 and t ≥ t 0 . Then we have, by Equation (3.5) and the fact that t → e Φ (t; −λ) is decreasing for any λ > 0, λ n e Φ (t; −λ n ) ≤ M (t 0 ). In particular, this implies +∞ n=0 λ n e Φ (t; −λ n )|g n Q n (y)| ≤ M (t 0 ) +∞ n=0 |g n Q n (y)|, and then the series of the derivative is totally convergent in any compact set separated from 0. Finally, we have −λ n e Φ (t; −λ n )g n Q n (y).
Arguing as we did for the series on the right-hand side of Equation (5.6), one can easily show that the series on the right-hand side of Equation (5.10) strongly converges in L 2 (m(x)dx) and that ∂ Φ t u(t, y) is strongly continuous in L 2 (m(x)dx) as t > 0. Now we have to show that we can exchange the operator G with the series on the right-hand side of Equation (5.6). Since G is a second order differential operator with polynomial coefficients, we only need to exchange the first and the second derivative with respect to y with the series. To do this, we need to argue differently depending on the process. Let us first consider the OU process. The differential recurrence relation between Hermite polynomials (see, for instance, [1, Formula 22.8.8]) becomes, after normalization Q ′ n (y) = √ nQ n−1 (y). Thus, combining this relation with (5.1) we have that, in any compact set [y 1 , y 2 ] ⊂ R containing y, it holds where the right-hand side is the summand of a convergent series, since both (n −3/4 ) n∈N and (g n ) n∈N belong to ℓ 2 . Concerning the second derivative, the Sturm-Liouville equation G Q n (y) = −λ n Q n (y) can be rewritten as Q ′′ n (y) = yQ ′ n (y) − nQ n (y), thus we have to study only the term in nQ n (y). For such term, recalling that λ n = θn, we know that n e Φ (t; −λ n ) ≤ M (t) for some M (t), hence the uniform convergence of +∞ n=0 Q n (y)g n on compact intervals implies the uniform convergence of +∞ n=0 e Φ (t; −λ n )nQ n (y)g n on compact intervals. Now let us consider the CIR process.
After normalization, [66, Formula 5.1.14] becomes where we denote with the indexes the dependence on the parameter b. By using (5.3), we obtain again equation (5.11) uniformly on compact intervals, thus obtaining the desired converges. Moreover, the Sturm-Liouville equation G Q n (y) = −λ n Q n (y) becomes yQ ′′ n (y) = (y − b)Q ′ n (y) − nQ n (y) thus, since λ n = θn, we can argue as in the OU process case to obtain theuniform convergence of +∞ n=0 e Φ (t; −λ n )nQ n (y)g n on compact intervals. Last case is the Jacobi process case. In this case, normalizing [1, Formula 22.8.1], we have Both terms lead to a series of the form +∞ n=0 n e Φ (t; −λ n )Q n (y)g n . However, this time we know that Q n (y) are bounded and λ n ∼ Cn 2 , thus we have that (by still using (3.5)) | e Φ (t; −λ n )nQ n (y)g n | ≤ C g n n where the right-hand side is the summand of a convergent series since (g n ) n∈N , (1/n) n∈N ∈ ℓ 2 . Concerning the second derivative, the Sturm-Liouville equation G Q n (y) = −λ n Q n (y) becomes (1 − y 2 )Q ′′ n (y) = −((b − a) − (a + b − 2)y)Q ′ n (y) − n(n + a + b + 1)Q n (y), thus, since n(n + a + b + 1) and λ n are both quadratic in n, the same argument as in the previous cases leads to uniform convergence on compact intervals of +∞ n=0 e Φ (t; −λ n )n(n + a + b + 1)Q n (y)g n . Concluding we have in general G u(t, y) = +∞ n=0 G e Φ (t; −λ n )Q n (y)g n and then also u(t, ·) ∈ C 2 (E). Equation (5.7) and the fact we can exchange the order of the operators with the series imply that the equations of (5.
To do this, let us use Jensen's inequality and Parseval's identity to achieve Now let us show that the derivative is actually a strong derivative in L 2 (m(x)dx), that is to say Let us argue for h > 0, as the case h < 0 is analogous. By triangular inequality and Equation (5.10) we have .

(5.13)
Let us argue with the second summand; by Bochner's theorem we have thus, taking the limit and recalling that, by strong continuity of u(t, ·) in L 2 (m(x)dx) the function t → u(τ, ·) − u(0+, ·) L 2 (m(x)dx) is continuous, we have (5.14) lim Now let us consider the first summand of the right-hand side of Equation (5.13). To handle this summand, let us first observe that Being | e Φ (t + h − τ ; −λ n ) − e Φ (t − τ ; −λ n )| ≤ 1 and +∞ n=1 |g n Q n (y)| convergent, we have that the series on the right-hand side of Equation (5.15) normally converges, thus we achieve Moreover, the same argument as before shows that the series on the right-hand side of Equation (5.16) is absolutely convergent. Being also the series in Equation (5.10) absolutely convergent, we have .

(5.17)
Now let us show that the series on the right-hand side of Equation (5.17) converges in L 2 (m(x)dx). To do this, let us define and observe that, for any N > M where we used Equation (3.5). Being (g n ) n≥0 ∈ ℓ 2 , we have that S N converges in L 2 (m(x)dx) by Cauchy's criterion. This implies that we can use Parseval's identity in the right-hand side of (5.17) to achieve Now let us observe that J Φ belongs to C 1 and thus, by Lagrange's theorem, we know there exists θ n (h) ∈ [t, t + h] such that Using the latter in Equation (5.18) we obtain (5.19) Now let us show that the series in the right-hand side of (5.19) normally converges. Indeed, considering h ∈ [0, δ] for some δ > 0, we have Thus we can take the limit as h → 0 + inside the summation sign, obtaining, since e Φ (t; −λ n ) is continuous, that is to say λ n e Φ (t; −λ n )g n Q n L 2 (m(x)dx) = 0.
Uniqueness of the solution follows from the fact that (Q n ) n∈N is a complete orthonormal system in L 2 (m(x)dx). Finally, the fact that p Φ (t, x; y) is the fundamental solution of (5.5) follows from the spectral decomposition Theorem 5.2.
Remark 5.5. Let us observe that if g ∈ C 2 0 (E) and Φ is a complete Bernstein function, we can use Theorem 3.3 to obtain that u(t, y) defined in Equation (7.1) is solution of (5.5) and can be extended to an analytic function in a sector C(α) for some α < π 2 .

The non-local forward Kolmogorov equation.
Here we want to focus on the following non-local Cauchy problem for suitable initial datum f . For the definition of strong solution, we still refer to Definition 5.3. To do this, we first need to show the following preliminary result, whose proof can be omitted since it follows by direct calculations and by using Pearson equation (2.4).
Lemma 5.6. Let X(t) be a Pearson diffusion of spectral category I with diffusion space E, associated family of orthonormal polynomials Q n with respective eigenvalues λ n , stationary density m(x) and Fokker-Planck operator F . Then F (m(x)Q n (x)) = −λ n m(x)Q n (x). Now that we have proven this Lemma, we can actually show existence, uniqueness and spectral decomposition of the strong solutions (in L 2 (m(x)dx)) of (5.21).
Theorem 5.7. Let X(t) be a Pearson diffusion of spectral category I with diffusion space E, associated family of orthonormal polynomials Q n with respective eigenvalues λ n , Fokker-Planck operator F and bounded stationary density m(x). Let Φ ∈ BF be a driftless Bernstein function. Let f : E → R be such that f /m ∈ L 2 (m(x)dx) is decomposed as where the series converges in L 2 (m(x)dx), absolutely for fixed x ∈ E and uniformly on compact intervals [x 1 , x 2 ] ⊂ E. Then the unique strong solution in L 2 (m(x)dx) of the problem (5.5) is given by Moreover, p Φ (t, x; y) is the fundamental solution of the problem (5.5), in the sense that, for any initial datum f such that f /m ∈ L 2 (m(x)dx) and f /m = +∞ n=0 f n Q n converges in L 2 (m(x)dx), absolutely for fixed x ∈ E and uniformly on compact interval [x 1 , Proof. Since we are supposing that m(x) is bounded, it is easy to see that if a function f is such that f /m ∈ L 2 (m(x)dx), then f ∈ L 2 (m(x)dx). Thus the convergence of the series in (5.22) in L 2 (m(x)dx) is ensured by the convergence of the series defined by v(t, x)/m(x). Moreover, the strong continuity of v(t, x) follows from the one of v(t, x)/m(x). For the single summand of (5.22), it holds, by using the previous Lemma, Finally, the exact same proof of Theorem 5.4 still works in this case.
Remark 5.8. If m is not bounded (i.e. in the cases of the CIR process as b ∈ (0, 1) and of the Jacobi process as a, b ∈ (−1, 0)), one can prove in the same way that v(t, x) is a strong solution in L 2 (m −1 (x)dx). If, additionally, f ∈ L 2 (m(x)dx) and the same convergence assumptions on f /m hold also for f = +∞ n=0 Q n f n , then v(t, x) is a strong solution in L 2 (m(x)dx) even if m(x) is unbounded. Finally, let us recall that if f ∈ C 2 0 (E) and Φ is a complete Bernstein function, then, by Lemma 4.4 and Theorem 3.3, we know that v(t, ·) can be analytically extended in a sectorial region C(α) where α < π 2 . 6. Spectral decomposition of non-local Pearson diffusions of spectral category II and III. Now we consider the case of Pearson diffusions of spectral category II and III. We discuss them together since to exploit strong solutions of Equations (5.5) and (5.21) we have to directly use semigroup theory, being the spectral decomposition of p Φ (t, x; y) more complicated. Let us first focus on the spectral decomposition of the transition densities.
6.1. Spectral decomposition of the transition density for spectral category II. Now let us show a spectral decomposition theorem for the transition density p Φ (t, x; y) of a non-local Pearson diffusion of spectral category II.
Theorem 6.1. Let Φ ∈ BF be a driftless Bernstein function, X(t) be a Pearson diffusion of spectral category II with diffusion space E, associated family of classical orthonormal polynomials Q n (x) for n ≤ N j , with respective eigenvalues in the discrete spectrum of − G given by λ n and m(x) stationary density. Then, for any t > 0 and x, x 0 ∈ E, it holds where j = 1, 2 and the involved functions are defined in subsection 2.2.
Proof. By Theorem 2.2 we already know that and Moreover, by Equation (4.2) we know that Concerning p d,Φ , we already know that it is given by equation ( Moreover, we have that and |∆ j (λ)| ≥ C j λ 1 2 as λ → +∞. Thus we finally get concluding the proof.
6.2. Spectral decomposition of the transition density for spectral category III. Concerning the spectral decomposition of p Φ (t, x; y) for spectral category III, we will not use asymptotics of the eigenfunctions. Instead, we will use the fact that the continuous part of the spectrum is of multiplicity two on our advantage. Theorem 6.2. Let Φ ∈ BF be a driftless Bernstein function, X(t) be a Student process with diffusion space E = R, associated family of classical orthonormal polynomials Q n (x) for n ≤ N 2 with respective eigenvalues in the discrete spectrum of − G given by λ n and m(x) Student density. Then, for any t > 0 and x, x 0 ∈ E, it holds Proof. As in Theorem 6.1, we only need to deal with p c,Φ . Let us denote .
Then it is easy to see, by simple complex analysis arguments, that Now set x = x 0 and let us first show the spectral decomposition in this case. Observe that by Fubini-Tonelli theorem, since the integrands are all non-negative. Thus we obtain, by Equation (4.2), Now let us fix x ∈ R and consider any scale function s : R → R, i.e. a continuous strictly increasing function on R such that for any a, b ∈ R with a < b and y ∈ (a, b) it holds where T b = inf{t > 0 : X(t) > b} and T a = inf{t > 0 : X(t) < b} (see [57, Proposition 3.2, Chapter V II]). Let us denote by s −1 its inverse function. Thus, define where sp(y) is the speed density of X(t). V x (t) is strictly increasing and continuous, thus so it is its inverse function V −1 x (t). Moreover, as V x (0) = 0 (by dominated convergence theorem, since the speed measure is finite) we also have V −1 x (0) = 0 and then V −1 x (t) is bounded in a neighbourhood of 0. Then, by [68, Theorem 2.1], we know that there exist two constants t 0 , C > 0 (depending on x) such that for any t < t 0 it holds In particular this means that there exists ε > 0 (for instance ε = t 0 /2) such that ε 0 p(t, x; x)dt < +∞. Now let us recall that f Φ (s; t) is continuous for fixed t and f Φ (0 + ; t) =ν Φ (t). Thus, for s ∈ [0, ε], there exists a constant M (t) > 0 such that f Φ (s; t) ≤ M (t). Moreover, from the spectral decomposition, we have that p(t, x; x) is decreasing in t, thus p(t, x; x) ≤ p(ε, x; x) for any t ≥ ε. Thus we have Thus, by Equation (6.4) and the fact that both the summands are non-negative, we have that +∞ Λ3 e Φ (t; −λ)u 2 3 (x; −λ)dλ < +∞.
Now we can consider the case x 0 = x. Indeed, in such case, we have thus we can use Fubini's theorem also in the case x = x 0 to conclude the proof. Theorem 6.3. Let X(t) be a Pearson diffusion of spectral category II or III with diffusion space E and generator G. Let Φ ∈ BF be a driftless Bernstein function and g ∈ C 2 0 (E). Then the strong solution of (5.5) with initial datum g is given by Theorem 6.4. Let X(t) be a Pearson diffusion of spectral category II with diffusion space E and Fokker-Planck operator F . Let Φ ∈ BF be a driftless Bernstein function. Let f ∈ C 2 c (E). Then the strong solution of (5.21) with initial datum f is given by x; y)f (y)dy.
Remark 6.5. For both previous Theorems, if Φ is a complete Bernstein function, then u(t, x) and v(t, x) can be extended by analyticity to a sectorial region C(α) with α < π 2 . 7. Stochastic representation of the solutions of the backward and forward Kolmogorov equations. Let us recall that in the previous two Sections we have characterized p Φ (t, x; y) as the fundamental solution of both problems (5.5) and (5.21). Recalling that p Φ (t, x; y), we obtain some stochastic representation results. Concerning the backward Kolmogorov equation (5.5) we have the following result. The fact that p Φ (t, x; y) is the fundamental solutions leads to a stochastic representation result for strong solutions of (5.5).
Proposition 7.1. Let X(t) be a Pearson diffusion with diffusion space E, associated family of orthonormal polynomials Q n with respective eigenvalues λ n , generator G and stationary density m(x). Let Φ ∈ BF be a driftless Bernstein function and g ∈ C 2 0 (E). Then the unique strong solution of (5.5) is given by Proof. By definition of conditional expectation and the fact that p Φ (t, x; y) is the transition density of X Φ (t), we have u(t, y) = E p Φ (t, x; y)g(x)dx concluding the proof by observing that p Φ (t, x; y) is the fundamental solution of (5.5).
Remark 7.2. If X Φ (t) belongs to spectral category I, then one can weaken the hypotheses on the initial datum, by asking the same hypotheses as in Theorem 5.4 in place of g ∈ C 2 0 (E). However, since we are using adjoint operators, we have to ask for more strict hypotheses for the stochastic representation of strong solutions of (5.21). Proposition 7.3. Let X(t) be a Pearson diffusion with diffusion space E, associated family of orthonormal polynomials Q n with respective eigenvalues λ n , Fokker-Planck operator F and stationary density m(x). Let Φ ∈ BF be a driftless Bernstein function and f ∈ L 1 (dx) ∩ C 2 c (E) with f ≥ 0 and f L 1 (dx) = 1. Then the unique strong solution of the problem (5.21) is given by where P f (·) is the probability measure obtained by P conditioning with the request that X Φ (0) admits f (x)dx as probability law.
Proof. Let us observe that (setting p Φ (t, x; y) as p Φ (t, x; y)χ E (x)) x; y)f (y)dy concluding the proof by the fact that p Φ (t, x; y) is the fundamental solution of (5.21).
Remark 7.4. As for Proposition 7.1, if X Φ (t) belongs to spectral category I, then one can weaken the hypotheses on the initial datum, by asking the same hypotheses as in Theorem 5.7 in place of f ∈ C 2 c (E). 8. Limit distributions, first order stationarity and correlation structure of non-local Pearson diffusions. Now let us use the previous results to establish some properties of non-local Pearson diffusions X Φ (t). Let us first show that for regular enough initial distributions the limit distribution of any non-local Pearson diffusion admits is given by the stationary distribution m(x)dx of the respective Pearson diffusion.
Theorem 8.1. Let X(t) be a Pearson diffusion with diffusion space E, associated family of orthonormal polynomials Q n with respective eigenvalues λ n and stationary distribution m(x)dx. Let Φ ∈ BF be a driftless Bernstein function. Then for fixed x, y ∈ E.
Moreover, for any f ∈ C 2 for any fixed x ∈ E.
Proof. Let us first argue for X Φ (t) of spectral category I. First of all, we have (since for any Pearson diffusion λ 0 = 0 and Q 0 (x) ≡ 1): Since we are sending t → +∞, we can suppose t ≥ t 0 > 0 and then e Φ (t; −λ n ) ≤ e Φ (t 0 ; −λ n ). By Lemma 5.1, we can use dominated convergence theorem to obtain Equation (8.1). Concerning X Φ (t) of spectral category II, arguing as before, we have (for j = 1, 2) where 0 n=1 ≡ 0. Arguing as before and observing as in Theorem 6.1 that we can use again dominated convergence theorem to achieve Equation (8.1).
where u 3 is defined in the proof of Theorem 6.2. Since we have shown by Cauchy-Schwartz inequality that | e Φ (t 0 ; −λ)u j (x; −λ)u j (y; −λ)| is integrable in λ, we can still use dominated convergence theorem to achieve Equation (8.1). Finally, arguing as in Proposition 7.3, we have x; y)f (y)dy and then we conclude the proof by dominated convergence theorem.
Remark 8.2. If X Φ (t) belongs to spectral category I, then we can substitute the hypothesis f ∈ C 2 c (E) with the ones of Theorem 5.7.
In particular, we can show that m(x) is actually the first-order stationary distribution of X Φ (t).
Proposition 8.3. Let X(t) be a Pearson diffusion with diffusion space E, associated family of orthonormal polynomials Q n with respective eigenvalues λ n and stationary distribution m(x)dx. Let Φ ∈ BF be a driftless Bernstein function. Then it holds for any t ≥ 0 and x ∈ E.
Proof. Let us observe that where we used Tonelli-Fubini's theorem (since the integrand is non-negative) and the fact that m(x) is the stationary measure of X(t).
Remark 8.4. An interesting problem we leave open consists in giving estimates on the convergence rates of non-stationary non-local Pearson diffusions, as done in the classical case for the Fisher-Snedecor diffusion in [37].
However, the process X Φ (t) is not second-order stationary, neither in wide sense, as we can see from the following result.
Proposition 8.5. Let X(t) be a Pearson diffusion with diffusion space E, associated family of orthonormal polynomials Q n with respective eigenvalues λ n and stationary distribution m(x). Let Φ ∈ BF be a driftless Bernstein function. Then it holds, for any t ≥ s ≥ 0, Proof. This is a consequence of the fact that Corr m (X(t), X(s)) = e −λ1|t−s| and [9, Theorem 2].
Since, if X Φ (0) admits m(x) as probability density function, the process X Φ (t) is first-order stationary, but not second-order stationary, we cannot exploit long-range or short-range dependence properties with the usual definitions. Thus, taking inspiration from [14, Lemmas 2.1 and 2.2], we give the following definitions.
Proof of Theorem 3.3. Let us first show that the family of operators (T Φ (t)) t≥0 is well-defined and uniformly bounded. To do this, we need to use Bochner's theorem. Indeed we have, for any u ∈ X, Now let us show strong continuity of the family (T Φ (t)) t≥0 at 0 + . Indeed we have that t ∈ R + 0 → T (t)u ∈ X is a bounded continuous function (by strong continuity and uniform boundedness of the semi-group T (t)) and so it is t ∈ R + 0 → T (t)u ∈ R + 0 . Thus it holds Taking the limit we obtain, lim t→0 + T Φ (t)u − u ≤ T (0)u − u = 0, since, by dominated convergence theorem, f Φ (s; t) weakly converges towards δ 0 as t → 0 + . Thus we have strong continuity at 0 + . Strong continuity at any t > 0 can be analogously proven.
To show strong analyticity of (T Φ (t)) t≥0 in a certain sectorial region, we will argue by using Laplace transforms.