A very partial answer.
Proposition 1. let $A$ be a P.D. matrix. Then there is a unique $B$ s.t. $B^2=A$ and every eigenvalue $\lambda$ of $B$ satisfies $Re(\lambda)>0$; moreover, if $A$ admits a P.D. square root, then necessarily it's $B$.
Proof. The key point is: if $U\in M_n$ is P.D., then every eigenvalue $\mu$ of $U$ satisfies $Re(\mu)>0$ (Beware, the converse is false!)
In particular, our $A$ has no $<0$ eigenvalues and, therefore, admits a unique square root $B$ s.t. every eigenvalue $\lambda$ of $B$ satisfies $Re(\lambda)>0$ (cf. Higham, functions of matrices). Thus $B$ is the only candidate that can be P.D.
Remark. That does not imply (despite Mathworld's article) that $A$ admits a P.D. square root.
EDIT 1. @Dap did a very pretty proof. I had thought of making such a recurrence, but I was sure that it would not work; which just shows that, in mathematics, you have to believe!
Moreover, using Dap's proof (mutatis mutandis), we can prove the following improvement
Proposition 2. Let $A\in M_n(\mathbb{R})$ be a P.D. matrix that satisfies $spectrum(A)\subset (0,+\infty)$. Then its principal square root (cf. Proposition 1.) is P.D.
Proof. Note that $B$ (the principal square root of $A$) has $>0$ eigenvalues and that it is triangularizable over $\mathbb{R}$ with a change of orthonormal basis.
Let $N'=\begin{pmatrix}N&x\\0&\alpha\end{pmatrix}$ (where $\alpha>0$) be the matrix $B$ after triangularization. We follow the Dap's proof.
We know that $N^2+N^{2T}>0,2\alpha^2-x^T(N+\alpha I)^T(N^2+N^{2T})^{-1}(N+\alpha I)x>0$ and we want to show that
$N+N^T>0,2\alpha-x^T(N+N^T)^{-1}x>0$.
It suffices to show that
$\Delta=x^T(N+\alpha I)^T(N^2+N^{2T})^{-1}(N+\alpha I)x- \alpha(x^T(N+N^T)^{-1}x)$ is non-negative.
We find $\Delta=N(N+N^T)N^T+2\alpha NN^T+\alpha^2(N+N^T)$ and we are done.
Remark 1. It remains to study the case when a P.D. matrix $A$ has non-real eigenvalues with positive real part.
Remark 2. Of course, Dap deserves the bounty.
EDIT 2. I just read the article ([1] by Johnson and all) cited by @Dap. Ewan will be happy; the result given by mathworld is true (what surprises me).
If $A\in M_n(\mathbb{C})$, let $F(A)=\{x^*Ax;||x||=1,x\in \mathbb{C}^n\}$ be its numerical range or its field of values. Note that $A$ is P.D. iff $F(A)\subset \{z;Re(z)>0\}$ iff $A+A^*$ is H.P.D..
[1] Theorem 7 (Kato,Masser, Neumann). If $A\in M_n(\mathbb{C})$ is s.t. $F(A)\cap (-\infty,0]=\emptyset$ (that is the case when $A$ is P.D.), then its principal square root is P.D..
[1] Corollary 8 (Johnson and all). If $A\in M_n(\mathbb{R})$ is s.t. $F(A)\cap (-\infty,0]=\emptyset$ (that is the case when $A$ is P.D. in the sense: for every $x\in \mathbb{R}^n\setminus\{0\}$, $x^TAx>0$), then its principal square root (which is a real matrix) is P.D.