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 that in the original equation $X$ must be symmetric and so does $\tilde{X}$. 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?