Power spectral density: Why are these two methods equal?
Regarding what we discussed in the comments, here is an explanation of the different concepts I was referring to:

The concrete usage of terms may vary, but my interpretation is this: For a stochastic process $x(t)$, the autocorrelation function is defined as $$r(\tau) = \mathbb{E}\{x(t) \cdot x(t+\tau)\}.$$ Note that the process needs to be stationary, otherwise the autocorrelation function depends on $t$ as well. Its power spectral density is given by $$R(f) = \mathbb{E}\{|X(f)|^2\}.$$ These two quantities are linked to each other via the Wiener-Khinchin-Theorem which shows that $$R(f) = \int_{-\infty}^{\infty} r(\tau) {\rm e}^{-\jmath 2\pi f \tau} {\rm d} \tau,$$ i.e., they are Fourier transform pairs. 

Since the expectation is hard to implement in practice (we typically do not have an infinite ensemble available), we need to resort to other methods. When the process is ergodic, ensemble average can be replaced by temporal average and we may, e.g., estimate the autocorrelation function from a window $T$ via $$r_T(\tau) = \frac{1}{T} \int_T x(t) \cdot x(t+\tau){\rm d}t.$$ Note that $r(\tau) = \lim_{T\rightarrow \infty} r_T(\tau)$. On the other hand, considering using $R_T(f) = |X_T(f)|^2$ as an estimate for the power spectral density (referred to as the periodigram), where $X_T(f)$ is the spectum of a $T$-windowed signal leads to an estimate that does not converge since $S(f) \neq \lim_{T\rightarrow \infty} S_T(f)$.

Now, I'm not sure what exactly is the question. Is it the explanation of the Wiener-Khinchin-Theorem for the stochastic signals? Or is it the deterministic equivalent as, e.g., discussed in this answer?

*edit: Regarding your clarification that you posted as a reply to your question(?): My feeling is this question has little to do with ACF and PSD being Fourier transform pairs, instead you seem to be worried how to normalize the result of a Discrete Fourier transform (carried out via an FFT) such that its units make sense physically. It's probably worth asking this in a separate question, but the short version is here. There are two things to be considered:


There are different conventions how to normalize DFT and IDFT. Let $\mathbf{F}$ be a DFT matrix with elements ${\rm e}^{-\jmath \frac{mn}{N}}$, then you can either (i) define DFT via $\mathbf{D} = \mathbf{F} \cdot \mathbf{d}$ in which case you need to define the IDFT via $\mathbf{d} = {\color{red}{ \frac{1}{N} }} \mathbf{F}^H \cdot \mathbf{D}$ or (ii) define the DFT via $\mathbf{D} = {\color{red}{ \frac{1}{N} }}\mathbf{F} \cdot \mathbf{d}$ which leads to the IDFT $\mathbf{d} =  \mathbf{F}^H \cdot \mathbf{D}$ or alternatively (iii) you let the DFT be $\mathbf{D} = {\color{red}{ \frac{1}{\sqrt{N}} }}\mathbf{F} \cdot \mathbf{d}$ which leads to the IDFT $\mathbf{D} = {\color{red}{ \frac{1}{\sqrt{N}} }} \mathbf{F}^H \cdot \mathbf{d}$. Different conventions are in place, Matlab uses the first by default, whereas Maple will prefer the third.
The Fourier transform changes the units: the ACF has a unit of power (say, Watt), the Fourier transform is an integral over time, hence it adds seconds, so that the unit is Watt seconds or better Watt/Hz. It's a density anyways (how do Watts distribute over Hertz?). If you want to approximate a Fourier transform via a Discrete Fourier transform, you have to be super careful, since the DFT/IDFT does not change the units. If your input to the DFT is a sequence with unit, say, Watt, its DFT coefficients will have the same unit. If you want to mirror what happens in a Fourier transform, you need to multiply the result with the sampling period of the samples in time domain. In other words: If you have a function function $u(t)$ with spectrum $U(f)$ and you want to approximate $U(f)$ with a DFT, you need to consider a finite window (say, length $T$) and a finite sampling interval (say, $t_0 = T/N$). Then, ${\color{red}{t_0}}\cdot D[\mu]$ is an approximation for $U(\mu f_0)$ where $f_0 = 1/T$. It's an approximation since sampling as well as truncation can cause errors (one of them will for sure).


I think the factor you were looking for can very likely be a combination of these two effects, your $1/(f_s N)$ is the same as my $t_0/N$.