I figure out something and would like to share and check whether it is correct.
It is perfectly analogous to real quadratic equations.
Since $L$ is symmetric and positive definite, they can be decomposed:
$L=U_L^{\prime} \Lambda_L U_L$ where $U_L$ is orthonormal, and $\Lambda_L$ is diagonal with positive entries.
Let $Q=U_L^{\prime} \Lambda_L^{-1/2} U_L$, which is the inverse of the square root of $L$, and it is symmetric.
Then the above equation is equivalent to
$\tilde{X}^{\prime} \tilde{X}-Q\tilde{X}=D$ where $\tilde{X}=Q^{-1}X$
Note thatSuppose, in the original equationaddition, $X$ must be symmetric$D$ and so does $L$ are simultaneously diagonalizable, I conjecture that $\tilde{X}$ is symmetric. Therefore, the above is equivalent to
$(\tilde{X}-\frac{1}{2}Q)^{\prime}(\tilde{X}-\frac{1}{2}Q)=\frac{1}{4}L^{-1}-D$
Note that the right-hand side is symmetric; hence, there exists a real solution of $\tilde{X}$ if and only if the right-hand side is positive semi-definite.
And after a bit of algebra, if solutions exist, it would be
$X=\frac{1}{2}L^{-1}[I \pm (I-4LD)^{1/2}]$ where $(I-4LD)^{1/2}$ is the real p.s.d square root of $(I-4LD)$, which exists because $(I-4LD)^{1/2}$ is p.s.d. and symmetric.
Hence, there are at most two solutions.
Are all of these arguments correct?