Why do we divide by $n-1$ when calculating sample correlation? I understand the rationale for dividing by $n-1$ when calculating the sample variance, i.e. that if we divide by $n$ we will have an estimate of population variance that is biased to be too low.
Buglear (2013, p. 57) states about the Pearson correlation:

We divide by $n - 1$ for the same reason as we do it when calculating
sample standard deviations – it gives us a better
estimator of the population equivalent.
[Buglear, J. (2013). Practical Statistics: A Handbook for Business Projects. Kogan Page Publishers]

However, I don't understand why this also applies to correlations. Why is it the case that dividing by $n$ would underestimate the population correlation coefficient?
 A: We do not need the Bessel correction "-1" to $n$ when we compute correlation, so I think the citated piece is wrong. Let me start by noticing that most of time we compute and use empirical $r$, or the $r$ of the sample, for both describing the sample (the statistic) and the population (the parameter estimate). This is different from variance and covariance coefficients where, typically, we introduce the Bessel correction to distinguish between the statistic and the estimate.
So, consider empirical $r$. It is the cosine similarity of the centered variables ($X$ and $Y$ both were centered): $r= \frac{\sum{X_cY_c}}{\sqrt{\sum X_c^2\sum Y_c^2}}$. Notice that this formula doesn't contain neither $n$ nor $n-1$ at all, we need not to know sample size to obtain $r$.
On the other hand, that same $r$ is also the covariance of the z-standardized variables ($X$ and $Y$ both were centered and then divided by their respective standard deviations $\sigma_x$ and $\sigma_y$): $r= \frac{\sum{X_zY_z}}{n-1}$. I suppose that in your question you are speaking of this formula. That Bessel correction in the denominator, which is called in the formula of covariance to unbias the estimate, - in this specific formula to compute $r$ paradoxically serves to "undo" the unbiasing correction. Indeed, recall that $\sigma_x^2$ and $\sigma_y^2$ had been computed using denominator $n-1$, the Bessel correction. If in the latter formula of $r$ you unwind $X_z$ and $Y_z$, showing how they were computed out of $X_c$ and $Y_c$ using those "n-1"-based standard deviations you'll find out that all "n-1" terms cancel each other from the formula, and you stay in the end with  the above cosine formula! The "n-1" in the "covariance formula" of $r$ was needed simply to take off that older "n-1" used.
If we prefer to compute those $\sigma_x^2$ and $\sigma_y^2$ based on denominator $n$ (instead of $n-1$) the formula for yet the same correlation value will be $r= \frac{\sum{X_zY_z}}{n}$. Here $n$ serves to take off that older "n" used, analogously.
So, we needed $n-1$ in the denominator to cancel out the same denominator in the formulas of variances. Or needed $n$ for the same reason in case the variances were computed as biased estimates. Empirical $r$ is itself not based on the information of the sample size.
As for a quest of better population estimate of $\rho$ than the empirical $r$, there we do need corrections, but there exist various approaches and a lot of different alternative formulas, and they use different corrections, usually not $n-1$ one.
