The use of $rad/\sqrt{\text{Hz}}$ suggests that this is phase noise specifically (a spectral density due to phase fluctuations), and typically in my use this has been described as a power spectral density (units of $rad^2/\text{Hz}$), so this is just the square root of that quantity.
The reason the DFT (of which the FFT computes) is divided by $N$ is to normalize the FFT to be the same units of the time domain signal, specifically using the following normalized form of the DFT:
$$X_1(k) = \frac{1}{N}\sum_{n=0}^{N-1}x[n]W_N^{nk}$$
Versus the typically version that is not normalized which the FFT returns:
$$X(k) = \sum_{n=0}^{N-1}x[n]W_N^{nk}$$
With such a normalization, the magnitude of $x[n]$ at any specific frequency will match the magnitude of $X(k)$ for that frequency. For example, if we had a time domain waveform of a sinusoidal phase error versus time given as:
$$\phi[n] = A\cos(\omega n) = \frac{A}{2}e^{j\omega n} + \frac{A}{2}e^{-j\omega n} \space \text{rad}$$
Then assuming $\pm\omega$ were exactly on a bin centers (for the DFT due to its circular nature $-\omega = N-\omega$), the resulting two bins in $X_1(k)$ would have a magnitude of $\frac{A}{2}$, matching the magnitudes of the time domain waveform.
As a power spectral density (meaning we are interested in the power over a given frequency range) the normalized power of each frequency index in the DFT (aka bin) is then:
$$|X_1(k)|^2 = \frac{|X(k)|^2}{N^2} \space \frac{\text{rad}^2}{\text{bin}}$$
(Where the units of $\text{rad}^2$ for the power quantity $|X_1(k)|^2$ only make sense if x[n] was the phase noise in units of radians).
$\frac{\text{rad}^2}{\text{bin}}$ is a power quantity per bin. To make this the recognized form of the power spectral density in power/Hz we recognize that $Nd = f_s$ where $N$ is the number of samples in the DFT, $f_s$ is the sampling rate, and $d$ is the spacing of each frequency index (bin as the OP used) in Hz resulting in the spectral width of each bin in Hz:
$$d = \frac{f_s}{N} \space \frac{\text{Hz}}{\text{bin}}$$
Thus
$$ \frac{|X(k)|^2}{N^2} \frac{\text{rad}^2}{\text{bin}} \times d^{-1} \frac{\text{bin}}{\text{Hz}} = \frac{|X(k)|^2}{N^2}\frac{N}{f_s} \frac{\text{rad}^2}{\text{Hz}} = \frac{|X(k)|^2}{N f_s} \frac{\text{rad}^2}{\text{Hz}}$$
This result would specifically be what we typically notate as $\scr{L}_{\phi}(f)$ as the two-sided power spectral density due to phase fluctuations (since the DFT contains both sides of the spectrum, in contrast to the one-sided PSD which is $S_\phi(f) = 2\scr{L}_{\phi}(f)$.).
Note we say "due to phase fluctuations" since the units here were phase. It is also interesting how the phase unit in radians when squared is the power unit relative to the carrier (often expressed as dBc/Hz). This is clear for small angles given the small angle approximation $sin(\theta) \approx \theta$, or geometrically the quadrature component being the noise as phase noise relative to the in-phase component being the carrier that has been rotated due to that phase, such that the ratio of the two is the phase unit in radians, for small angles!) This is why when phase noise is dominant this computation will match the actual power measurement we see under test with a spectrum analyzer.
Further update:
The OP clarified in his comments that his question is specific to the peak at 30 Hz offset as shown in this plot:
It isn't specified but assuming this is a two-sided spectral density, the peak of a single tone would have a total power independent of density, so we would typically report its result as $\text{rad}^2$ and not $\text{rad}/\text{Hz}$ (or the magnitude quantity as the square root $\text{rad}$ as used in this plot, meaning this plot is $\sqrt{\scr{L}_{\phi}(f)}$). The paper also incorporates a moving average of 5 and suggests in a foot note that the peak would be $\approx 1.6 \text{nrad}/\sqrt{1\text{kHz}}/5$, and the plot was scaled (moved up or down) such that the level of the tone landed on this expectation.
I suggest that the peak would be at either $\approx 1.6 \text{nrad}/20$ or $\approx 1.6 \text{nrad} \sqrt{2}/20$ depending on if the spectrum is intended to be double-sided or single-sided which should be specified. The sampling rate does not change the value of the tone on the spectral density when the units are already in nrad, so there should also be no $\sqrt{1\text{kHz}}$ in that answer - The sine wave theoretically occupies zero bandwidth, or for practical reasons we can assume we integrated that power over a small bandwidth in order to measure the peak we see. Either way the density becomes a single figure for the tone independent of bandwidth. Any windowing applied in the time domain prior to the FFT (other than the rectangular window) will also shift the value of the tone differently from the values for the noise. Further details below.
To confirm that assumption, here is my prediction of where such a tone would be:
The 1.6 nrad oscillation is specified as the peak to peak value and thus is of the form:
$$\phi(t) = \frac{1.6}{2} \cos(2\pi f t) \space\space \text{nrad}$$
with $f=30e3$
If the spectrum is two-sided (as $\sqrt{\scr{L}_\phi(f)}$ rather than one-sided as $\sqrt{S_{\phi}(f)}$), then the spectrum is only showing the upper half of this two-sided spectrum, with both sides given by:
$$\phi(t) = \frac{1.6}{2} \cos(2\pi f t) = \frac{1.6}{4}e^{j 2\pi f t} + \frac{1.6}{4}e^{-j 2\pi f t} \space\space\text{nrad}$$
Thus prior to the effect of the moving average filter (MAF), I would predict the tone shown on a double-sided spectrum to be at:
$$\frac{(1.6e-9)}{4} = (4e-10) \space \text{rad}$$
Notice the units are $\text{rad}$ and not $\text{rad}/\sqrt{\text{Hz}}$ as the standard deviation of the tone itself is not a density spread across frequency, unlike that of the noise.
I assume the moving average filter that is mentioned was done on the frequency domain samples. If in the time domain there would be an additional loss of 0.963 but I don't see evidence of such a moving average response in the plot, in which case with a moving average of frequency samples, the tone is reduced by a factor of 5 as the author had done, resulting in $(4e-10)/5 = (8e-11)$.
If the plot was supposed to be a single-sided spectrum $\sqrt{S_{\phi}(f)}$, then the result would be $\sqrt{2}$ larger or $1.13e-10$, which is consistent with the standard deviation of $\phi(t)$ reduced by the MAF.
Neither of these results match the plot, but this is where I would expect a 30 Hz tone after a moving average of 5 samples when sampled at 1 KHz if the units of the spectral density are $\text{nrad}/\sqrt{\text{Hz}}$, for either case of a single-sided or double-sided spectral density. Also note that my computation was independent of the bin size or number of samples since as the author of the paper was intending to do (and perhaps did if I made an error in my prediction) was to predict the expected value of that tone and then scale the plot accordingly. My earlier answer shows how I would scale the result from the DFT directly in which case the bin size and number of samples would be involved.
As a further note since these spectrums are being derived from FFT's and since the OP is interested ultimately in assessing noise: We must also be careful to account for the equivalent noise bandwidth due the effect of windowing especially if we are normalizing the plot based on the power of a tone. (and other effects such as scalloping loss etc which have been minimized by choosing a tone at or near a bin center as was done). Any windowing done on the time domain signal other than the rectangular window will widen the bandwidth of each bin beyond the single bin as given by the rectangular window, which means that the noise measured will be larger than the actual noise! Further the window has a loss reducing the signal from the tone and the noise, but because of the effectively wider noise bandwidth of each bin the noise will go down less than the tone (the tone only occupies one bin)! The effect of the moving average in frequency on SNR is also affected by the window since the adjacent noise bins are no longer uncorrelated. I detail this further in this post: Find the Equivalent Noise Bandwidth
