anuttarasammyak said:
(X+E)^2=A+E
X=(A+E)^{1/2}-E
But I am not certain that any matrix has its "square root".
I think you can just take
\begin{align*}
X = \pm (I + A)^{1/2} - I
\end{align*}
and expand the ##(I + A)^{1/2}##. It convergences because ##A## is nilpotent.
Then you can use that the eigenvalues of ##-(I + A)^{1/2} - I## are non-zero implying it is not nilpotent. The nilpotent solution, for ##A^k=0##, is:
\begin{align*}
X & = \sum_{r=1}^{k-1} (-1)^{r-1} \dfrac{(2r)!}{(2r-1) (2^r r!)^2} A^r
\nonumber \\
& = \frac{A}{2} - \frac{A^2}{8} + \frac{A^3}{16} - \frac{5A^4}{128} + \cdots + (-1)^k \dfrac{(2k-2)!}{(2k-3) (2^{k-1} (k-1)!)^2} A^{k-1}
\end{align*}
You can check this. Say ##A^2=0##, so that ##X=a_0 I + a_1 A##, then:
\begin{align*}
(a_0 I + a_1 A)^2 + 2 (a_0 I + a_1 A) = A
\end{align*}
or
\begin{align*}
(a_0^2+2a_0) I + 2 a_1 (a_0 + 1) A = A
\end{align*}
so that
\begin{align*}
(a_0^2+2a_0) =0 , \quad 2 a_1 (a_0 + 1) = 1 .
\end{align*}
Implying ##a_0=-2 , a_1= -\frac{1}{2}## or ##a_0=0 , a_1= \frac{1}{2}##.
Say ##A^3=0##, so that ##X=a_0 I + a_1 A + a_2 A^2##, then:
\begin{align*}
(a_0 I + a_1 A + a_2 A^2)^2 + 2 (a_0 I + a_1 A + a_2 A^3) = A
\end{align*}
or
\begin{align*}
(a_0^2+2a_0) I + 2 a_1 (a_0 + 1) A + (a_1^2 + 2a_0 a_2 + 2 a_2) A^2 = A
\end{align*}
so that
\begin{align*}
(a_0^2+2a_0) =0 , \quad 2 a_1 (a_0 + 1) = 1 , \quad a_1^2 + 2a_0 a_2 + 2 a_2= 0 .
\end{align*}
Implying ##a_0=-2 , a_1= -\frac{1}{2} , a_2 = \frac{1}{8}## or ##a_0=0 , a_1= \frac{1}{2} , a_2 = -\frac{1}{8}##.
The conditions you get on ##a_0,a_1, \dots, a_{k-1}## are the same conditions you would get by plugging ##\alpha = \sum_{r=0}^\infty a_r x^r## into
\begin{align*}
\alpha^2 + 2 \alpha = x
\end{align*}
and considering terms of powers up to ##x^{k-1}## on the LHS. This is why the approach of plugging ##X=\sum_{r=0}^{k-1} a_r A^r## into ##X^2 + 2X = A## will give the same answer as expanding ##\pm (I + A)^{1/2} - I##.