Is there a closed formula for the positive-definite matrix square root?

Let $$S_{>0}$$ be the space of symmetric positive definite real $$n \times n$$ matrices.

Is there a closed-form formula for the positive square root function $$\sqrt \cdot :S_{>0} \to S_{>0}$$?

I want an explicit formula (using roots, algebraic operations, and perhaps some other natural operators) for $$\sqrt A$$ in terms of the entries of $$A$$.

So, using the fact that $$A$$ is orthogonally diagonalizable is not explicit enough for me:

We can write $$A=Q \Sigma Q^T$$, and so $$\sqrt A=Q \sqrt\Sigma Q^T$$, but there are no explicit formulas for the eigenvalues - since they are roots of a high-order polynomial (when $$n \ge 5$$).

However, a-priori the fact we do not have a formula for the eigenvalues does not rule out that the specific combination $$Q \sqrt\Sigma Q^T$$ could be expressed in a formula.

A related question is whether there is a formula for the positive factor of a matrix:

Given an $$n \times n$$ matrix $$A$$, it can be written as $$A=OP$$, where $$O$$ is orthogonal and $$P$$ is symmetric positive definite. $$P$$ is given by $$P=\sqrt{A^TA}$$.

Is there a formula for $$P$$?

(If $$A=U\Sigma V^T$$ is the SVD of $$A$$, then $$P=V\Sigma V^T$$ but again this is not really explicit).

Of course, if we had a formula for the square root, we had also a formula for $$P$$, but maybe more can be said on this particular case - in particular $$P(A)$$ behaves homogeneously with $$A$$, i.e. $$P(\lambda A)=\lambda P(A)$$, while the square root behaves differently of course.

I guess I am rather pessimistic about the existence of nice formulas, but I haven't seen an explicit discussion of these anywhere. (Except for when $$n=2$$).

$$\textbf{Proposition}.$$ Let $$A\in S_n^{++}(\mathbb{Q})$$. Then $$\chi_A$$, the characteristic polynomial of $$A$$, is not solvable by radicals IFF $$\sqrt{A}$$ is not calculable by radicals.

We consider a generic matrix $$A$$; in particular, its eigenvalues $$(\lambda_i)$$ are distinct and the Galois group of $$\chi_A$$ is the total permutation group. Note that $$\sqrt{A}$$ is a polynomial in $$A$$ and decomposes uniquely on $$\{I,\cdots,A^{n-1}\}$$. Thus, if we can calculate $$\sqrt{A}$$, then we can calculate the coefficients $$(u_i)$$ in $$\sqrt{A}=\sum_{j=0}^{n-1}u_jA^j$$.

First point of view. $$\sqrt{A}=P(A)$$ where $$P$$ is the Lagrange interpolating polynomial (of degree $$n-1$$) that sends the $$\lambda_i$$'s on the $$\sqrt{\lambda_i}$$'s. We see that the $$(u_j)$$'s are in $$\mathbb{Q}(\sqrt{\lambda_1},\cdots,\sqrt{\lambda_n})$$, an algebraic extension of $$\mathbb{Q}$$ of degree $$2^nn!$$.

In particular, if $$\chi_A$$ is solvable, then we can calculate $$\sqrt{A}$$ by radicals.

Otherwise, one may wonder if these elements are in a solvable extension of $$\mathbb{Q}$$, for example, of dimension $$2^n$$. The answer is no as we will see.

Second point of view. We write the system of $$n$$ equations in the $$n$$ unknowns $$(u_j)$$; $$(\sum_{j=0}^{n-1}u_jx^j)^2=x \;mod( \;\chi_A(x))$$.

This system reduces (in general) to solving a polynomial of degree $$2^n$$ with Galois group containing $$2^nn!$$ elements, that is not solvable when $$n\geq 5$$.

That follows is a randomly chosen example for $$n=5$$

$$A = \begin{pmatrix}21& 2& -5& -9& 1\\2& 10& -5& -6& 6\\-5& -5& 14& -6& -9\\-9& -6& -6& 24& 9\\1& 6& -9& 9& 17\end{pmatrix}>0$$.

$$\chi_A(x)=x^5-86x^4+2491x^3-27698x^2+102394x-61296$$ has $$S_5$$ as Galois group. $$u_0$$ is some root of a polynomial of degree $$32$$ with a Galois group containing $$2^55!$$ elements and $$u_1,u_2,u_3,u_4$$ are polynomials in $$u_0$$ with rational coefficients.