Process with covariance $E[Y_{t}Y_{s}]=a_{1}-a_{2}|t-s|$ We have a centered Gaussian process $X_{t}$ where we have exact equality $$E[X_{t}X_{s}]=a_{1}-a_{2}|t-s|$$ for $|t-s|<\epsilon_{0}\ll \frac{a_{1}}{a_{2}}$ and $a_{i}>0$.  

Q: I am curious if there is any other concrete Gaussian process $(Y_{s})_{s\in [0,\epsilon_{0}]}$  out there with the same exact covariance when $|t-s|<\epsilon_{0}$ for some $\epsilon_{0}>0$ (not asymptotical behaviour with error term, but exact equality).

It will be interesting if $Y_{t}$ is in terms of some known process like a functional of Brownian motion or a stationary solution of some SDE.
We are not concerned with $Y_{t}$ having different distribution than $X_{t}$(even though they do looked as Gaussian processes over $t\in [0,\epsilon']$ for $\epsilon'$ small enough). Our main concern is if such covariances have been studied in the literature or if we can devise one. 
Some idea: start from $Y_{t}=\int_{0}^{t}f(r,t)dW_{r}$ and try to find a deterministic $f(r,t)$ with the desired covariance: by Ito isometry $\int_{0}^{s}f(r,s+h)f(r,s)ds=a_{1}-a_{2}h$.
Our process
Let $X_{\epsilon}(x)\sim N(0,\ln\frac{1}{\epsilon})$ with covariance:

For simplicity above we suppressed the $\epsilon$ and just let $X_{t}:=X_{\epsilon}(t)$.
Our particular process. Consider the hyperbolic measure $\lambda:=\frac{1}{y^{2}}dx dy$ in the upper half-plane and a White noise process W indexed by Borel sets of finite hyperbolic area: 
$$\{A\subset \mathbb{H}: \lambda(A)<\infty; \sup_{(x,y),(x',y')\in A}|x-x'|<\infty\}$$
with covariance:
$$E[W(A_{1})W(A_{2})]:=\lambda(A_{1}\cap A_{2}).$$
Then let $X_{t}=W(V_{\epsilon}+t)$ for 
$$V_{\epsilon}:=\{(x,y)\in \mathbb{H}: x\in [-1/4,1/4]\text{ and }max(2|x|,\epsilon)\leq y<1/2\}.$$
 A: By Pólya’s theorem, any even real-valued function $f$ on $\mathbb R$ with $f(\infty-)=0$ which is convex on $[0,\infty)$ is positive definite. So, any such function is the (auto)covariance function of a stationary Gaussian process; see e.g. Section "Properties of the Autocovariance Function", page 2. 
Now just take any two different functions, $f_1$ and $f_2$, of the Pólya class such that $f_2(t)=1-|t|=f_2(t)$ for $|t|\le1/2$. Then the corresponding stationary Gaussian processes, say $(X_{1,t})$ and $(X_{2,t})$, with the covariance functions $f_1$ and $f_2$ will have different distributions. Therefore, these two processes will be different from each other. 

To be more specific, note first here that, 
by vertical and horizontal re-scaling, without loss of generality $a_1=a_2=1$, so that
$$EX_sX_t=1-|t-s|\quad\text{if}\quad|t-s|\le u, \tag{1}$$
where $u\in(0,1)$. 
Let then 
$$Y_t:=B_{t+1}-B_t=\int_t^{t+1}dB_s,$$
where $(B_t)_{t\in\mathbb R}$ is the standard Brownian motion with $B_0=0$. Then 
$$EY_sY_t=1-|t-s|\quad\text{if}\quad|t-s|\le 1$$
(with $EY_sY_t=0$ if $|t-s|>1$),
so that 
$$EY_sY_t=EX_sX_t\quad\text{if}\quad|t-s|\le u,$$
as desired. 
For more examples, take any $h\in(0,1)$ and let 
$$U_t:=\frac1{\sqrt2}\,(Y_{(1-h)t}+Z_{(1+h)t}),$$
where $(Z_t)$ is an independent copy of the Gaussian process $(Y_t)$. Then 
$$EU_sU_t=1-|t-s|=EY_sY_t \quad\text{if}\quad|t-s|\le1/(1+h)$$
and hence 
$$EU_sU_t=EX_sX_t \quad\text{if}\quad|t-s|\le\min[u,1/(1+h)],$$
as desired.
A: By vertical and horizontal re-scaling, without loss of generality $a_1=a_2=1$, so that
$$EX_sX_t=1-|t-s|\quad\text{if}\quad|t-s|\le u, \tag{1}$$
where $u\in(0,1)$. 
Take any $h\in(0,1)$ and let 
$$U_t:=\frac1{\sqrt2}\,(X_{(1-h)t}+Y_{(1+h)t}),$$
where $(Y_t)$ is an independent copy of your Gaussian process $(X_t)$. Then 
$$EU_sU_t=1-|t-s|=EX_sX_t \tag{2}$$
if $|t-s|\le u/(1+h)$, as desired. 
Moreover, letting $u$ take the largest possible value such that (1) still holds, we will see that the first equality in (2) will fail to hold if $|t-s|=u$, which shows that the distribution of $(U_t)$ is different from that of $(X_t)$. Therefore, the process $(U_t)$ is different from the process $(X_t)$. 
