Number of microphone limitation for Spherical Harmonics based DOA estimation
First, I do not find any information on what is the order of the sound field. Can anyone elaborate please? 


$Y_{lm}(\Omega)$ is the spherical harmonics with mode $m$, order $l$ and frequency $\Omega=(\theta,\phi)$ (dictating angles of arrival, ex: azimuth, spherical, elevation, etc.). According to [1], we have the following formula:

$$Y_{lm}(\Omega) = \sqrt{\frac{2l + 1}{ 4\pi} \frac{(l-m)!}{(l+m)!} }P_l^m(\cos \theta)e^{j m \phi}$$
where $P_l^m(\cos \theta)$ is the associated Legendre functions of degree (mode) $m$ and order $l$. The above serves as a basis of spherical harmonics (similar to Fourier expansion).


  Second, does that mean that I can do very little with 4 microphones? Is there a lower limit on the number of microphones for a DOA estimation and for a number of active sources estimation?


The classical rule of thumb is that given $q$ sources, you need at least $q+1$ antennas to resolve uniquely your sources. However, this is really algorithmic dependent. For example, the paper in [2] (which formulate a 2D AoA DoA estimation) describes shift invariant methods with different lower bounds. So, in principle, the more antennas you have the better resolvability and lower MSE you will have, since the Cramer-Rao bound of DoA estimation is proportional to $\frac{1}{N}$ where $N$ is the number of antennas.


  What other, none neural network, methods can I use to estimate the number of active sources and DOA with (no necessarily by the same model)?


You can use subspace methods like MUSIC, or co-variance fitting. Some generic methods aim at directly optimizing the Maximum likelihood problem like Expectation-Maximization or IQML for specific antenna configurations. A 2D IQML approach was derived in [3].

[1] Rafaely, Boaz. Fundamentals of spherical array processing. Vol. 8. Berlin: Springer, 2015.

[2] Bazzi, Ahmad, Dirk TM Slock, and Lisa Meilhac. "Single snapshot joint estimation of angles and times of arrival: A 2D Matrix Pencil approach." 2016 IEEE International Conference on Communications (ICC). IEEE, 2016.

[3] Bazzi, Ahmad, Dirk TM Slock, and Lisa Meilhac. "Efficient maximum likelihood joint estimation of angles and times of arrival of multiple paths." 2015 IEEE Globecom Workshops (GC Wkshps). IEEE, 2015.