I refer to the following structure of a phaser (from this article):
Let's call the depth factor $a$, and the feedback factor $f$ ($|f|<1$).
I assume that there are $2N$ first order allpass filters, each with transfer function
$$A_0(z)=\frac{z^{-1}-b}{1-bz^{-1}},\quad |b|<1\tag{1}$$
The allpass frequency response is obtained by evaluating (1) on the unit circle $z=e^{j\omega}$:
$$A_0(e^{j\omega})=\frac{e^{-j\omega}-b}{1-be^{-j\omega}}=
\frac{e^{-j\omega}-2b+b^2e^{j\omega}}{|1-be^{-j\omega}|^2}=\\
=\frac{(1+b^2)\cos(\omega)-2b+j(b^2-1)\sin(\omega)}{|1-be^{-j\omega}|^2}\tag{2}$$
And since $A_0(z)$ is an allpass function
$$A_0(e^{j\omega})=e^{j\phi(\omega)}\tag{3}$$
From (2) and (3) the phase $\phi(\omega)$ of the allpass satisfies
$$\tan(\phi(\omega))=\frac{(b^2-1)\sin(\omega)}{(1+b^2)\cos(\omega)-2b}\tag{4}$$
If you cascade $2N$ of these allpass filters, the frequency response of this cascade is
$$[A_0(e^{j\omega})]^{2N}=e^{j2N\phi(\omega)}$$
Let us first consider a system without feedback, i.e. $f=0$. The total frequency response of such a system is
$$H(e^{j\omega})=1+ae^{j2N\phi(\omega)}=
1+a\cos(2N\phi(\omega))+ja\sin(2N\phi(\omega))$$
and its squared magnitude is
$$|H(e^{j\omega})|^2=1+2a\cos(2N\phi(\omega))+a^2\tag{5}$$
From (5) it is obvious that the locations of the maxima of $|H(e^{j\omega})|^2$ are independent of the depth paramter $a$, and that they only depend on $N$ and on the allpass phase response $\phi(\omega)$ (i.e. on the paramter $b$). The maxima occur where the argument of the cosine in (5) equals $2\pi k$, $k=0,1,\ldots,N$, which is equivalent to
$$\phi(\omega_k)=\frac{k\pi}{N},\quad k=0,1,\ldots,N$$
where $\omega_k$ are the peak frequencies. With (4) this condition becomes
$$\frac{(b^2-1)\sin(\omega_k)}{(1+b^2)\cos(\omega_k)-2b}=\tan\left(\frac{k\pi}{N}\right),\quad k=0,1,\ldots,N\tag{6}$$
Equation (6) can be used to determine the peak frequencies $\omega_k$. Note that this needs to be done numerically since this equation cannot be solved for $\omega_k$ analytically. There are just some special indices $k$ for which an explicit solution is possible:
$$\begin{array}{rcl}\omega_0&=&0\\
\omega_{N/2}&=&\arccos\left(\frac{2b}{1+b^2}\right)\quad\textrm{if } N \textrm{ is even}\\
\omega_N&=&\pi
\end{array}$$
The following figure is an example for $N=4$, $a=0.7$, and $b=0.6$ which shows $|H(e^{j\omega})|$
and the locations of the maxima (shown by 'x') which I've computed using equation (4).
If there is feedback ($f\neq 0$) then the frequency response looks more complicated:
$$H(e^{j\omega})=\frac{1+(a-f)e^{j2N\phi(\omega)}}{1-fe^{j2N\phi(\omega)}}$$
and its squared magnitude is
$$|H(e^{j\omega})|^2=\frac{1+2(a-f)\cos(2N\phi(\omega))+(a-f)^2}{1-2f\cos(2N\phi(\omega))+f^2}$$
Luckily it turns out that feedback doesn't change the location of the peaks, it just makes them 'peakier'. In the following figure the system the frequency response shown above is compared to the response with feedback ($f=0.5$):
