Thank you again for the reply. I have tried to ask clarifications about the correctness of the answer provided in the second link. Unfortunately, the way I managed things were not suitable with the site's policies and so I was forced to ask a new question (the one asked here). In any case I would like to hear Carlo's derivation of that formula in order to be sure that we are not missing anything at all here and, if needed, to adjust the reply appeared in the the second link. By the way Igor: essentially you are telling me that the answer in the second link is wrong, isn't it?

@Igor: you are right about missing information in the present post and I apologize for the inconvenience. The problem I need to solve is (hopefully) well stated in the first link and is the following: I need an expansion of $(A+n^{-1/2}B)^{1/2}$ in terms of solely $A$ and $B$. The Lyapunov equation just came out as a device to solve the problem. As I am not a mathematician and because Carlo's reply (the one in the second link) meet my requirements, I was trying to figure out if there exists any relation between the solution to $B=A^{1/2}D+DA^{1/2}$ and Carlo's answer.

@Carlo: as I expected the integral does not have an explicit expression. So I was wondering how it occurred that in the second link you came up with the answer $(A+B)^{1/2}=A^{1/2}(1+A^{-1}B/2)$? I mean which is the path that you follow to derive such a result? Because mine above has the disadvantage to provide an useless (at least for my purposes) solution.

@Igor: thank you for the reply! Your solution involves the eigen values and vectors of $A^{1/2}$, Am I wrong? If not, regrettably your solution is not feasible for my purposes as it does not allow me to develop further calculations.

Thank you Carlo again for the reply. What it is not clear to me is how from that integral we obtain $D=A^{-1/2}B/2$, or equivalently, the final expansion for the square root matrix $A^{1/2}+n^{-1/2}A^{-1/2}B/2$.

Thank you Carlo for the reply. As far as I understand, that solution coincides with the one I posted above, i.e. $\frac{1}{2}\underline\gamma^{ra}c_{as}$ and which, regrettably, is not working (in the sense that it does not solve the equation $\underline\gamma_{ra}\underline c_{as}+\underline c_{ra}\underline\gamma_{as}=c_{rs}$). Any ideas?

That's the whole point. Of course I realize that if I plug back the (wrong) solution in $\underline\gamma_{ra}\underline c_{as}+\underline c_{ra}\underline\gamma_{as}$ I do not get $c_{rs}$, but why symmetry of all the matrices involved does not help to simplify the calculations?

All matrices involved are positive definite and real valued. As you told me and as I have verified the solution is guaranteed to be symmetric as well, therefore the equation should reduce to $2\underline\gamma_{ra}\underline c_{as}=c_{rs}$ and the solution should be $\underline c_{rs}=\frac{1}{2}\underline\gamma^{ra}c_{as}$, where $\underline\gamma^{rs}$ is the inverse of $\underline\gamma_{rs}$. Am I right or Am I missing something here? I'm asking because this was my first (and raw) solution to the problem, but whenever I use such a result to develop further expansions I get strange results.

Thank you Federico for the reply, I am going to see the Lyapunov equation right now! Yes, $c$ is symmetric (I am 100% sure about it, as I have an explicit expression). What do you mean by definite? If I understand correctly, I can tell you that I have the expressions for both $\gamma$ and $c$.