Why can adding delay improve the phase fit in fitting complex transfer functions?
Even though an IIR filter of a given order can in principle approximate any given phase response, you will get a better fit if you choose the total delay such that it "matches" the chosen filter order. The problem is that there is no formula to help you choose such an ideal delay; you simply have to find out by trial and error. As you've already mentioned, filter stability is also determined by the choice of the total delay. You have to choose the total delay such that the given desired phase response can be reasonably approximated by a stable IIR filter of the chosen order. Also here there are no formulas for how to choose the total delay. With a bit of trial and error you will gain some intuition for it.

In sum, if $\phi(\omega)$ is your desired phase response, you should specify the phase to be approximated as

$$\tilde{\phi}(\omega)=\phi(\omega)-\tau\omega\tag{1}$$

where the delay parameter $\tau$ (which can be positive or negative) is an additional design variable.

The reason why an additional constant group delay term can make such a difference is that each pole and each zero contribute to the total group delay. The closer to the unit circle they are the bigger is their contribution. So for a given filter order, if the total delay is too small, some poles will move outside the unit circle causing a negative contribution to the group delay. This can be viewed as a non-causal approximation of the desired response. However, since we can only implement causal IIR filters, the corresponding filter is unstable. Increasing the total delay will move all poles inside the unit circle, hence the designed causal filter will be stable. Then the approximation error can be minimized by choosing the ideal total delay such that it best matches the chosen filter order. Note that knowing the filter order is not sufficient for choosing an optimal total delay, because the latter will depend on the exact positions of all poles and zeros.

The following example hopefully illustrates these points. The desired frequency response is low pass with a constant group delay in the pass band. The pass band edge is $\omega_p=0.2\pi$, and the stop band edge is $\omega_s=0.25\pi$. So the desired frequency response is

$$D(\omega)=\begin{cases}e^{-j\omega\tau},&amp;0\le\omega\le\omega_p\\
0,&amp;\omega&gt;\omega_s\end{cases}$$

where $\tau$ is a design parameter. I chose a numerator order $M=16$ (i.e., $16$ zeros), and a denominator order $N=6$ (i.e., $6$ poles away from the origin of the complex plane). It turns out that for $\tau&lt;17$ samples, the designed filter is unstable, i.e., there are poles outside the unit circle. For $\tau=17$ samples, we get a stable filter, but with two poles very close to the unit circle (pole radius $r=0.99576$). This filter is very hard to implement, and its frequency response has a large peak close to the pass band edge (see blue curve in the figure below). Choosing a slightly larger desired pass band group delay of $\tau=19$ samples decreases the pole radius to $r=0.98725$, and hence removes the frequency response peak (see green curve below). Finally, further increasing the desired group delay will further decrease the maximum pole radius, but the approximation error will increase. For $\tau=30$ samples the maximum pole radius is $r=0.97253$, but the approximation error is considerably larger (red curve below).



Also have a look at these answers to related questions: link 1, link 2, link 3, link 4.

You might also be interested in this article, which presents an algorithm that automatically chooses an optimal delay.