The scaling of $\frac{1}{N}$ satisfies Parseval's theorem. The proof is pretty straightforward. We wish to show
\begin{equation}
\sum_{n=0}^{N-1} \lvert x[n] \rvert^{2} = \frac{1}{N} \sum_{k=0}^{N-1} \lvert X[k] \rvert^{2}
\end{equation}
We can then say
\begin{align}
\sum_{k=0}^{N-1} \lvert X[k] \rvert^{2} &= \sum_{k=0}^{N-1}\sum_{m=0}^{N-1}x[m]e^{-j2\pi\frac{km}{N}}\sum_{l=0}^{N-1}x^{*}[l]e^{j2\pi\frac{kl}{N}} \\ &= \sum_{m=0}^{N-1}x[m]\sum_{l=0}^{N-1}x^{*}[l]\sum_{k=0}^{N-1}e^{j2\pi\frac{(l-m)k}{N}}
\end{align}
The last term is a geometric series that is equal to
\begin{align}
\sum_{k=0}^{N-1}e^{j2\pi\frac{(m-l)k}{N}} &= \frac{1 - e^{j2\pi(l-m)}}{1-e^{j2\pi\frac{(l-m)}{N}}} \\ &= \delta_{l,m}
\end{align}
Thus, computing
\begin{equation}
\sum_{m=0}^{N-1}\sum_{l=0}^{N-1}\sum_{k=0}^{N-1}e^{j2\pi\frac{(m-l)k}{N}} = N\delta_{l,m}
\end{equation}
Plugging this back in we get
\begin{align}
\sum_{k=0}^{N-1}\lvert X[k] \rvert^{2} &= N\sum_{m=0}^{N-1}x[m]x^{*}[m] \\ &= N\sum_{m=0}^{N-1}\lvert x[m] \rvert^{2}
\end{align}
Parseval's theorem is very closely related to the periodogram estimate of the power spectral density (PSD). The power spectral density is defined as
\begin{equation}
\phi(\omega) = \lim_{N\to\infty}E\{\frac{1}{N}\lvert\sum_{t=0}^{N}y(t)e^{-j\omega t}\rvert^{2}\}
\end{equation}
The periodogram simply drops the expected value to give (in the discrete case)
\begin{equation}
\phi[k] = \frac{1}{N}\lvert\sum_{n=0}^{N-1}x[n]e^{-j2\pi\frac{kn}{N}}\rvert^{2}
\end{equation}
The $\frac{1}{N^{2}}$ comes into play with something called the power spectrum. There is much confusion over the difference between the two. The power spectrum (PS) $S$ is related to the PSD by
\begin{equation}
S(\omega_{1-2}) = \frac{1}{\pi}\int_{\omega_{1}}^{\omega_{2}}\phi(\theta)d\theta
\end{equation}
What this means is that the PSD gives the average power (note the expected value in the definition) in the signal at frequency $\omega$, the PS accumulates the average power over a finite bandwidth to give the band-limited power. In the discrete world, there is no continuous spectrum, so the PSD accumulates the power from a small bandwidth and then normalizes by the bandwidth to estimate the average power at the center frequency $2\pi\frac{k}{N}$. The PS undoes this normalization. For a digital filter with temporal aperture $N$, the bandwidth $\beta$ is roughly $\frac{1}{N}$. So, we have
\begin{equation}
S[k] = \beta\phi[k] = \frac{1}{N}\frac{1}{N}\lvert X[k]\rvert^{2} = \frac{1}{N^{2}}\lvert X[k]\rvert^{2}
\end{equation}
If there is windowing or overlapping, this will affect the bandwidth of the filter and requires an alternate scaling, which the reference mentioned by @Jdip in the comments gives. It's also important to note that a one-size-fits-all scaling doesn't apply to adaptive-bandwidth or parametric spectral estimators.
The PS, in my opinion, is also not best represented as a "spectrum" (it's not really a continuous function), but that's a discussion for another time, as in the discrete world it's really just splitting hairs.
EDIT: Updated the geometric series formula, which was wrong. Jeeahmuhtry hard. It should be the correct formula now.
