Theory
------
In the limit $f/f_s \to 0,$ where $f$ is the sine wave frequency and $f_s$ is the sampling frequency, the answer to the question will be found in the [Fourier series][1] of the continuous-time piece-wise constant quantized waveform. We can construct the waveform as a sum of components like the red curve here illustrated for sine wave amplitude $N = 7$:
[![enter image description here][2]][2]<br>*Figure 1. limit $f/f_s \to 0$ waveform, of amplitude $N=7$, and its quantization in blue, with a component highlighted in red.*
We are working with cosine instead of a sine because the math is nicer this way. Because of the symmetries, each component has only odd harmonics and only cosine terms in its Fourier series. The non-zero coefficients of the Fourier series are given by:
$$a_n = \frac{4}{\pi}\int_{x_0}^{x_1}\cos(nx)\,k\,dx = \frac{4k}{\pi n}\big(\sin(x_1n) - \sin(x_0n)\big),\quad n\text{ odd},\tag{1}$$
with integer amplitude $k \in 1\ldots N$ in range $x_0 < x < x_1$ in the first quarter-period of the cosine. I did not bother to explicitly write the contributions from the symmetrically arranged pieces in the other quarters, because they will contribute identically to the odd harmonic cosine terms in the Fourier series. Instead I simply included the implicit factor 4 in the equation.
The Fourier series of the full piece-wise constant waveform is a sum of the Fourier series of the components. The boundaries of the pieces that we need to include in the sum are:
$$\begin{gather}0 < x < \operatorname{acos}\left(\frac{N - 0.5}{N}\right),\quad \text{if }k = N,\\
\operatorname{acos}\left(\frac{k + 0.5}{N}\right) < x < \operatorname{acos}\left(\frac{k - 0.5}{N}\right),\quad \text{if }k \in 1\ldots N-1,\end{gather}\tag{2}$$
where $\operatorname{acos}\left(\frac{s}{N}\right)$ comes from solving $N\cos(x) = s$ in the first quadrant. The non-zero coefficients of the Fourier series of the full piece-wise constant waveform are then given by:
$$\begin{gather}a_n = \frac{4N}{\pi n}\sin\bigg(\operatorname{acos}\left(\tfrac{N - 0.5}{N}\right)n\bigg) +\\\frac{4}{\pi n}\sum_{k=1}^{N-1} \Bigg(k\sin\bigg(\operatorname{acos}\left(\tfrac{k - 0.5}{N}\right)n\bigg) - k\sin\bigg(\operatorname{acos}\left(\tfrac{k + 0.5}{N}\right)n\bigg)\Bigg),\quad n\text{ odd}.\end{gather}\tag{3}$$
16-bit quantization
-------------------
It is not out of reach to compute Eq. 3 for $N=2^{23} - 1$ and $n \in {5, 7},$ as parameterized in the question, using something like Python's `mpmath`. Let's try that but with a 16-bit sine wave first:
import mpmath as mp
def a_n(N, n):
return 4*N/(mp.pi*n)*mp.sin(mp.acos((N-0.5)/N)*n) + 4/(mp.pi*n)*mp.nsum(lambda k: k*( mp.sin(mp.acos((k-0.5)/N)*n) - mp.sin(mp.acos((k+0.5)/N)*n) ), [1, N-1])
N = 2**15-1 # amplitude (integer)
n = 5 # number of the harmonic (odd integer)
mp.mp.prec = 53 # default precision
20*mp.log10(mp.fabs(a_n(N, n))/N) # Fourier series term amplitude in dB
mp.mp.prec = 106 # Compute again with double the default precision, see if we get the same result
20*mp.log10(mp.fabs(a_n(N, n))/N) # Fourier series term amplitude in dB
mp.mp.prec = 53 # Restore default precision
I'm only reporting as many digits as are agreed about by repeated computation using two different `mpmath` precision settings, including also the digits that got changed by $\pm1$ due to larger changes in later digits. For a 16-bit full-scale $N = 2^{15}-1$ sine wave and $f/f_s \to 0$, the result is **-154.66449 dBFS** for the 5th harmonic and, after setting `n = 7` in the script, **-154.66506 dBFS** for the 7th harmonic, in about 15 seconds of computation on my PC for a single result with the higher precision setting. As a sanity check, I also computed the the amplitude of the 1st harmonic, which is interestingly 0.0000001605530 dBFS (or 1.605530E-7 dBFS in scientific notation) compared to the peak value $2^{15}-1$ at 0 dBFS. I think the rounding at the top of the sine wave is "pulling it up", giving it a higher amplitude of the fundamental.
Let's have a closer look at the decay of the amplitudes of the harmonics:
for m in range(1, 20):
n = mp.mpf(10)**m+1 # number of the harmonic
mp.mp.prec = 70 # precision
a_n_low_prec = a_n(N, n)
mp.mp.prec = 90 # precision
a_n_high_prec = a_n(N, n)
print(str(n)+","+str(20*mp.log10(mp.fabs(a_n_low_prec)/N))+","+str(20*mp.log10(mp.fabs(a_n_high_prec)/N)))
[![enter image description here][3]][3]<br>*Figure 2. With $N=2^{15}-1$, a selection of amplitudes of harmonics of the limit $f/f_s \to 0$ waveform, in dBFS.*
At very large numbers of the harmonic, the factor $\frac{1}{n}$ in Eq. 3 will start to take effect in the amplitude of the harmonics, making the asymptotic decay 20 dB/decade (~6 dB/octave) as function of the number of the harmonic, $n$. The graph is not very smooth, and it would probably look worse if we were to plot all points. I think this is because the waveform is similar to a sawtooth riding on a sine wave, and near the zero crossing the sawtooth frequency stays the same over a large portion of the waveform. Actually, the sawtooth frequency is $2\pi N$ times the fundamental frequency, or approximately the 2.05881E5th harmonic. Indeed, its amplitude turns out as -135.49405259 dBFS, probably much higher than that of a "random" early harmonic.
Aliasing
--------
As a side note, what we have learned from the results already is that for such a high $N$ the amplitudes of the early harmonics do not really decay as function of the number of the harmonic, so with rational $f/f_s,$ higher harmonics aliased to a lower harmonic are likely to be a majority contributor to the amplitude of the lower harmonic. As an example of such aliasing, with $f_s =$ 48 kHz and $f =$ 1000 Hz, the 5th harmonic gets added to by the 43th, 53th, 91th, 101th, 139th, 149th, etc. harmonic:
$$\begin{eqnarray}
&& \ldots\\
&=& 48000\text{ Hz}\times -3 + 1000\text{ Hz}\times 149\\
&=& 48000\text{ Hz}\times -2 + 1000\text{ Hz}\times 101\\
&=& 48000\text{ Hz}\times -1 + 1000\text{ Hz}\times 53\\
&=& 1000\text{ Hz}\times 5\\
&=& 48000\text{ Hz}\times 1 - 1000\text{ Hz}\times 43\\
&=& 48000\text{ Hz}\times 2 - 1000\text{ Hz}\times 91\\
&=& 48000\text{ Hz}\times 3 - 1000\text{ Hz}\times 139\\
&=& \ldots\end{eqnarray}\tag{4}$$
The effect of the addition of the aliases depends also on the phase of the sinusoid, because the phase of the harmonics controls whether there is constructive or destructive interference. I will not analyze the effect of phase in detail, because such a computation is better done by a less theoretical approach.
24-bit quantization
-------------------
Then back to limit $f/f_s \to 0$ free of aliasing. The following continuation of the earlier Python script computes in high precision your desired numbers for $N=2^{23} - 1$, in the limiting case $f/f_s \to 0$:
N = 2**23-1 # amplitude (integer)
n = 5 # number of the harmonic (odd integer)
mp.mp.prec = 70 # Precision 1
20*mp.log10(mp.fabs(a_n(N, n))/N) # in dBFS
mp.mp.prec = 140 # Precision 2
20*mp.log10(mp.fabs(a_n(N, n))/N) # in dBFS
n = 7 # number of the harmonic (odd integer)
mp.mp.prec = 70 # Precision 1
20*mp.log10(mp.fabs(a_n(N, n))/N) # in dB
mp.mp.prec = 140 # Precision 2
20*mp.log10(mp.fabs(a_n(N, n))/N) # in dB
mp.mp.prec = 53 # Restore default precision
After a total of a few hours of computation on my PC, the result is that with $N=2^{23} - 1$ and $f/f_s \to 0,$ the amplitude of the 5th harmonic is **-226.91150085 dBFS** and the amplitude of the 7th harmonic is **-226.9115030 dBFS**. I also computed the amplitude of the 1st harmonic, 3.9195785E-11 dBFS.
[1]: https://en.wikipedia.org/wiki/Fourier_series
[2]: https://i.sstatic.net/1a4gD.png
[3]: https://i.sstatic.net/t29Mg.png