Rather than the frequency domain, let's look at this in the time domain and particularly, the characteristic equation associated with a linear homogeneous 2nd order differential equation for some system:
\$r^2 + 2 \zeta \omega_n r + \omega^2_n = 0\$.
If the roots of the characteristic equation are real (which is the case if \$\zeta \ge 1\$), the general solution is the sum of real exponentials:
\$Ae^{\sigma_1 t} + Be^{\sigma_2t} \$
where
\$\sigma_1 = -\zeta \omega_n + \sqrt{(\zeta ^2 - 1)\omega^2_n} \$
\$\sigma_2 = -\zeta \omega_n - \sqrt{(\zeta ^2 - 1)\omega^2_n} \$
Since these are real exponentials, there is no oscillation in these solutions.
If the roots are complex conjugates (which is the case if \$\zeta < 1\$), the general solution is the sum of complex exponentials:
\$e^{\sigma t}(Ae^{j\omega t} + Be^{-j\omega t})\$
where
\$\sigma = -\zeta \omega_n\$
\$\omega = \sqrt{(1 - \zeta ^2)\omega^2_n}\$
This solution is a sinusoid with angular frequency \$\omega\$ multiplied by a real exponential. We say the system has a "natural frequency" of \$\omega\$ for a reason that I think is obvious.
Finally, setting \$\zeta = 0\$ (an undamped system) , this solution becomes:
\$Ae^{j\omega_n t} + Be^{-j\omega_n t}\$
which is just a sinusoid of angular frequency \$\omega_n\$.
In summary, a system may or may not have an associated natural frequency. Only systems with \$\zeta < 1\$ have a natural frequency \$\omega\$ and only in the case that \$\zeta = 0\$ will the natural frequency \$\omega = \omega_n\$, the undamped natural frequency.
The specific formula applies only for a first order RC low pass filter. This is derived from its frequency response:
$$H(j\omega)=\frac{1}{1+j\omega RC}$$
The cutoff frequency is defined as the frequency where the amplitude of \$H(j\omega)\$ is \$1\over\sqrt2\$ times the DC amplitude (approximately -3dB, half power point).
$$|H(j\omega_c)|=\frac{1}{\sqrt{1^2+\omega_c^2R^2C^2}}=\frac{1}{\sqrt{2}}\cdot|H(j0)|=\frac{1}{\sqrt{2}}$$
Solve it for \$\omega_c\$ (cutoff angular frequency), you'll get \$1\over RC\$. Divide that by \$2\pi\$ and you get the cutoff frequency \$f_c\$.
If you know the frequency response of your filter, you can apply this method (given that the cutoff frequency is defined as above). Obviously, for high-pass filters for example, you calculate with the value for \$\omega\to \infty\$ as opposed to the DC value (always the maximum of the amplitude response, relative to which there is a 3dB decrease in amplitude at the cutoff frequency.)
Best Answer
This type of passive circuit can be easily solved and expressed in a so-called low-entropy format using the fast analytical circuits techniques or FACTs. Without writing a single line of algebra, you can "inspect" the circuit and determine the transfer function. In this approach, you determine the natural time constants of the circuit by reducing the stimulus \$V_{in}\$ to 0 V. When you do that, the left terminal of \$R_1\$ is grounded. In this configuration, remove the capacitors and "look" at the resistance from their terminals. The obtained resistance multiplied by the capacitance forms the time constant \$\tau\$ we need. Here, we have two energy-storing elements (with independent state variables) so this is a 2nd-order circuit obeying the following expression for the denominator \$D(s)\$:
\$D(s)=1+s(\tau_1+\tau_2)+s^2\tau_1\tau_{12}\$
We start with \$s=0\$ for which you open all caps. The transfer function is simply:
\$H_0=\frac{R_3}{R_3+R_1+R_2}\$
If you now apply the technique consisting of "looking" at the resistance offered by the capacitors terminals while \$V_{in}\$ is 0 V, you should find:
\$\tau_1=C_1(R_1||(R_2+R_3))\$
\$\tau_2=C_2(R_3||(R_1+R_2))\$
\$b_1=\tau_1+\tau_2=C_1(R_1||(R_2+R_3))+C_2(R_3||(R_1+R_2))\$
Then, consider shorting \$C_1\$ while you look at the resistance offered by \$C_2\$ terminals in this mode. You have
\$\tau_{12}=C_2(R_2||R_3)\$
\$b_2=\tau_1\tau_{12}=C_1(R_1||(R_2+R_3))C_2(R_2||R_3)\$
Assembling these expressions, we have the complete transfer function as there is no zero in this network.
\$H(s)=H_0\frac{1}{1+s(C_1(R_1||(R_2+R_3))+C_2(R_3||(R_1+R_2)))+s^2(C_1(R_1||(R_2+R_3))C_2(R_2||R_3))}\$
This is a second-order polynomial form obeying:
\$H(s)=H_0\frac{1}{1+\frac{s}{\omega_0Q}+(\frac{s}{\omega_0})^2}\$
in which \$Q=\frac{\sqrt{b_2}}{b_1}\$ and \$\omega_0=\frac{1}{\sqrt{b_2}}\$
If \$Q\$ is sufficiently low (low-\$Q\$ approximation) you can replace the second-order polynomial form by two cascaded poles. All appears in the below picture:
If you look at the raw expression transfer function \$H_{ref}(s)\$ (using Thévenin), it perfectly matches the low-entropy version. The difference is that you now have a well-ordered transfer function letting you calculate the values for all components depending on how you want to tune this filter. What truly matters is the low-entropy well-ordered form which tells you what terms contribute gains (attenuation), poles and zeros. Without this arrangement, there is no way you can design your circuit to meet a certain goal. To my opinion, the FACTs are unbeatable to obtain these results in one clean shot. Furthermore, as you can see, I have not written a single line of algebra. All I did was inspecting the network (through small individual sketches if necessary).
You can discover FACTs further here
http://cbasso.pagesperso-orange.fr/Downloads/PPTs/Chris%20Basso%20APEC%20seminar%202016.pdf
and also through examples published in the introductory book
http://cbasso.pagesperso-orange.fr/Downloads/Book/List%20of%20FACTs%20examples.pdf