Evaluating discrete spectral density at only a few frequencies
I believe the OP's question is simplified to the following (confirming I didn't actually miss the salient question):


  Given we can compute a power spectral density from the DFT using the
  conjugate product as:
  
  $$S(k) = |X(k)|^2 = X(k)X^*(k)$$
  
  Can we use a more efficient algorithm, such as the Goertzel, to
  compute a subset of $S(f)$ when only a few points are needed. The OP
  also wanted to understand if there was a limitation with using the
  Goertzel specifically, especially when only the lowest frequencies are of interest.


The answer is yes we can use a more efficient algorithm and specifically we can approximate DFT bins directly using the Goertzel. A windowing function is also recommended in either case (full DFT) or with the Goertzel to minimize aliasing issues (spectral leakage between bins). 

The Goertzel is more efficient than the FFT when $M&lt;log_2(N)$ where N is the total number of samples and $M$ is the total number of bins to compute. Other advantages are we can center each bin of interest at exact frequency values. 

The Goertzel algorithm does have an error bound proportional to $N^2$ that is most pronounced near $\omega = 0$ and $\omega= 2\pi$. An alternate approach that would be identical to FFT results is to compute the DFT bins directly using the DFT formula, with the windowing suggested above also shown:

With windowing applicable to any approach: $x_{win}(n) = w(n)x(n)$, where $w(n)$ is a window, such as Kaiser, Hamming, Blackman, etc.

$$X(k) = \sum_{n=0}^{N-1}x_{win}(n)e^{-j2\pi nk/N} \tag{1}\label{1}$$

Which assuming a complex $x_w(n)$ requires a total of $4N$ real multiplies and $4N$ real adds for each bin computed.  In comparison, the FFT would compute all $N$ bins with $2Nlog_2(N)$ real multiplies and $2Nlog_2(N)$ additions.  

Therefore the total number of bins where it there is savings in computing directly using the DFT is given by $M$ for:

$$4NM &lt; 2Nlog_2(N)$$

Resulting in:

$$M &lt; \frac{log2(N)}{2}$$

So for example, if we had a 1024 point DFT, if less than $log_2(1024)/2 = 5$ bins were needed, it would be more efficient to compute them with the DFT equation directly rather than use the FFT. The Goertzel extends this threshold to 10 bins.

Further, as we can in the Goertzel, we can modify the DFT equation to center the bin on any frequency $\omega_o$ for $\omega_o \in [0, 2\pi)$ as follows:

$$X(\omega_o) = \sum_{n=0}^{N-1}x_w(n)e^{-j\omega_o n} \tag{2} \label{2}$$

In the DFT equation, $\omega$ is limited to discrete frequencies given by $k\omega_o$ where $\omega_o = 2\pi n/N$, thus it should be clear how $\ref{2}$ and $\ref{1}$ are related and how the form in $\ref{2}$ allows us to place the bin on any frequency along the continuous $\omega$ axis rather than the discrete axis for the DFT.

If windowing, please also refer to the last paragraph in this post to properly compensate noise and signal terms after using a window:

Blackman-Tukey Autopower equation

The FFTW page also provides a useful reference on Pruned FFTs and specifically an approach when the interest is in only the first K outputs which may be applicable to the OP's challenge:

http://www.fftw.org/pruned.html