Let $\{X(t): -\infty < t < \infty\}$ denote a low-pass band-limited WSS white noise process of bandwidth $B$ Hz and total noise power $N_0B$. Thus, the $X(t)$'s are a collection of zero-mean random variables, and the power spectral density of the process is $$S_X(f) = \dfrac{N_0}{2}\operatorname{rect}\left(\dfrac{f}{2B}\right)\tag{1}$$ while its autocorrelation function is $$R_X(\tau) = E[X(t)X(t+\tau)] = N_0B \operatorname{sinc}(2Bt).\tag{2}$$ Note that all the $X(t)$'s are zero-mean random variables with variance $N_0B$ and that $X(t)$ and $X(t+\tau)$ are uncorrelated if and only if $\tau$ is a nonzero integer multiple of $\frac{1}{2B}$. For other values of $\tau$, $X(t)$ and $X(t+\tau)$ are correlated random variables.

Now suppose that $\{Y(t): -\infty < t < \infty\}$ is a related random process such that $Y(t) = X^2(t)$ for all $t, -\infty < t < \infty.$ The OP asks: what is $S_Y(f)$,the PSD of $\{Y(t)\}$? Well, let's determine $R_Y(\tau)$ first. We have that $$R_Y(\tau) = E[Y(t)Y(t+\tau)] = E[X^2(t)X^2(t+\tau)]$$ but generally cannot proceed further because we don't know the joint distribution of $X(t)$ and $X(t+\tau)$. However, when $\{X(t)\}$ is assumed to be a Gaussian white noise process, we have that the random variables $X(t)$ and $X(t+\tau)$ enjoy a jointly Gaussian distribution and it is known (see the section Properties: Higher moments in the Wikipedia link provided) that \begin{align}E[X^2(t)X^2(t+\tau)] &= \operatorname{var}\big(X(t)\big)\operatorname{var}\big(X(t+\tau)\big)+2\operatorname{cov}\big(X(t),X(t+\tau)\big)\\ &= (N_0B)^2 + 2N_0B \operatorname{sinc}(2Bt).\end{align} Applying $(1)$ and $(2)$, we see that $$S_Y(f) = (N_0B)^2\delta(f) + N_0\operatorname{rect}\left(\dfrac{f}{2B}\right)\tag{3}$$ with an impulse at $f=0$ as mentioned in my comment on the main question. Note that the total power has increased from $N_0B$ to $(N_0B)^2 + 2N_0B$ but the PSD still has support $[-B,B]$, that is, squaring has not increased the bandwidth.

Let $\{X(t): -\infty < t < \infty\}$ denote a low-pass band-limited WSS white noise process of bandwidth $B$ Hz and total noise power $N_0B$. Thus, the $X(t)$'s are a collection of zero-mean random variables, and the power spectral density of the process is $$S_X(f) = \dfrac{N_0}{2}\operatorname{rect}\left(\dfrac{f}{2B}\right)\tag{1}$$ while its autocorrelation function is $$R_X(\tau) = E[X(t)X(t+\tau)] = N_0B \operatorname{sinc}(2Bt).\tag{2}$$ Note that all the $X(t)$'s are zero-mean random variables with variance $N_0B$ and that $X(t)$ and $X(t+\tau)$ are uncorrelated if and only if $\tau$ is a nonzero integer multiple of $\frac{1}{2B}$. For other values of $\tau$, $X(t)$ and $X(t+\tau)$ are correlated random variables.

Now suppose that $\{Y(t): -\infty < t < \infty\}$ is a related random process such that $Y(t) = X^2(t)$ for all $t, -\infty < t < \infty.$ The OP asks: what is $S_Y(f)$,the PSD of $\{Y(t)\}$? Well, let's determine $R_Y(\tau)$ first. We have that $$R_Y(\tau) = E[Y(t)Y(t+\tau)] = E[X^2(t)X^2(t+\tau)]$$ but generally cannot proceed further because we don't know the joint distribution of $X(t)$ and $X(t+\tau)$. However, when $\{X(t)\}$ is assumed to be a Gaussian white noise process, we have that the random variables $X(t)$ and $X(t+\tau)$ enjoy a jointly Gaussian distribution and it is known (see the section Properties: Higher moments in the Wikipedia link provided) that \begin{align}E[X^2(t)X^2(t+\tau)] &= \operatorname{var}\big(X(t)\big)\operatorname{var}\big(X(t+\tau)\big)+2\operatorname{cov}\big(X(t),X(t+\tau)\big)\\ &= (N_0B)^2 + 2N_0B \operatorname{sinc}(2Bt).\end{align} Applying $(1)$ and $(2)$, we see that $$S_Y(f) = (N_0B)^2\delta(f) + N_0\operatorname{rect}\left(\dfrac{f}{2B}\right)\tag{3}$$ with an impulse at $f=0$ as mentioned in my comment on the main question. Note that the total power has increased from $N_0B$ to $(N_0B)^2 + 2N_0B$ but the PSD still has support $[-B,B]$, that is, squaring has not increased the bandwidth. This is in sharp contrast to what happens with (finite energy) deterministic signals for which squaring generally increases the bandwidth, and so take all the above with a large grain of salt!

Let $\{X(t): -\infty < t < \infty\}$ denote a low-pass band-limited WSS white noise process of bandwidth $B$ Hz and total noise power $N_0B$. Thus, the $X(t)$'s are a collection of zero-mean random variables, and the power spectral density of the process is $$S_X(f) = \dfrac{N_0}{2}\operatorname{rect}\left(\dfrac{f}{2B}\right)\tag{1}$$ while its autocorrelation function is $$R_X(\tau) = E[X(t)X(t+\tau)] = N_0B \operatorname{sinc}(2Bt).\tag{2}$$ Note that all the $X(t)$'s are zero-mean random variables with variance $N_0B$ and that $X(t)$ and $X(t+\tau)$ are uncorrelated if and only if $\tau$ is a nonzero integer multiple of $\frac{1}{2B}$. For other values of $\tau$, $X(t)$ and $X(t+\tau)$ are correlated random variables.

Now suppose that $\{Y(t): -\infty < t < \infty\}$ is a related random process such that $Y(t) = X^2(t)$ for all $t, -\infty < t < \infty.$ The OP asks: what is $S_Y(f)$,the PSD of $\{Y(t)\}$? Well, let's determine $R_Y(\tau)$ first. We have that $$R_Y(\tau) = E[Y(t)Y(t+\tau)] = E[X^2(t)X^2(t+\tau)]$$ but cannot evaluate the expectation furthering the general case. However, when $\{X(t)\}$ is assumed to be a Gaussian white noise process, the random variables $X(t)$ and $X(t+\tau)$ enjoy a bivariate jointly Gaussian distribution and it is known (see section Properties: Higher moments) that \begin{align}E[X^2(t)X^2(t+\tau)] &= \operatorname{var}\big(X(t)\big)\operatorname{var}\big(X(t+\tau)\big)+2\operatorname{cov}\big(X(t),X(t+\tau)\big)\\ &= (N_0B)^2 + 2N_0B \operatorname{sinc}(2Bt).\end{align} Applying $(1)$ and $(2)$, we see that $$S_Y(f) = (N_0B)^2\delta(f) + N_0\operatorname{rect}\left(\dfrac{f}{2B}\right)\tag{3}$$ with an impulse at $f=0$ as mentioned in my comment on the main question. Note that the total power has increased from $N_0B$ to $(N_0B)^2 + 2N_0B.$

Let $\{X(t): -\infty < t < \infty\}$ denote a low-pass band-limited WSS white noise process of bandwidth $B$ Hz and total noise power $N_0B$. Thus, the $X(t)$'s are a collection of zero-mean random variables, and the power spectral density of the process is $$S_X(f) = \dfrac{N_0}{2}\operatorname{rect}\left(\dfrac{f}{2B}\right)\tag{1}$$ while its autocorrelation function is $$R_X(\tau) = E[X(t)X(t+\tau)] = N_0B \operatorname{sinc}(2Bt).\tag{2}$$ Note that all the $X(t)$'s are zero-mean random variables with variance $N_0B$ and that $X(t)$ and $X(t+\tau)$ are uncorrelated if and only if $\tau$ is a nonzero integer multiple of $\frac{1}{2B}$. For other values of $\tau$, $X(t)$ and $X(t+\tau)$ are correlated random variables.

Now suppose that $\{Y(t): -\infty < t < \infty\}$ is a related random process such that $Y(t) = X^2(t)$ for all $t, -\infty < t < \infty.$ The OP asks: what is $S_Y(f)$,the PSD of $\{Y(t)\}$? Well, let's determine $R_Y(\tau)$ first. We have that $$R_Y(\tau) = E[Y(t)Y(t+\tau)] = E[X^2(t)X^2(t+\tau)]$$ but generally cannot proceed further because we don't know the joint distribution of $X(t)$ and $X(t+\tau)$. However, when $\{X(t)\}$ is assumed to be a Gaussian white noise process, we have that the random variables $X(t)$ and $X(t+\tau)$ enjoy a jointly Gaussian distribution and it is known (see the section Properties: Higher moments in the Wikipedia link provided) that \begin{align}E[X^2(t)X^2(t+\tau)] &= \operatorname{var}\big(X(t)\big)\operatorname{var}\big(X(t+\tau)\big)+2\operatorname{cov}\big(X(t),X(t+\tau)\big)\\ &= (N_0B)^2 + 2N_0B \operatorname{sinc}(2Bt).\end{align} Applying $(1)$ and $(2)$, we see that $$S_Y(f) = (N_0B)^2\delta(f) + N_0\operatorname{rect}\left(\dfrac{f}{2B}\right)\tag{3}$$ with an impulse at $f=0$ as mentioned in my comment on the main question. Note that the total power has increased from $N_0B$ to $(N_0B)^2 + 2N_0B$ but the PSD still has support $[-B,B]$, that is, squaring has not increased the bandwidth.