Explanation of Spectral Line Splitting when AR process is overmodeled
One way of doing power spectrum estimation of random processes actually requires finding a set of ARMA model parameters that best approximates the given observed random process $x[n]$.

In your case, you selected an all-pole signal model with $p$ poles as the roots of the denominator polynomial $A_p(z)$ with coefficients $a_p[k]$, and a $b(0)$ related with power scaling of the process. Each pole represents a peak in the spectrum for $0 \leq \omega \leq \pi$.

So, PSD estimation is reduced to finding (estimating) all-pole model parameters $a_p[k]$ and $b[0]$, from the given data by first estimating the auto-correlation sequence values $r_x[k]$ for $k=0,1,2,...,p$  , and then solving a set of linear equations (Prony or Yule-Walker etc) to find those $a_p[k]$ and $\epsilon_p$.

Note that there are a number of approaches in estimating the ACS $r_x[k]$ of $x[n]$, such as Auto-Correlation, Covariance, Burg's, Modified Covariance, etc...

Certain of these techniques (such as the autocorrelation and Burg's method) produce an artifact known as spectral splitting when they are forced to fit a high order model to an actually low order data; aka a model-data mismatch problem. There is not much you can do to prevent it. I don't know an analytic derivation of the relation between the chosen model order $p$ and the number of splitted peaks you will observe. That's also dependent on the observed noise instance being processed. So that number of peaks is also random. 

However, it's stated that modified covariance method of ACS estimation can overcome the spectral splitting problem. 

Furthermore, if you continue reading, you will see in the section named Selecting the Model Order (p. 445, Eq.8.23) that by trying to make a good guess of model order you may reduce the chances of generating that line splitting problem.