How to find the square root of a matrix Could anyone help me to understand how the author was able to take the square root of the following matrix:

$$
\rho  = {1 \over Z}\left( {\matrix{
   {e^{ - 2\beta B} } & 0 & 0 & 0  \cr 
   0 & {\cosh 2\beta \delta } & { - e^{i\theta } \sinh 2\beta \delta } & 0  \cr 
   0 & { - e^{ - i\theta } \sinh 2\beta \delta } & {\cosh 2\beta \delta } & 0  \cr 
   0 & 0 & 0 & {e^{2\beta B} }  \cr 

 } } \right)
$$

to be in the following form: 

$$
\sqrt \rho   = {1 \over {\sqrt Z }}\left( {\matrix{
   {e^{ - \beta B} } & 0 & 0 & 0  \cr 
   0 & {\cosh \beta \delta } & { - e^{i\theta } \sinh \beta \delta } & 0  \cr 
   0 & { - e^{ - i\theta } \sinh \beta \delta } & {\cosh \beta \delta } & 0  \cr 
   0 & 0 & 0 & {e^{\beta B} }  \cr 

 } } \right)
$$

 A: You already got the simple math answer from @Doug M 's comment. But, beyond that, you are meant to recognize the familiar physics structure of it. 
That is, the evident factorization to two commuting matrices, one of them diagonal,
$$
    \rho  = {1 \over Z}\left( {\matrix{
       {e^{ - 2\beta B} } & 0 & 0 & 0  \cr 
       0 & {\cosh 2\beta \delta } & { - e^{i\theta } \sinh 2\beta \delta } & 0  \cr 
       0 & { - e^{ - i\theta } \sinh 2\beta \delta } & {\cosh 2\beta \delta } & 0  \cr 
       0 & 0 & 0 & {e^{2\beta B} }          } } \right)= {1 \over Z}\left( {\matrix{
       {e^{ - 2\beta B} } & 0 & 0 & 0  \cr 
       0 &  1 &  0 & 0  \cr 
       0 &  0&  1 & 0  \cr 
       0 & 0 & 0 & {e^{2\beta B} }          } } \right)  \left( {\matrix{
     1 & 0 & 0 & 0  \cr 
       0 & {\cosh 2\beta \delta } & { - e^{i\theta } \sinh 2\beta \delta } & 0  \cr 
       0 & { - e^{ - i\theta } \sinh 2\beta \delta } & {\cosh 2\beta \delta } & 0  \cr 
       0 & 0 & 0 &  1            } } \right).
$$
The first diagonal factor is trivial to take the square root of, diag$(e^{-\beta B}, 1,1,e^{\beta B})/\sqrt{Z}$.
So all you need to do is recognize the middle (unimodular) 2×2 block as a phased- equivalent hyperbolic rotation,
$$
M=  \begin{pmatrix}
               {\cosh 2\beta \delta } & { - e^{i\theta } \sinh 2\beta \delta }  \\
 { - e^{ - i\theta } \sinh 2\beta \delta } & {\cosh 2\beta \delta }  \end{pmatrix}=  
   \begin{pmatrix}
               1&0\\
  0 &  -e^{-i\theta}\end{pmatrix}  \begin{pmatrix}
               {\cosh 2\beta \delta } &   \sinh 2\beta \delta   \\
 \sinh 2\beta \delta  & {\cosh 2\beta \delta }  \end{pmatrix} \begin{pmatrix}
                1 &  0  \\
  0 &    - e^{ i\theta }\end{pmatrix}. 
$$
But you know the hyperbolic rotation matrix angles (rapidity) add upon composition, as you may verify here through the standard double-angle hyperbolic function identities, so the middle hyperbolic rotation resolves to just
$$
  \begin{pmatrix}
               {\cosh \beta \delta } &   \sinh \beta \delta   \\
 \sinh \beta \delta  & {\cosh \beta \delta }  \end{pmatrix} ^2,
$$
and, of course, you may insert the two phases in-between these two factor matrices to preserve the equivalence in the product. 
This is a standard calculation you  are expected to be repeating again and again, by now in your head, as you go on,  ever mindful of the interlocking pieces...
