Can concepts like "critical damping" or "resonant frequency" be applied to more complex systems than just a spring and damper in parallel?

The best way to understand this is through Transfer Functions and Pole-Zero Plots, these are often used in Control Theory when trying to understand/manipulate the behavior of RLC circuits, which are described by the same Differential Equations which describe Oscillators and Dampers. Consider the general form of the $text{n}^{th}$ order linear differential equation.

$$ frac{d^ny}{dt^n} + a_1frac{d^{n-1}y}{dt^{n-1}} + dotsb + a_{n-1}frac{dy}{dt} + a_ny = x(t) $$

Taking Laplace transforms on both sides yields: $$frac{Y(s)}{X(s)}=H(s)$$

Where $H(s)=frac{1}{s^n + a_1s^{n-1}+...+a_{n-1}s^1+a_n}$, now if we find the points where $H(s)$ is undefined (where the denominator is zero), those are referred to as "poles of the transfer function". The imaginary component of these poles are the damping frequencies of oscillation. So in short, you are correct, because we have a higher order differential equation, we will have more poles and as a result have more frequencies of oscillation. Note: In the second order case we have 2 frequencies, but because they are symmetric, they are $+omega$ and $-omega$, so as far as we are concerned they are the same.

This might be more than what you asked for, but Control Theory is literally an entire field devoted to Transfer Functions and understanding how its poles effect the response of a system. Moreover I only considered the case with 0 as all of its Initial Conditions, but that is quite common with these types of problems. There is a lot you can tell about a system by analyzing its Transfer Function, stability, damping, oscillating.

Here are some general links on transfer functions: https://www.cds.caltech.edu/~murray/courses/cds101/fa04/caltech/am04_ch6-3nov04.pdf https://web.mit.edu/2.14/www/Handouts/PoleZero.pdf https://lpsa.swarthmore.edu/Representations/SysRepTransformations/TF2SDE.html#:~:text=To%20find%20the%20transfer%20function,often%20called%20H(s).

this is your differential equation: $$σ + (frac{η_1}{E_1} + frac{η_2}{E_2})dot{σ} + frac{η_1η_2}{E_1E_2}ddot{σ}= (η_1+η_2)dot{ϵ} + frac{η_1η_2(E_1+E_2)}{E_1E_2}$$ Transfer this Ode to Laplace domain with:

$$sigma+ left( {frac {eta_{{1}}}{E_{{1}}}}+{frac {eta_{{2}}}{E_{{2 }}}} right) Ssigma+{frac {eta_{{1}}eta_{{2}}{S}^{2}sigma}{E_{{1} }E_{{2}}}}= left( eta_{{1}}+eta_{{2}} right) Sepsilon +{frac { eta_{{1}}eta_{{2}} left( E_{{1}}+E_{{2}} right) }{E_{{1}}E_{{2}}}} $$ solve this equation for $~sigma(S)~$ where $epsilon(S)~$ is the input

$$sigma={frac {eta_{{1}}eta_{{2}}E_{{1}}+eta_{{1}}eta_{{2}}E_{{2}}+S epsilon ,E_{{1}}E_{{2}}eta_{{1}}+Sepsilon ,E_{{1}}E_{{2}}eta_{{2 }}}{SE_{{2}}eta_{{1}}+Seta_{{2}}E_{{1}}+eta_{{1}}eta_{{2}}{S}^{2}+ E_{{1}}E_{{2}}}} $$ take the partial fraction you obtain begin{align*} &sigma={frac {-{eta_{{1}}}^{2}eta_{{2}}E_{{1}}-{eta_{{1}}}^{2}eta_{{2}}E _{{2}}+{E_{{1}}}^{2}epsilon ,E_{{2}}eta_{{1}}+{E_{{1}}}^{2} epsilon ,E_{{2}}eta_{{2}}}{ left( eta_{{2}}E_{{1}}-E_{{2}}eta_{{ 1}} right) left( Seta_{{1}}+E_{{1}} right) }}&+{frac {eta_{{1}}{ eta_{{2}}}^{2}E_{{1}}+eta_{{1}}{eta_{{2}}}^{2}E_{{2}}-{E_{{2}}}^{2} epsilon ,E_{{1}}eta_{{1}}-{E_{{2}}}^{2}epsilon ,E_{{1}}eta_{{2}} }{ left( Seta_{{2}}+E_{{2}} right) left( eta_{{2}}E_{{1}}-E_{{2} }eta_{{1}} right) }}tag 1 end{align*}

the roots of the denominators are the pols of your equation. to have oscillation for $~sigma~$ the denominator of Eq. (1) must contained at list the term $~S^2~$, thus with this equation you can't find critical frequency.