I am trying to draw random numbers
$$Z_{l,m} = \int_{-1}^1 dx \, P_l^m(x)W(x)$$
Here $P_l^m(x)$ are the associated Legendre polynomials with integers $l\geq0$ and $-l\leq m \leq l$. The variable $W(x)$ corresponds to white noise with zero mean and variance
$$W(x)W(x')=\delta(x-x')$$
Here $\delta(x-x')$ is the delta distribution.
I noticed already, that $P_l^m$ is proportional to $P_l^{-m}$, hence I only need to draw numbers for $m\geq0$.
However, I am not sure at all if the remaining random numbers $Z_l^m$ are independent. Does anyone have ideas on either how to show the independence or on how to draw the random numbers?
I think it is better to avoid thinking in terms of bona fide white noise, as physicists often think. You should instead think in terms of computing
$$I_l^m=\int_0^2 P_l^m(x-1) dW(x).$$
This is a stochastic integral. (I shift everything to start at zero for consistency with standard notation in math.)
Approach 1:
Rewrite the stochastic integral by "formal integration by parts":
$$I_l^m=P_l^m(1) W(2)-P_l^m(-1)W(0)-\int_0^2 W(x) (P_l^m)'(x-1) dx \\= P_l^m(1) W(2) - \int_0^2 W(x) (P^m_l)'(x-1) dx.$$
The second term is now a regular integral, not a stochastic integral anymore. You can now approximate this integral as:
$$\sum_{k=1}^N (P_l^m)'\left ( \frac{2k}{N}-1 \right ) \sum_{j=1}^k N_j$$
where $N_j$ are iid normal random variables with mean zero and variance $\frac{2}{N}$. This is essentially the rectangle rule for the regular integral, with the inner sum serving to approximate a sample path of the Wiener process.