Having problems with understanding the meaning of PSD for colored Gaussian noise
Think in terms of noise realizations. When your head is full to the brim of covariances and multivariate Gaussian variables, generate a Gaussian white noise realization of length N, and create, using this sequence, a sequence $y(n) = x(n) - x(n-1)$. Let us call the latter a blue Gaussian noise realization. I plot histograms of these sequencies which show how many values fall into each bin, those bins uniformly partitioning a range of $y, x$ values.

In compliance with your convention, a blue colored graph represents the white noise data, a red colored, the blue noise data. This histograms visually agree with your observation that

$w_n \sim N(0,2σ_n^2)$, 2x higher variance than $x_n$

Each histogram plotted here is an approximate representation of the distribution of noise signal data, id est, the Probability Density Function of noise signal. The realization sequence's length N is 220 samples. The longer sequence length gets and the greater bin size gets, the closer these numerical PDF graphs resemble theoretical curves.
Go on to Power Spectral Density. Here are discrete Fourier transforms (actually, FFTs) of our white-noise and blue-noise realizations, squared and plotted. Color scheme according to your convention (white noise in blue, blue noise in red):

A white-noise PSD (blue graph) is near-constant, a blue-noise PSD graph (red graph) is two symmetrical parabolic segments sewn at 1/4 and 3/4 of the entire range (the shapes are revealing $\omega$ squared function and FFT mirroring). The high-frequency end is in the middle of horizontal axis.
No relationship between variance and PSD is revealed in our calculations and pictures. It appears that these concepts are things that can only be related indirectly, using formulas that have yet to be derived at this stage of learning.
Variance 'measures how far a set of numbers is spread out from their average value' (Wikipedia). Variance is defined in the time/space domain, and this definition never mentions frequencies or wavenumbers or Fourier transform. Although variance can be measured for any &quot;set of numbers&quot;, in many signal processing scenarios said set of numbers is almost synonymous with sampled noise values.
A PSD thing, on the contrary, belongs in the frequency domain. The article in Wikipedia states from the beginning

According to Fourier analysis, any physical signal can be decomposed
into a number of discrete frequencies

and continues

PSD of the signal describes the power present in the signal as a
function of frequency. PSD of noises is a specific case of PSD of
signals in general.

So variance and PSD are sort of &quot;orthogonal&quot; concepts; whatever, variance cannot be a function of frequency in any meaningful way, while PSD is.
Compare also the histograms of the white-noise realization PSD (top) and the blue-noise realization PSD (bottom). The white-noise PSD shape replicates, as it should, a normal distribution of source signal values. A peak of the blue-noise PSD shape shows greater share of the high-frequency components in the blue noise signal.

PS There is a number of Guassian Process tutorials on github and elsewhere, mostly Jupyter notebooks. See, for example, https://peterroelants.github.io/posts/gaussian-process-kernel-fitting/.