Inconsistency between the units of power spectral density and the definition that people often give
The OP is correct in their dimensional analysis
$|X(f)|^2$ is NOT the power spectral density, despite what other authors might claim. Other authors probably call this the power spectral density because it is close to right and it captures most of the important features without having to delve into technicalities.


Power has dimensions of $[\text{signal}^2]$. Energy has dimensions of $[\text{power}\cdot\text{time}] = [\text{signal}^2\cdot\text{time}]$.
The spectral density of anything has dimensions of $[\text{thing}\cdot \text{frequency}^{-1}]$. Thus, power spectral density has dimensions of $[\text{signal}^2 \cdot \text{frequency}^{-1}] = [\text{signal}^2\cdot \text{time}]$. Note that it is coincidental that power spectral density has the same dimensions as energy and it should be understood that power spectral density is power per frequency. Also note that the Fourier transform of anything always has dimensions of $[\text{thing}\cdot\text{frequency}^{-1}]$.

The power spectral density is more nicely defined as follows. We define the windowed signal

$$
x_{\Delta t}(t) = \begin{cases}
x(t) \text{ for } |t|&lt; \frac{\Delta t}{2}\\
0 \text{ for } |t| \ge \frac{\Delta t}{2}
\end{cases}
$$

The windowed Fourier transform is then

$$
X_{\Delta t}(f) = \int_{t=-\infty}^{+\infty} x_{\Delta t}(t) e^{-i2\pi f t} dt = \int_{t=-\frac{\Delta t}{2}}^{\frac{\Delta t}{2}} x(t) e^{-i2\pi f t} dt
$$

The power spectral density is then defined by

$$
S_{xx}(f) = \lim_{\Delta t\rightarrow \infty} \frac{1}{\Delta t} |X_{\Delta t}(f)|^2
$$

More properly when dealing with random signals one might take an expectation value of the squared windowed transform.

This can be expressed another way. We can define a window function

$$
w_{\Delta t}(t) = \frac{1}{\sqrt{\Delta t}} \theta\left(t-\frac{\Delta t}{2}\right)\theta\left(\frac{\Delta t}{2} - t\right)
$$

Here $\theta$ is the Heaviside function.
And a windowed version of $x(t)$ given by

$$
x_{w_{\Delta t}}(t) = x(t)w_{\Delta t}(t)
$$

Note that this is the exact same as the windowed function defined above but with a factor of $\frac{1}{\sqrt{\Delta t}}$ built in.
The Power spectral density can then be defined equivalently as

$$
S_{xx}(f) = \lim_{\Delta t \rightarrow \infty} |X_{w_{\Delta t}}(f)|^2
$$

The reason we must work with $x_{w_{\Delta t}}(t)$ rather than $x(t)$ is that $x(t)$ is that, if $x(t)$ has constant power or at least finite power for infinite time, then $x(t)$ has infinite energy. However, even if $x(t)$ has infinite energy, $x_{w_{\Delta t}}(t)$ has finite energy. Note that the window function is not dimensionless but acts so that the finite total energy in $x_{w_{\Delta t}}(t)$ given by $\int |x_{w_{\Delta t}}(t)|^2 dt$ is in fact the average finite energy in $x(t)$.

We also have the fact that infinite length signals do not have well behaved Fourier transforms, for example, the Fourier transform of a pure tone $e^{+i2\pi f_0 t}$ is a dirac delta function, i.e. not well behaved.
The windowed version of this will have a well-behaved Fourier transform.

@Dan Boschen expresses some confusion about reconciling the dimensions of $S_{xx}(f)$ with the Fourier transform of the autocorrelation function. There is no need for confusion. The units agree.

$$
S_{XX}(f) = \tilde{R}_{xx}(f) = \int R_{xx}(t) e^{-i2\pi ft} dt = \int \langle x(t)x(0)\rangle e^{-i2\pi ft}dt
$$

The expression on the right has dimensions of $[\text{signal}^2\cdot \text{time}]$ which is the same as the units of power spectral density expressed above. This should hint that the Fourier transform of the auto-correlation function is NOT given by $|X(f)|^2$...

$R_{xx}(t)$ (for stationary $x(t)$) is defined as

ensemble average:
\begin{align}
R_{xx}(t) = \langle x(t)x(0) \rangle = \int yz f_{x(t),x(0)}(y,z) dy dz
\end{align}

$f_{x(t),x(0)}(y,z)$ is the joint probability density function for the random variables $x(t)$ and $x(0)$ so it has dimensions of $[\text{signal}^{-2}]$.

time average:
\begin{align}
R_{xx}(t) = \langle x(t)x(0) \rangle = 
\lim_{\Delta t \rightarrow \infty} \frac{1}{\Delta t} \int_{t'=-\frac{\Delta t}{2}}^{\frac{\Delta t}{2}} x(t'+t)x(t') dt'
\end{align}