## 2 Theory  
Let us consider a scalar, quasi-monochromatic, and statistically stationary optical beam, propagating along the z-direction. In the space-frequency domain, the second-order statistics of such GPCB are typically characterized by the CSD function \(W(\rho_{1},\rho_{2},z)=\langle E^{*}(\rho_{1},z)E(\rho_{2},z)\rangle\), where \(E\) denotes the random electric field, \(\rho_{1}\) and \(\rho_{2}\) are two arbitrary position vectors in the transverse plane, the asterisk is the complex conjugate, and the angular brackets is the ensemble average over the field fluctuations. Assume that the GPCB we consider is a Schell-model beam whose DOC function is only dependent on the separation of two points. In this situation, the CSD function in the source plane (z = 0) can be written as [15; 16]  
\[W_{0}(\mathbf{r}_{1},\mathbf{r}_{2})=\tau_{0}(\mathbf{r}_{1})\tau_{0}^{*}( \mathbf{r}_{2})\mu_{0}(\mathbf{r}_{1}-\mathbf{r}_{2}), \tag{1}\]where \(\mathbf{r}_{1}\) and \(\mathbf{r}_{2}\) are two arbitrary position vectors in the source plane, \(\tau_{0}\) is a complex amplitude function, and \(\mu_{0}\) is the DOC function. For a physically realizable CSD function, the DOC can be alternatively represented the following integral form [19]  
\[\mu_{0}(\mathbf{r}_{1}-\mathbf{r}_{2})=\int W_{s}(\mathbf{v}_{1},\mathbf{v}_{2 })\exp[-i2\pi(\mathbf{r}_{1}-\mathbf{r}_{2})\cdot\mathbf{v}]d^{2}\mathbf{v}, \tag{2}\]  
with  
\[W_{s}(\mathbf{v}_{1},\mathbf{v}_{2})=\sqrt{p(\mathbf{v}_{1})}\sqrt{p(\mathbf{ v}_{2})}\delta(\mathbf{v}_{1}-\mathbf{v}_{2}). \tag{3}\]  
The function \(p(\mathbf{v})\) in Equation (3) must be non-negative for any \(\mathbf{v}\) to ensure that the CSD is genuine; \(\delta\) denotes the Dirac function. In order to introduce our numerical approach, we rewrite the \(\delta\) function in the follow form  
\[\delta(\mathbf{v}_{1}-\mathbf{v}_{2})=\langle C_{n}(\mathbf{v}_{1})C_{n}^{*} (\mathbf{v}_{2})\rangle, \tag{4}\]  
where \(C_{n}(\mathbf{v})\) (\(n=1,2,3,\dots\) ) denotes the random complex function whose probability density functions of amplitude and phase of \(C_{n}\) obey the negative exponent and uniform distribution, respectively [39]. The angular bracket denotes the ensemble average. Here, \(C_{n}(\mathbf{v})\) can be viewed as the field of the white noise. By substituting Equation (4) into Equations (1)-(3) and after some mathematical operations, Equation (1) turns out to be  
\[W_{0}(\mathbf{r}_{1},\mathbf{r}_{2})=\left\langle T_{n0}(\mathbf{r}_{1})T_{n0 }^{*}(\mathbf{r}_{2})\right\rangle\approx\frac{1}{N}\sum_{n0=1}^{N}T_{n0}( \mathbf{r}_{1})T_{n0}^{*}(\mathbf{r}_{2}), \tag{5}\]  
with  
\[T_{n0}(\mathbf{r}_{i})=\tau_{0}(\mathbf{r}_{i})\times F_{T}\bigg{[}\sqrt{p( \mathbf{r}_{i})}C_{n}(\mathbf{r}_{i})\bigg{]},i=1,2. \tag{6}\]  
where \(F_{T}\) stands for the symbol of Fourier transform, \(T_{n0}\) (\(\mathbf{r}\)) (\(n=1,2,3\dots\) ) denotes the instantaneous electric fields, each of which represents one realization of a GPCB. In Equation (5), we use the incoherent superposition of instantaneous random realizations instead of ensemble averaging over the random fields, since we assume the beam is statistically stationary. If the number of realizations \(N\) is large enough, the incoherent realization superposition is a good approximation. The first term \(\tau_{0}\) in Equation (6) is the complex amplitude, as we described above, and the second term contains the information of the DOC of the source. To synthesize the realization of \(T_{n0}\) (\(\mathbf{r}\)) numerically, one can use commercial software such as MATLAB to generate the random complex function \(C_{n}\) and then take the Fourier transform of the product of the function \(C_{n}\) and the square root of the \(p(\mathbf{v})\) function. Finally, one instantaneous electric filed is acquired by multiplying the Fourier result by the function \(\tau_{0}\). Note that the DOC is determined by the function \(p(\mathbf{v})\) in the simulation.