## 2 Theory
Once each instantaneous electric field in the source plane is known, the relation between such fields in the source plane and in the output plane passing through a paraxial ABCD optical system is linked to the Collins integral formula, given by  
\[T_{nz}(\mathbf{\rho})=\left(-\frac{i}{\lambda B}\right)\int T_{n0}(\mathbf{r}) \exp\left[\frac{ik}{2B}\big{(}A\mathbf{r}^{2}-2\mathbf{\rho}\cdot\mathbf{r}+D \mathbf{\rho}^{2}\big{)}\right]d^{2}\mathbf{r}, \tag{7}\]  
where \(k=2\pi/\lambda\) is the wavenumber associated with \(\lambda\), which is the wavelength of the light, and \(\mathbf{\rho}\) denotes the arbitrary position vector in the output plane. \(A,B,D\) are the elements of the transfer matrix of an optical system. On Substituting Equation (6) into (7), we obtain the expression  
\[T_{nz}(\mathbf{\rho})=\left(-\frac{i}{\lambda B}\right)\exp\left(\frac{ikD \mathbf{\rho}^{2}}{2B}\right)\int F_{T}\bigg{[}\sqrt{p(\mathbf{r})}C_{n}( \mathbf{r})\bigg{]}\tau_{0}(\mathbf{r})\exp\left(\frac{ikA\mathbf{r}^{2}}{2B }\right)\times\exp\left(\frac{ik\mathbf{\rho}\cdot\mathbf{r}}{B}\right)d^{2} \mathbf{r}. \tag{8}\]After some rearrangement, Equation (8) can be written in the simple convolution form  
\[T_{nz}(\rho)=\left(-\frac{i}{\lambda B}\right)\exp\left(\frac{ikD\rho^{2}}{2B} \right)\left(\sqrt{\rho}C_{n}\right)\otimes\overline{U}\left(-\frac{\rho}{ \lambda B}\right), \tag{9}\]  
where \(\overline{U}\) is the inverse Fourier transform of the new defined function \(U(\mathbf{r})=\tau_{0}(\mathbf{r})\exp\left(ikAr^{2}/2B\right)\), and \(\otimes\) denotes the convolution symbol. Equation (9) establishes the relation of an instantaneous field between the source plane and the output plane with the convolution form, which is the main result of this paper. It provides one convenient way to simulate the second-order statistics of a GPCB in the output plane. Here, we describe the method for numerical simulating the instantaneous electric field \(T_{nz}(\rho)\) using MATLAB. It follows from Equation (9) that the key point is to calculate the convolution term \(\left(\sqrt{\rho}C_{n}\right)\otimes\overline{U}\). According to the convolution theorem, such term can be expressed as \(F_{T}^{-1}\left[U\times F_{T}\left(\sqrt{\rho}C_{n}\right)\right]\), where \(F_{T}^{-1}\) denotes the symbol for inverse Fourier transform. The operations for the \(F_{T}^{-1}\) and \(F_{T}\) in MATLAB are based on the fast Fourier transform (FFT) algorithm. In numerical simulation, each realization of any function is represented as a matrix with \(N_{1}\times N_{1}\) number of elements. We first assume the function \(U\) has the physical dimensions \(L\times L\), where \(L\) is the length of the coordinate. The corresponding coordinate \(\mathbf{r}=(x,y)\) for \(x\) or \(y\) is from \(-L/2\) to \(L/2\) with step \(L/N_{1}\). In order to match the coordinate \((x,y)\), the function \(p\) (related to the DOC) must sample in the spatial frequency domain having dimensions \(N_{1}/L\times N_{1}/L\), and the corresponding coordinate \(f_{x}\) or \(f_{y}\) is from \(-N_{1}/2L\) to \(N_{1}/2L\), with step \(1/L\)[41]. The random complex function \(C_{n}\) is generated by \([\text{randn}(N_{1})+\text{i}\times\text{randn}(N_{1})]/\text{sqrt}(2)\), where \(\text{randn}(N_{1})\) is the build-in function returning the \(N_{1}\times N_{1}\) matrix consisting of random numbers of normal distribution. Then, the Fourier transform function \(\sqrt{\rho}C_{n}\) is multiplied by the function \(U\), and, finally, \(T_{nz}(\rho)\) is obtained by multiplying the inverse transform of the above results \(U\times F_{T}\left(\sqrt{\rho}C_{n}\right)\) by the coefficient before the convolution shown in Equation (9). The corresponding coordinates for \(\rho_{x}\) or \(\rho_{y}\) in the out plane will be from \(-\lambda BN_{1}/2L\) to \(\lambda BN_{1}/2L\), with step \(\lambda B/L\)[41].