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.

It is not very difficult to find an unbiased estimator for $P(\theta > T)$ : you want to find some function $g(\theta ; T)$ that has the property that $\mathbb E[ g(\theta;T)] = P(\theta > T)$. But $P(\theta > T) = \mathbb E[ 1(\theta > T)]$, so clearly $g(\theta;T) =  1(\theta > T)$ does the job.
Furthermore this is the only unbiased estimator that can be constructed with one sample:
you can write the above unbiasedness condition as the equality
$$ \int d\theta \left[ g(\theta;T) -  1(\theta > T) \right ]f(\theta|\mu,\sigma^2) = 0$$
which need to hold for all values of $\mu$ (where $f(\theta|\mu,\sigma^2)$ is the pdf of the normal distribution). Clearly this can happen only if the integrand itself is zero.
Of course this is not a very good estimator : the estimate is one if $\theta > T$ and zero otherwise. But this just demonstrate the general potential difficulties when dealing with frequentist quantities - unbiasedness doesn't guarantee anything other than unbiasedness.
Note that if you want instead a confidence interval for  $P(\theta > T)$, you can immediately get it from the standard confidence interval for $\mu$, since $P(\theta > T) = \Phi(\frac{T-\mu}{\sigma})$ is a monotonically decreasing function of $\mu$. So just mapping the interval using the normal CDF will maintain the coverage property.
