Prediction interval: but instead, the probability that the next datum is above a fixed threshold? I've been struggling with this problem, and I think I must be missing some important conceptual step. Imagine we observe $\theta_1 \sim N(\mu, \sigma^2)$, with unknown $\mu$ and known $\sigma^2$ (for simplicity). We will draw a second datum, $\theta_2 \sim N(\mu, \sigma^2)$ (for now, unobserved). (Point of clarification—$\theta$ is just a random variable here, not a stand-in for the parameter $\mu$, it might have been more clear if I wrote $X_1,X_2 \sim N(\mu, \sigma^2)$.)
It is straightforward to compute a "prediction interval" for $\theta_2$. By definition, this is a frequentist procedure that guarantees that $\theta_2$ is captured by the interval the correct proportion of the time. Note that $\theta_2 - \theta_1 \sim N(0, 2\sigma^2)$. Thus, $P(\theta_2 > \theta_1 - z_\alpha(\sqrt{2} \sigma) )= \alpha $ (where $z_\alpha$ is the z-score, e.g. $z_{.10} = 1.28$ gives a 90% one-sided prediction interval). (See Wikipedia for details).
However, I am interested in a slightly different property. For some fixed threshold $T$, I want to estimate the probability that $\theta_2 > T$. That is, rather than a procedure that creates an interval with "90% coverage", I want a procedure that estimates the probability that the next datum will fall within a fixed interval, $\theta_2 \in [T, \infty]$. Presumably, the "guarantee" I want is that the average probability of this estimate matches the true probability, $P(\theta_2 > T)$.
I'm stumped! Given that $\theta_2 - \theta_1 \sim N(0, 2\sigma^2)$, my instinct is just to compute  $P(Z > (T - \theta_1)/(\sqrt{2}\sigma))$, where $Z \sim N(0,1)$ (i.e. just use the normal CDF). But under simulation, that isn't right. I assume the issue is that $\theta_1$ and $\theta_2 - \theta_1$ aren't independent? But it feels like this problem is doable, and I can't figure it out.

A maximum likelihood estimate of the probability $p_T = P(Y\geq T)$ can be derived using the estimate of the mean $\hat{\mu}$
$$\hat p_T = 1-\Phi\left( \frac{T-\hat\mu}{\sigma}\right) $$
(You have a factor $\sqrt{2}$ which I find not intuitive, with this formula you just reparameterize the distribution and describe it in terms of $\sigma$ and $p_T$ instead of $\sigma$ and $\mu$)
The expectation value of this is
$$\begin{array}{}E\left[\hat p_T)\right] &=& E\left[1-\Phi\left( \frac{T-X_1}{\sigma}\right)\right]\\
&=&1- \int_{-\infty}^{\infty} \Phi\left(\frac{T-\mu}{\sigma} - x\right)\phi(x) dx \\&= &1 - \Phi\left(\frac{T-\mu}{\sqrt{2}\sigma}\right)
\end{array}$$
(For that final equality see https://en.m.wikipedia.org/wiki/List_of_integrals_of_Gaussian_functions )
and it will thus be biased.
This bias depends on the value of $\mu$ and can not be corrected.
