What does $P(|X_n - X| \geq \epsilon)$ represent intuitively? I get that $P(|X_n - c| \geq \epsilon)$ represents the probability that the random variable $X_n$ is outside the interval of $(c - \epsilon, c + \epsilon)$ but I am not sure how it works with a random variable where $P(|X_n - X| \geq \epsilon)$. How can a random variable take only a single value like a constant? How would you look at this intuitively with the constant $c$ case?
 A: I think this is a case where thinking of $X_n$ and $X$ as functions is helpful.
Suppose we had two functions $f,g : \mathbb R\to\mathbb R$ and maybe they look like this:

In that plot I marked on the x-axis the regions where $|f-g| > .75$, where $0.75$ was chosen arbitrarily. We can then think of the total length of the regions of disagreement as a sense of how different $f$ and $g$ are. If there are long intervals where $f$ and $g$ are more than $0.75$ apart, we'll have a big disagreement, while if there are only small regions of disagreement then $f$ and $g$ are pretty similar.
We are doing the exact same thing in $\newcommand{\e}{\varepsilon}P(|X_n-X|>\e)$: by definition,
$$
P(|X_n-X|>\e) = P\left(\{\omega\in\Omega : |X_n(\omega) - X(\omega)| > \e\}\right).
$$
Here, the sample space $\Omega$ is playing the role of the domain in my example with $f$ and $g$, and instead of just looking at the length of the region of disagreement (which would be using the Lebesgue measure), we measure it with $P$. But fundamentally the big idea is that we're measuring the region of the domain where $X_n$ and $X$ disagree by more than our threshold of $\e$.
If one of the functions happens to be constant it's the exact same idea, just simpler since only one thing is varying now.

Here's the code for that plot if you want to experiment. I sampled $f$ and $g$ iid from a Gaussian process with a squared exponential kernel to get smooth arbitrary-looking functions.
set.seed(2)
n <- 1000
xseq <- seq(-5, 5, length=n)
K <- exp(-.5 * as.matrix(dist(xseq))^2)
sims <- MASS::mvrnorm(2, rep(0,n), K)
plot(sims[1,] ~ xseq, type="l", col=2, lwd=2, xlab="x", ylab="y",
     main="f and g with |f-g|>.75 marked", ylim=c(-2,2))
lines(sims[2,] ~ xseq, col=4, lwd=2)
legend("topleft", c("f", "g"), lwd=2, col=c(2,4), bty="n")

thresh <- .75
rug(xseq[abs(sims[1,] - sims[2,]) > thresh])

