Square root of a matrix appears in massive gravity. How to solve $\sqrt{A+B}$ perturbatively? $A=\text{diag}\{\lambda_1,...,\lambda_n\}$, where $\lambda_i$ can be any number and not necessarily a small number, $\lambda_i>0$, $B$ is a positive definite symmetric matrix, and $\text{max}\{B_{ij} \}\ll \text{min}\{\lambda_i\}$.
Note that the perturbative calculation of square root of $I+B$ is very easy, where $B$ is a small matrix.
$$\sqrt{(I+B)}= I+\frac{1}{2}B-\frac{1}{8}B^2\cdots$$
In general how to calculate the square root of $A+B$ perturbatively?
This question is nontrivial because $\sqrt{A}\sqrt{B}\not=\sqrt{B}\sqrt{A}$, $\sqrt{AB}\not=\sqrt{A}\sqrt{B}$ and $\sqrt{AB}\not=\sqrt{B}\sqrt{A}$.
If we write $C=A+B-I$ and $A+B=I+C$, $C$ is not a small quantity, therefore the perturbative expansion fails.
 A: Hints:

*

*The square root function has a Taylor expansion around $a>0$
$$ \sqrt{a+b}~=~\sum_{n=0}^{\infty} \begin{pmatrix}\frac{1}{2} \cr n\end{pmatrix}a^{\frac{1}{2}-n}b^n,  \qquad |b| ~<~a. \tag{1} $$


*One may show that a possible non-commutative generalization reads
$$\begin{align}  \sqrt{A+B}~=~&\sqrt{A}+\sum_{n=1}^{\infty} \begin{pmatrix}\frac{1}{2} \cr n\end{pmatrix}\left( \left(L_B +{\rm ad} A \right)^{n-1}B \right) A^{\frac{1}{2}-n},  \cr  A~>~&0, \end{align} \tag{2}$$
and $B$ sufficiently small. Here
$$ L_B(A):=BA \quad\text{and}\quad ({\rm ad} A)B~:=~[A,B] \tag{3}$$
are left composition and adjoint action, respectively.
A: An elementary way to proceed is as follows. Let's put an explicit factor of $\epsilon$ in $B$. The problem is then to solve the following equation for the matrix $X$:
$$A + \epsilon B = X^{2}$$
We want to do this perturbatively, so we assume that $X$ can be represented as:
$$X = \sum_{n=0}^{\infty}\epsilon^{n}X^{(n)}$$
We can then write:
$$
\begin{split}
X^2 &= \sum_{n=0}^{\infty}\epsilon^{n}\sum_{k=0}^{n}X^{(k)}X^{(n-k)}= \left(X^{(0)}\right)^2  + \epsilon\left(X^{(0)}X^{(1)}+X^{(1)}X^{(0)}\right) \\&+\sum_{n=2}^{\infty}\epsilon^{n}\left(X^{(0)}X^{(n)}+X^{(n)}X^{(0)} + \sum_{k=1}^{n-1}X^{(k)}X^{(n-k)}\right)
\end{split}
$$
We then see that:
$$X^{0} = A^{\frac{1}{2}}$$
In components we can write this as:
$$X^{(0)}_{i,j} = \sqrt{\lambda_{i}}\delta_{i,j}$$
where the summation convention is modified so that contraction only happens when that would also happen absent any factors of the $\lambda$'s in an equation.
We can then simplify the expression $Y^{(n)} \equiv X^{(0)}X^{(n)} + X^{(n)}X^{(0)}$:
$$Y^{(n)}_{i,j} =  \left(\sqrt{\lambda_{i}}+\sqrt{\lambda_{j}}\right)X^{(n)}_{i,j}$$
We thus have:
$$X^{(1)}_{i,j} = \frac{B_{i,j}}{\sqrt{\lambda_{i}}+\sqrt{\lambda_{j}}}$$
And the higher order terms are then obtained recursively from:
$$X^{(n+1)}_{i,j} = -\frac{\sum_{r=1}^{n}X^{(r)}_{i,k}X^{(n+1-r)}_{k,j}}{\sqrt{\lambda_{i}}+\sqrt{\lambda_{j}}}$$
