This is a fun question. I thought a bit about it and the specific $D_4^-$ example, and here is what I got (apologies for a longish answer).
Refined count of singularities in unfoldings:
First, let me note that the conditions that Morse indices have to add up can be reformulated in terms of a refinement of the Milnor number. From the work of Eisenbud-Levine and Khimshiashvili, there is a symmetric bilinear form on the local algebra $Q$, given as follows: choose a linear form $\phi:Q\to \mathbb{R}$ such that the image of the Hessian matrix $(\frac{\partial^2 f}{\partial x_i\partial x_j})_{i,j}$ has positive image under $\phi$; then consider the symmetric bilinear form $\langle f,g\rangle:=\phi(fg)$. The rank of this form is the Milnor number for the complexification, and the signature provides refined information about the singularities which appear in unfoldings of $f$.
More precisely, consider an unfolding $f(x_1,\dots,x_n)+\sum a_ix_i$ (only adding linear terms!) and assume it has only Morse (i.e. nodal) singularities. To each of the nodal singularities, one can associate a similar bilinear form, and the direct sum of those must be isometric to the Eisenbud-Levine-Khimshiashvili form of the singularity of $f$. I learned this from the fantastic work of Kass and Wickelgren:
- J. Kass and K. Wickelgren. The class of Eisenbud-Khimshiashvili-Levine is the local $\mathbb{A}^1$-Brouwer degree. arXiv:1608.05669v1. (link to paper on arXiv)
Also check out the references in there. Their results work for varieties over general fields, but if you're only interested in working over the real numbers, there are earlier versions of refinements for Milnor numbers e.g. in the work of Wall on topological invariance of Milnor numbers mod 2:
Singularity counting in case $D_4^-$:
In the example of the singularity $D_4^-$ in the question, we have
$$
\det\left(\begin{array}{ccc}
2y&2x&0\\2x&-6y&0\\0&0&1\end{array}\right)=-24y^2.
$$
So we choose the linear form on $Q$ which maps $1,x,y$ to $0$ and $y^2$ to $1$. In the basis $x,y,\frac{y^2+1}{2},\frac{y^2-1}{2}$, the symmetric bilinear form is given by the diagonal matrix ${\rm diag}(-1,-1,-1,1)$.
Note that if $n$ is odd the quadratic form isn't an invariant of the singularity but of the equation. Multiplying the equation in the question by $-1$ produces the diagonal matrix $(1,1,1,-1)$. (This works better in even dimensions where the quadratic form is really an invariant of the singularity.)
In any case, most of the preliminary exclusion of possible singularities formulated in the question can be rephrased in terms of the quadratic form above. It should be noted that arguments with the form also exclude that the singularity bifurcates into two pairs of complex conjugate nodes -- at least two of the nodes must be defined over $\mathbb{R}$. However, this reasoning doesn't exclude the possibility of $A_2$-singularities appearing.
Space of unfoldings: Ok, now we can also look at the space of all unfoldings $u=f+\sum a_ig_i$ where $g_i$ runs through representative polynomials of the local algebra $Q$. If we fix a parameter tuple $(a_1,\dots,a_{\mu-1})$ we can compute the zeros of $\nabla u$. The vanishing loci for the components of $\nabla u$ are hypersurfaces in the space $\mathbb{R}^{\mu-1}$ of all unfoldings. The simple zeroes (i.e. where the hypersurfaces for the components of $\nabla u$ intersect transversely) will correspond to Morse singularities, the multiple zeroes to more complicated ones. The system of equations defining the locus of multiple zeroes is an analogue of the discriminant of a one-variable polynomial (but here we have several variables $a_1,\dots,a_{\mu-1}$). Nevertheless, the complement of this locus of multiple zeroes will consist of several connected components. Now, one of the possible definitions of "essentially different" for Morse unfoldings is: "lying in different connected components". (Note that this only counts essentially different unfoldings where only nodal/Morse singularities appear. The more complicated singularities that in the multiple-zero-locus of $\nabla u$ have to be investigated separately.)
So here we have a more precise formulation of the question: count the connected components of the complement of the discriminant locus in the space of all unfoldings. I don't think that this question can be answered by general theory (at this point):
determining the discriminant locus is not so easy. This is the locus where the scheme-theoretic intersection of the (singular) hypersurfaces where the components of $\nabla u$ vanish is non-reduced. Already working out the degree would require more elaborate methods, and to get explicit formulas is more complicated than formulas for the discriminant of one-variable polynomials. (The discriminant-resultant book of Gelfand-Kapranov-Zelevinsky could be very useful here.)
Determining the number of connected components of the complement seems a pretty complicated task. In real algebraic geometry, there are generally only bounds for the number of connected components of the real points of a variety in terms of information like degree, genus etc. (e.g. Harnack's bound) but the exact number can vary in families.
Addendum: An estimate for the number of connected components of the complement can be found in the answers to this MO-question. Papers that study statistics for Hilbert's 16th problem seem to say that the expected number of components for random hypersurfaces are close to the upper bound (which for degree 4 in $\mathbb{R}^3$ for the $D_4$-singularity would be 216) -- but then the discriminant isn't just any random polynomial.
- at this point I am not even sure if the number of connected components will be the same for all possible realizations of the singularity.
This is similar to other problems in real enumerative geometry: the space of solutions to the enumeration problem has several connected components and the number of solutions is only locally constant and varies between different components. All you can hope for is some signature-type information (as in the refined Milnor number discussed in the beginning). There seems to be no general theory which would answer questions concerning the number of components in the solution space, the possible number of real solutions in each of the components or the more precise information on the topology of the components of the solution space.
Two examples for the structure of the space of unfoldings are discussed below.
A simpler example $A_3$: The hypersurface is given by $x^2+y^2+z^4=0$. In this case, the local algebra is generated by $1,z,z^2$. So the unfoldings would be of the form $g(x,y,z)=x^2+y^2+z^4+a_1z+a_2z^2$. Then
$$
\nabla g=(2x,2y,4z^3+a_1+2a_2z).
$$
Therefore, the singularities always have $x=y=0$ and the $z$-coordinates of the singularities satisfy $4z^3+2a_2z+a_1=0$. The discriminant hypersurface for this degree 3 polynomial is defined by $-8a_2^3-27a_1^2=0$. This is a curve which separates the $(a_1,a_2)$-plane into two connected components. On the one component where the discriminant is positive, there are three distinct real Morse singularities, and on the component where the discriminant is negative, there is one real Morse singularity (and the other two are complex conjugates). In this case we can confidently say that there are two types of nodal/Morse unfoldings, determined by their numbers of real Morse singularities. The types of more complicated singularities that appear on the discriminant hypersurface have to be investigated separately.
Space of unfoldings for example $D_4^-$: the hypersurface in question is $x^2y-y^3+z^2$. The local algebra is generated by $1,x,y,y^2$ and the unfoldings are of the form $u(x,y,z)=x^2y-y^3+z^2+a_1x+a_2y+a_3y^2$. Then
$$
\nabla u=(2xy+a_1,x^2-3y^2+a_2+2a_3y,2z).
$$
All the singularities will have $z=0$. There is a case distinction:
Now we can write out the discriminant for the degree 4 polynomial which determines the locus where we have multiple zeroes. The hypersurface has equation
$$
a_1^2\left(108a_1^4-72a_1^2a_2^2-108a_1^2a_2a_3^2-12a_2^4-27a_1^2a_3^4-4a_2^3a_3^2\right)=0.
$$
The component $a_1=0$ is sort of artificial: outside the line $a_1=a_2=0$, there are always four distinct zeroes of $\nabla u$, given by $(0,\frac{a_3}{3}\pm\frac{\sqrt{a_3^2+3a_2}}{3})$ and $(\pm\sqrt{-a_2},0)$. On the line $a_1=a_2=0$, there are multiple singularities, these are probably still the original $D_4^-$ singularities.
The other component is the interesting one, a hypersurface of degree $6$. I don't quite know how many components the complement of this hypersurface has. Some plots in the $(a_1,a_3)$-plane for various fixed values of $a_2$ suggest that it's more than three. So I would expect that the classification of the different types is a bit more complicated than suggested in the question. Also, for now I haven't looked at the singularities that appear on for parameters $(a_1,a_2,a_3)$ on the discriminant hypersurface. I would expect that these are non-Morse singularities, maybe with the $A_2$ suggested in the question.
Well, that's all for now, and it already suggests that this is going to be more complicated for other singularities (like $E_8$-type). Let me know if this goes in the direction you were interested in, or if this doesn't answer the question at all. I apologize again that the length of this answer has gotten out of hand...
