Since the Spearman correlation is just the Pearson correlation of the ranks of the values, an intermediate step in finding the distribution of this statistic is to obtain the distribution of the rank vector for the multivariate normal distribution arising in the Gaussian $\text{AR}(1)$ process. I am not aware of an analytic result for this problem, and I strongly suspect that the distribution would not have a closed form solution. In any case, let's examine this problem.
Suppose we let $\mathbf{X} = (X_1,...,X_n)$ be a time-series vector from a stationary Gaussian $\text{AR}(1)$ process with auto-correlation parameter $\phi$, mean $\mu$ and error variance $\sigma^2$. We have the distribution $\mathbf{X} \sim \text{N}(\mu \mathbf{1}, \sigma^2 \mathbf{\Phi})$ where $\mathbf{\Phi} \equiv [\phi^{|i-j|}]_{i,j}$ is the exponentially decaying correlation matrix. Since the distribution of the ranks of the time-series values is invariant to the location $\mu$ and the scale $\sigma$, we can set $\mu = 0$ and $\sigma = 1$ without loss of generality in the analysis. (For the remainder of the analysis we assume these parameter values.)
From the underlying time-series vector we can define the ranking vector $\mathbf{R} = (R_1,...,R_n)$, which is distributed over all permutations of $\{ 1,...,n \}$. The probability mass function for this random vector is:
$$\begin{align}
p_{\mathbf{R}}(\mathbf{r})
&\equiv \mathbb{P}(\mathbf{R} = \mathbf{r}) \\[18pt]
&= \mathbb{P}(X_{r_1} \leqslant \cdots \leqslant X_{r_n}) \\[18pt]
&= \underset{x_{r_1} \leqslant \cdots \leqslant x_{r_n}}{\int \cdots \int} \ \text{N}(\mathbf{x} | \mathbf{0}, \mathbf{\Phi}) \ d \mathbf{x} \\[6pt]
&= \frac{1}{(2 \pi)^{n/2} |\mathbf{\Phi}|^{1/2}} \underset{x_{r_1} \leqslant \cdots \leqslant x_{r_n}}{\int \cdots \int} \ \exp \Bigg( - \frac{1}{2} \cdot \mathbf{x}^\text{T} \mathbf{\Phi}^{-1} \mathbf{x} \Bigg) \ d \mathbf{x}. \\[6pt]
\end{align}$$
This is a complicated integral and it does not appear to have any closed-form representation. Consequently, you would need to compute this integral numerically. The simplest way to deal with this problem is probably just to generate simulations of the random vector $\mathbf{X}$ (e.g., using the rGARMA function in the ts.extend package) and then compute the rank vector for each simulation, using the empirical distribution of the rank vector to approximate its true distribution.
Now, in your problem you want to find the Spearman correlation between two disjoint parts of the time-series. Suppose you use disjoint parts with a gap of size $\ell \geqslant 0$ and denote the rank vectors for the two parts as $\mathbf{R} = (R_1,...,R_n)$ and $\mathbf{R'} = (R'_1,...,R'_n)$, with the latter being the later vector. The joint distribution of these two rank vectors has probability mass function:
$$\begin{align}
p_{\mathbf{R}, \mathbf{R}'}(\mathbf{r}, \mathbf{r}')
&\equiv \mathbb{P}(\mathbf{R} = \mathbf{r}, \mathbf{R}' = \mathbf{r}') \\[18pt]
&= \mathbb{P}(X_{r_1} \leqslant \cdots \leqslant X_{r_n}, X_{n + \ell + r'_1} \leqslant \cdots \leqslant X_{n + \ell + r_n}) \\[18pt]
&= \underset{\begin{matrix} x_{r_1} \leqslant \cdots \leqslant x_{r_n} \\ x_{n + \ell + r'_1} \leqslant \cdots \leqslant x_{n + \ell + r_n} \end{matrix}}{\int \cdots \int} \ \text{N}(\mathbf{x} | \mathbf{0}, \mathbf{\Phi}) \ d \mathbf{x} \\[6pt]
&= \frac{1}{(2 \pi)^{n/2} |\mathbf{\Phi}|^{1/2}} \underset{\begin{matrix} x_{r_1} \leqslant \cdots \leqslant x_{r_n} \\ x_{n + \ell + r'_1} \leqslant \cdots \leqslant x_{n + \ell + r_n} \end{matrix}}{\int \cdots \int} \ \exp \Bigg( - \frac{1}{2} \cdot \mathbf{x}^\text{T} \mathbf{\Phi}^{-1} \mathbf{x} \Bigg) \ d \mathbf{x}. \\[6pt]
\end{align}$$
This distribution also does not appear to have a closed form representation, so the distribution of the resulting Spearman correlation with the time index is likely to be highly complicated and also would not have a closed form representation. Again, the simplest way to determine the distribution would be via simulation.
Simulating the distribution of the Spearman/Pearson correlation: For completeness, I will give an example of simulating the distribution of the Spearman and Pearson correlations for disjoint parts of the time-series in this kind of process. Here is a function to simulate from the stationary Gaussian $\text{AR}(1)$ process and compute the Spearman/Pearson correlations of the simulations. The output is a list containing the simulations, rank vectors and correlations.
#Create function to simulate the distribution of Spearman/Pearson correlation
simulate <- function(sims = 10^7, n, gap, phi) {
#Set rank matrix and Spearman/Pearson vectors
RANKS1 <- matrix(0, nrow = sims, ncol = n)
RANKS2 <- matrix(0, nrow = sims, ncol = n)
SPEARMAN <- rep(0, sims)
PEARSON <- rep(0, sims)
#Generate simulations of AR(1) process and compute rank vectors and correlation
SIMS <- ts.extend::rGARMA(n = sims, m = 2*n+gap, ar = phi)
for (i in 1:sims) {
RANKS1[i, ] <- order(SIMS[i, 1:n])
RANKS2[i, ] <- order(SIMS[i, n + gap + 1:n])
SPEARMAN[i] <- cor(RANKS1[i,], RANKS2[i,])
PEARSON[i] <- cor(SIMS[i, 1:n], SIMS[i, n + gap + 1:n]) }
#Give output
list(simulations = SIMS, ranks1 = RANKS1, ranks2 = RANKS2,
spearman = SPEARMAN, pearson = PEARSON) }
Running this function for a sufficiently large number of simulations (relative to the vector size $n$) allows us to simulate the true distribution of the Spearman correlation. I will give an example for time-series parts with $n=10$ consecutive elements with a gap of $\ell = 6$ elements, using the auto-regression parameter $\phi = 0.8$. If you want to simulate for higher values of $n$ of $\ell$ then this may give you an indication of how the respective distributions of the correlations change as these values increase.
#Set parameters
n <- 10
gap <- 6
phi <- 0.8
#Simulate the Spearman correlation
set.seed(1)
SIMULATIONS <- simulate(n = n, gap = gap, phi = phi)
#Show histograms of the simulated Spearman/Pearson correlation
par(mfrow = c(2, 1))
hist(SIMULATIONS$spearman, col = 'blue', xlim = c(-1, 1), freq = FALSE,
main = NULL, xlab = 'Spearman Correlation', ylab = 'Probability')
hist(SIMULATIONS$pearson, col = 'red', xlim = c(-1, 1), freq = FALSE,
main = NULL, xlab = 'Pearson Correlation', ylab = 'Probability')
