Gaussian process with Matérn kernel on a finite domain with periodic boundary conditions I'm concerned with a Gaussian process $f(x)$ on a finite domain $x\in[0,L)$ with periodic boundary conditions. Naively using the distance after accounting for periodic boundary conditions to evaluate the covariance kernel $K=k\left(\left\vert x-x'\right\vert\right)$ fails: $K$ is not positive (semi)-definite. I'd like to find a kernel $k$ that is "like" the Matern kernel in the sense that it reduces to the standard Matern kernel in the limit of the correlation length being much smaller than the domain. Any ideas on how to address this problem? I've managed to solve the problem for squared exponential kernels but am struggling with the Matern equivalent (see below).
We know that the squared exponential kernel is a solution to the heat equation, where the correlation length corresponds to the square root of time. We can solve the heat equation on the finite domain and impose periodic boundary conditions to obtain the "heat" kernel. It's not a pretty solution, but it's doable. This kernel is positive definite (energy can't be negative).
I was planning to use the same approach for the Matern kernel (find the PDE for which $k$ is the solution and solve it on the finite domain). Unfortunately, I cannot figure out the right PDE to solve. To be more concrete, what are the PDEs for which
$$
\begin{align}
k_{\nu=3/2}&=\left(1+\frac{\sqrt{3} x}{\ell}\right)\exp\left(-\frac{\sqrt{3} x}{\ell}\right)\\\\
\text{and }k_{\nu=5/2}&=\left(1+\frac{\sqrt{5} x}+\frac{5x^2}{3\ell^2}\right)\exp\left(-\frac{\sqrt{5} x}{\ell}\right)\\
\end{align}
$$
are solutions? I think what I'm after are second-order linear PDEs similar to the heat equation, although it appears there are some fractional sPDEs whose solution is the Matern kernel (but my skills are probably not sufficient to solve them on finite domains).
Thank you for your time!

Indeed it is a fractional SPDE, in particular $u$ is a Matern distributed if it solves:
$$
\tau(\kappa^2-\Delta)^\frac{\alpha}{2}u = \mathcal{W}
$$
where $\mathcal{W}$ is Gaussian white noise, $\alpha=\nu+\frac{d}{2}$ where ($d$ is dimension), $\Delta$ is the Laplacian operator, $\kappa^2$ is the lengthscale (so your $\ell$), and $\tau$ is a constant that control $y$'s variance with other parameters fixed (in particular $\sigma^2 = \frac{\Gamma(\nu)}{\Gamma(\alpha)(4\pi)^\frac{d}{2}\kappa^{2\nu}\tau^2}$).
To get more information about this process, I would look into the INLA approach to spatial modeling, which relies on speedy PDE evaluations to solve low-dimensional GPs. I got this information in particular from this review article: https://arxiv.org/pdf/2111.01084.pdf.
