Frequency response of a digital state variable (Chamberlin) filter
The main reason you're seeing these results is because your transfer function is in the $z$-domain. The unique frequency response is limited to $[-\pi, \pi]$ and is periodic. In your case, you are using invalid values of $f$ to define where the $z$-domain spectrum exists. This is why you get those spikes. Using your frequency vector, I've replicated your result below

What you're generating is somewhat actually correct, but you're looking at it in the wrong way.
If you choose the range to be between $[0.1, \pi]$, you'll get the following frequency response

It's a little cleaner and you can play around with how you define your function and frequency axis to make it look more ideal. If you look at the data marker, the 3-dB point is still at that pesky $2 \times10^{-5}$ Hz. This is still correct! It's just that this is a &quot;discrete&quot; frequency.
In order to yield the continuous-time frequency $f$, you must convert from the discrete-domain frequency $f'$ by using the formula
$$f = f'fs  \space \text{Hz}$$
Using the 3-dB cutoff frequency we get
$$f = (2 \times 10^{-5})(500 \times 10^3 )= 10 \space \text{Hz}$$
Which is exactly what you expect. There's a related answer I have regarding continuous vs discrete frequencies along with MATLAB-specific details here.
EDIT: Using the filter
To show that the filter works, we generate and filter the signal
$$x(t) = \cos(2\pi(5)t) + \cos(2\pi(10)t) + \cos(2\pi(50)t)$$
Which are sinusoids are $5 \text{ Hz}$, $10 \text{ Hz}$, and $50 \text{ Hz}$.
I used MATLAB's filter() function which takes the transfer function's coefficients and generates a difference equation to perform the filtering. Identifying the coefficients is trivial given that your function is already in rational form. Python should have an equivalent.
Below is the signal $x(t)$ before and after filtering

You can see that the $5 \text{ Hz}$ component survives, the $10 \text{ Hz}$ component is partially attenuated since it's at the 3-dB cutoff, and the $50 \text{ Hz}$ tone is suppressed. Below is the MATLAB code I used to generate these results.
%% Sampling and constants

fc = 10;
fs = 500e3;
fn = fs/2;

q = 1/sqrt(2);
Q = 1/q;
Fc = 2*sin(pi*fc/fs);

%% Manually define the transfer function. Uncomment to generate and manually plot the frequency response.
% f = logspace(-10, pi, 1e5);
% w = 2.*pi.*f;
% z = exp(-1i.*w);
% 
% freqResponse = (Fc.^2.*z)./(1 - z.*(2 - Q.*Fc - Fc.^2) + z.^2.*(1 - Q.*Fc));
% 
% figure;
% semilogx(f, 20*log10(abs(freqResponse)));
% xlabel(&quot;Normalized Frequency (Hz/sample)&quot;);
% ylabel(&quot;Magnitude (dB)&quot;);
% axis tight;
% ylim([-120 10]);

%% Using built-in function filter()

b = [0 Fc.^2];
a = [1 -(2 - Q.*Fc - Fc.^2) (1 - Q.*Fc)];
[h, w] = freqz(b, a, 1e5);

figure;
semilogx(w./(2*pi), 20*log10(abs(h)));
axis tight;
ylim([-120 10]);

%% Use the filter to process a signal

t = 0:1/fs:2;

x = cos(2*pi*(5).*t) + cos(2*pi*(10).*t) + cos(2*pi*(50).*t);

nfft = 10*numel(x);
f = fs.*(-nfft/2:nfft/2-1)./nfft;

figure;
subplot(2, 1, 1);
plot(f, abs(fftshift(fft(x, nfft)./nfft)).^2);
xlim([-100 100])
xlabel(&quot;Frequency (Hz)&quot;);
ylabel(&quot;Magnitude&quot;);
title(&quot;Original Signal&quot;);

subplot(2, 1, 2);
plot(f, abs(fftshift(fft(filter(b, a, x), nfft)./nfft)).^2);
xlim([-100 100])
xlabel(&quot;Frequency (Hz)&quot;);
ylabel(&quot;Magnitude&quot;);
title(&quot;Filtered Signal&quot;);