I'm trying to implement higher order digital filter from 2nd order design (without using some ready to use method as like MIM for the job):

H(s) = w0^2 / (s^2 + w0 s/Q + w0^2) , where w0 = 2 * pi * fc

Is there a way to manipulate this tf for higher order implementation the way that 'characteristics' of 2nd order filter are preserved (or should the design be started from the scratch and do it by the selected filter order)?
If there is a way to do it then, how it's done?

This is an easy-ish problem if you know what kind of filter you want to make (IE, have an algebraic relationship between filter parameters and pole locations). 

For example, take an Nth order butterworth filter. Recall that an Nth order butterworth has N poles located in a semicircle of radius w0 around the origin of the s plane, or poles at 

s = w0*exp(j*k*pi/N) and s = w0*exp(-j*k*pi/N), k = [1,...,N/2]. 

That means your denominator for the entire filter can be represented as the product of

s^2 + 2*w0*cos(k*pi/2) + w0^2, k = [1,...,N/2]

let Q = 1/(2*w0*cos(k*pi/2))

and you can see how your filter prototype can turn into any even order butterworth by putting N/2 of them in series and setting Q = 1/(2*w0*cos(k*pi/2) where k is the index of the filter in the chain. if you want odd orders, you can extend this process but it's a little bit trickier. 

You can also see how the trick here is that you get higher orders by putting many of your filter prototype in series and the only parameter you need to vary is Q. The catch there is you need an algebraic way of finding the appropriate Q value for the kth filter stage for the total Q to be the correct. I haven't done the math on it but I imagine it's tedious. You can start with the butterworth example and see if that is a good basis. Remember the Q of a butterworth filter is sqrt(2)/2. 

Thanks, I'll check your example case. 

I was planning to have this working with all filter types RBJ's Audio EQ Cookbok presents.

You should ask yourself, "why do I want a higher order _ filter with identical parameters, and does it justify the extra  computation?" 

For lowpass/highpass/bandpass filters I can see an argument but what's the practical difference between a 4th order parametric EQ filter and 2nd order with a higher Q? 

I'd also stick to just a lowpass because you can turn it into a highpass or bandpass really easily. And if you so choose you can take LP/HP/BP outputs and turn them into those other filter types if you wanted to.

AFAIK, higher order filter gives more options for to match the magnitude response follow closer the analog prototype response (I've faced the limits of 2nd order filter in this work).

Here's Octave script for to make 4th order CD-emphasis filter from known parameters 

T1 = 50µs; 
T2 = 15µs;
fs = 44100;

% MZT ------------------------------------
fs=44100;
samplingtime = 1/fs;
p1 = exp(-1.0/(fs*(0.000015)));
z1 = exp(-1.0/(fs*(0.000050)));
p2 = -exp(-1.0/(fs*0.00000567)) reduces the magnitude error from +0.4dB to +/- 0.06dB;
z2 = 0;
p3 = 0;
z3 = 0;
p4 = 0;
z4 = 0;

a0 = 1.0;
a1 = -p1-p2-p3-p4;
a2 = p1*p2 + p1*p3 + p1*p4 + p2*p3 + p2*p4 + p3*p4;
a3 = -p1*p2*p3 - p1*p2*p4 - p1*p3*p4 - p2*p3*p4;
a4 = p1*p2*p3*p4;

b0 = 1.0;
b1 = -z1-z2-z3-z4;
b2 = z1*z2 + z1*z3 + z1*z4 + z2*z3 + z2*z4 + z3*z4;
b3 = -z1*z2*z3 - z1*z2*z4 - z1*z3*z4 - z2*z3*z4;
b4 = z1*z2*z3*z4;

Bmzt = [b0 b1 b2 b3 b4];
Amzt = [a0 a1 a2 a3 a4];

DsMZT = tf(Amzt, Bmzt, samplingtime);
DsMZT = DsMZT/dcgain(DsMZT);

I have additional value there in p2 slot for to improve the magnitude response. There are still few p/z slots available for tweaking but I don't know what to put there programmically (giving values manually I have tried something and got the error down to -0.02-+0.03dB). I know that above implementation is not good example for what I'm after with this thread but shows that one can improve the filter by using additional poles/zeros or a/b coefficients.

I guess that intuitively makes sense, but the pole/zero locations at infinity are still going to be mapped to Nyquist no matter what filter order you use.

Why not just solve for the desired magnitude gain at Nyquist as a function of w0 and scoot your zeros back accordingly? iirc that's what Massberg did. Granted it turns into more of a shelf at Nyquist, but I'd be interested to see if humans could even perceive the difference between a Massberg style filter and a normal LPF running at an absurd sample rate. 

edit: i'm wrong about massberg. this paper has a nice summary of analog matching techniques.

Which kind of filter do you want ? IIR or FIR ? The fir2 function in MATLAB may be able to satisfy your requirement.

IIR

This really boils down to the question of what characteristics you want to preserve.  If you want to preserve all the filter characteristics, then you just want an identical filter that takes more work -- I don't think that's what you meant.

So -- what is it that you like and want to preserve about the filter behavior, and what is it you don't like and want to improve?

I suppose preserving roll-off and slope are the main targets ... kind of 1st or 2nd order filter but few additional poles/zeros (or a/b coefficients) for to tweak the filter magnitude response match closer the analog model (there are those warping and aliasing issues at low sample rates ...). 

ATM, I'm able to do this by using MIM method which takes analog model of filter as input and transforms it to selected (1st... nth) order digital filter ... quite well. These transfer functions results same filter (higher order is more accurate of course):

        0.465 z - 0.08588
 H(z):  -----------------
         z - 0.6208


        0.4602 z^8 + 1.077 z^7 + 0.6194 z^6 - 0.2768 z^5 - 0.3722 z^4 - 0.08696 z^3 + 0.01111 z^2 + 0.004812 z + 0.0002547
 H(z):  ------------------------------------------------------------------------------------------------------------------
        z^8 + 1.904 z^7 + 0.2502 z^6 - 1.39 z^5 - 0.7071 z^4 + 0.1701 z^3 + 0.1752 z^2 + 0.03267 z + 0.001412

Some good advice given herein. However, I would like to also recommend ASN FilterScript as the symbolic math interface has been specifically designed to aid designers faced with these types of issues.