Do all complex operators $T: \mathbb{C}^n \rightarrow \mathbb{C}^n$ admit a square root $S$ such that $S \circ S = T$?
It is clear that all normal operators admit this: Write the matrix $T$ in the eigenbasis (that is guaranteed from spectral decomposition of normal operators), which gives us a diagonal representation, with $t_i$ as the diagonal elements. That is, $T = \text{diag}(t_1, t_2, \dots t_n)$.
Let $S$ be the diagonal matrix (in the eigenbasis of $T$) with entries $\sqrt t_i$. is always possible since $\mathbb{C}$ is algebraically closed. That is, $S = \text{diag}(\sqrt t_1, \dots, \sqrt t_n).$ Clearly, $S \circ S = T$.
I assumed the existence of an eigenbasis of $T$ for this. What happens in general, for all operators? What about other fields? More general an answer, the better :)
For the case $T$ singular, see the Marc van Leeuwen's post in
$A = B^2$ for which matrix $A$?
In the sequel, we assume that $K$ is a perfect field.
When $T\in M_n(K)$ is invertible, if you want explicitly (that is, not an approximation) a square root, then practically the Jordan form is useless.
EDIT 1. There is a method, with polynomial complexity, that reduces the problem to the case where $T$ is semi simple. Yet, to obtain an explicit square root in this last case, there remains a gap because we must work in a field that contains the square roots of the eigenvalues of $T$.
Calculate the Jordan Chevalley decomposition: $T=D+N$ where $D$ is semi simple, $N$ is nilpotent, $DN=ND$ and $D,N$ are polynomials of $T$ of degree $<n$ with coefficients in $K$. There exists a pseudo-Newton algorithm to do that; cf. (in french)
https://www.math.u-bordeaux.fr/~jesterle/Jordan-Chevalley.pdf
Clearly, $D$ is invertible and $D^{-1}N$ is nilpotent; now $T=D(I+D^{-1}N)$ and, (formally and rigorously)
$T^{1/2}=D^{1/2}(I+1/2(D^{-1}N)-1/8(D^{-1}N)^2+\cdots+(?)(D^{-1}N)^{n-1})$.
Note that the second term of RHS is a matrix $M\in K[T]$. We want only square roots that are polynomials in $T$; yet the coefficients are no more in $K$ but in an algebraic extension $L$ of $K$ (cf. below).
Note that $D=Pdiag(\lambda_1I_{n_1},\cdots,\lambda_pI_{n_p})P^{-1}$ over $\overline{K}$, the algebraic closure of $K$ (the $(\lambda_i)$ being distinct). Then there are $2^p$ values for $D^{1/2}\in L[T]$: $D^{1/2}=Pdiag(\sqrt{\lambda_1}I_{n_1},\cdots,\sqrt{\lambda_p}I_{n_p})P^{-1}$.
Let $q$ be the minimal polynomial of $D$ (of degree $p$). Then a solution as above for $D^{1/2}$ can be written $r(D)$ where $r(x)=a_0+\cdots+a_{p-1}x^{p-1}$ satisfies $r(x)^2=x \;mod(q(x))$; it's a system of $p$ equations of degree $2$ in the $p$ unknowns $(a_i)$.
EDIT 2. Using Grober basis theory, we can eliminate $p-1$ among the $p$ unknowns; in fine, the software gives a polynomial $s\in K[x]$ of degree $2^p$ that is canceled by the last unknown (note that its coefficients are very large for $n\geq 5$). If $L$ is the decomposition field of $s$, then $r\in L[x]$ and $D^{1/2}\in M_n(L)$. Note that $r_{\epsilon}$ is the Lagrange interpolation polynomial that sends the $(\lambda_i)$ to the $(\epsilon_i\sqrt{\lambda_i})$; that explains why $L$ contains all the $(\sqrt{\lambda_i})_i$; then the Galois group of $s$ has (in the generic case where $p=n$) $n!2^n$ elements. Finally, the problem "find explicitly some square roots of $D$" is solvable (by radicals) IFF the characteristic polynomial of $T$ is solvable. For a generic $T$, that can be done only if $n\leq 4$.
$\textbf{Proposition.}$ Let $T\in GL_n(K)$ be s.t. the minimal polynomial $m_D$ of $D$ has degree $p$. Then there are $2^p$ square roots of $T$ that are polynomials in $T$; each square root is in the form $r(D)M$ where $M\in K[T]$ is calculable and $r(D)\in M_n(L)$ where $L$ is the decomposition field of a polynomial $s\in K[x]$ of degree $2^p$ which is explicitly calculable. Moreover, $s$ is solvable iff $m_D$ is solvable.
$\textbf{Remark.}$ We can do similar calculations for the $p^{th}$ roots of $T$.