Assume $x[n]$ and $y[n]$ represent two jointly WSS (wide-sense stationary) discrete-time random processes. Further assume that they have zero means, without losing the generality.
If $x[n]$ and $y[n]$ are independent, then they are also uncorrelated, and since they have zero means, they are also orthogonal; i.e.
$$ E\{ x[n]y^*[n-k]\} = 0 \tag{1}$$
Let $$z[n] = x[n] + y[n] \tag{2} $$ be the sum process, then PSD (power spectral density) of $z[n]$ is
$$ S_z(\omega) = FT\{ r_z[m] \} \tag{3} $$
where $r_z[m]$ is the Auto-Correlation sequence of $z[n]$:
$$ r_z[m] = E\{ z[n]z^*[n-m] \} \tag{4} $$
Based on Eq.1 and Eq.2, expand Eq.4 as:
$$
\begin{align}
r_z[m] &= E\{ z[n]z^*[n-m] \} \\
& = E\{ (x[n]+y[n])(x^*[n-m]+y^*[n-m]) \} \\
& = E\{ x[n]x^*[n-m]\} + E\{ y[n]y^*[n-m]\} + \\
&+ E\{ x[n]y^*[n-m]\} + E\{ y^*[n]x[n-m]\} \\
& = r_x[m] + r_y[m] + 0 + 0 \tag{5}
\end{align}
$$
The zeros on the last line follow because of Eq.1, ($x[n]$ and $y[n]$ are othogonal).
Hence what follows is that the PSD of $z[n]$ is the sum of PSDs of $x[n]$ and $y[n]$ :
$$ S_z(\omega) = S_x(\omega) + S_y(\omega) \tag{6}$$
Everything up to this point is theory based on mathematical (probabilty) model of random processes. What you are after, however, is about the practical computation of the PSD estimation and how the independence of $x[n]$ and $y[n]$ is reflected into tho Fourier calculations.
Assuming we use the periodogram estimate for PSD :
$$ \hat{S}_x(\omega) = \frac{1}{N} |X(\omega)|^2 \tag{7}$$
then what we want to show is, as you have also worked out,
$$ \hat{S}_z(\omega) = \hat{S}_x(\omega) + \hat{S}_y(\omega) \tag{8}$$
which implies that
$$ |FT(x+y)|^2 = (X + Y)(X^*+Y^*) = |X|^2 + |Y|^2 + XY^* + X^*Y = |X|^2 + |Y|^2 \tag{9} $$
and Eq.9 implies that:
$$ XY^* + X^*Y = 0 \tag{9}$$
Eventually it's convincing to show that $XY^* = 0$.
To show this is you can use the ergodicty assumption. If $x[n]$ and $y[n]$ are ergodic in their cross-correlation, then the theoretical expectation can be replaced by the practical sums as:
$$E\{x[n]y^*[n-m] \} = \lim_{N \to \infty} \frac{1}{N} \sum_{n=0}^{N-1} x[n]^*y[n-m] \tag{10}
$$
then it can be shown that the estimate is equivalent to convolution of the sequences :
$$
\frac{1}{N} \sum_{n=0}^{N-1} x[n]y^*[n-m] = x[n] \star y^*[-n]\tag{11}
$$
But from Eq.1 we know that the $x[n]$ and $y[n]$ are orthogonal implying that
$$ x[n] \star y^*[-n] = 0 \tag{12}$$
finally, if this sum is zero then it's Fourier transform also is zero, completing the proof.
$$ x[n] \star y^*[n] = 0 \implies X[k]Y^*[k] = 0 $$
Note that, as MArcus has already commented, the equality holds only for the ideal case of infinite length process realization and for any practical finite length sample sequence, the Eq.s (6,7,8,9,10,12) won't hold true;i.e., they will be approximate).