Classical spectrum analyzer model
I believe the issue is the OP is simulating a spectrum analyzer with a &quot;Zero-IF&quot; architecture so is seeing the image issues as the sweep signal approach the input signal (since both are real signals). The product is down-converted to baseband with a real LO and passed through a low-pass Gaussian filter, so the resulting signal as it approach the DC center as a very low frequency real sine-wave has increasingly fewer cycles (and ultimately partial cycles) which then contributes to high variability in the power estimate. By using a complex LO instead, the envelope of the resulting signal is constant eliminating that variability (or by observing the two real components of the baseband complex signal $e^{j\omega t} = \cos(\omega t) + j\sin(\omega t)$ we see that as one partial cycle is decreasing the other is increasing to maintain the constant power at the detector).
I confirmed this was the case by changing the LO to a complex signal in the OP's code which would then match the traditional architecture of an actual Zero-IF receiver:
product(:,i) = (signal.*exp(-j*2*pi*freqs(i)*time))';

and use a complex conjugate product for the power (which also could be done using abs(x).^2)
values(i) = 10*log10(sum(product((order+1):end,i).*(conj(product((order+1):end,i)))))-power;

Resulting in:

The other option (that I would NOT recommend for simulation! - unless you need to prove to yourself that it is identical) is to move the filter from being a low-pass filter to a band-pass filter and offset the LO such that an IF frequency passes through the filter (matching traditional spectrum analyzer architectures).
For simulation purposes I would always simulate the equivalent baseband signal using complex signals rather than actual carriers which would be significantly more time consuming to process.
Further, this process can be done in one operation using an FFT which would be significantly faster and for a stationary input such as this produce an identical result. The Gaussian filter could be implemented by multiplying the time domain waveform by the Gaussian window prior to taking the FFT. For non-stationary signals further additions could be made to properly emulate the effect of the sweep time.
As far as the bandwidth, the Gaussian window used in the OP's code with
$\alpha = 10$ has an equivalent noise BW of 5.64 bins as computed from:
$$\text{ENBW} = N\frac{\sum (w[n]^2)}{(\sum w[n])^2} \tag{1} \label{1}$$
That said for a sampling rate of 10e6, and a target resolution BW of 10KHz, this would require a time capture of
$$N = 5.64 \times (10e6/10e3) = 5640 \text{ samples}$$
This would then also match the duration for the Gaussian filter for a 10KHz resolution BW in the OP's code by selecting the coefficients directly from the window function of that length, as in:
b =  gaussian(5640, 10)

I use equivalent noise bandwidth (ENBW) as resolution BW instead of the 3 dB BW, as it can be directly computed from the window as in \ref{1} and accurately indicates the same power that would be in a brick-wall filter of that bandwidth for a white noise source. The ENBW is slightly higher than the 3 dB bandwidth (on the order of 5% but depends on the window).
As recently reiterated by fred harris at the DSP Online Conference (https://www.dsponlineconference.com/), the Gaussian window is NOT the best choice. As he detailed, the reason Gaussian is often selected for spectral analysis is that it is believed to have the minimum time-bandwidth product, so minimum localization in time and frequency. However this is only true for a true-Gaussian which requires an infinite time support and only in that case when $t \rightarrow \pm\infty$ results in the minimum achievable time-bandwidth product BT of $BT = 0.5$. When we truncate the Gaussian response for the time-limited window (as required for practical implementation) the truncated Gaussian $BT&gt;0.5$ and is no longer the minimum BT window function available for the case of finite time duration.
A window that has minimum time BW product under finite time support is the DPSS (Digital Prolate-Spheroid Sequence) window and very close to that is the Kaiser Window, both of which are available in MATLAB, Octave and Python scipy.signal. These windows also have decreasing stop-band attenuation which is also attractive for many applications. For more info on this see: https://ccrma.stanford.edu/~jos/sasp/Slepian_DPSS_Window.html , https://www.dsprelated.com/freebooks/sasp/Kaiser_DPSS_Windows_Compared.html#fig:dpsstest and https://ieeexplore.ieee.org/document/543677