A Simplified Model for Optical Systems with Random Phase Screens

A first-order optical system with arbitrary multiple masks placed at arbitrary positions is the basic scheme of various optical systems. Generally, masks in optical systems have a non-shift invariant (SI) effect; thus, the individual effect of each mask on the output cannot be entirely separated. The goal of this paper is to develop a technique where complete separation might be achieved in the common case of random phase screens (RPSs) as masks. RPSs are commonly used to model light propagation through the atmosphere or through biological tissues. We demonstrate the utility of the technique on an optical system with multiple RPSs that model random scattering media.


Introduction
An optical system with multiple arbitrary masks placed at arbitrary positions is the core scheme of a variety of optical systems. A scheme where the masks are random phase screens (RPSs) is an excellent model for the analysis, simulation, and interpretation of optical systems that involve optical scattering due to refractive index fluctuations [1][2][3][4]. RPS models were first introduced to represent the fading of radio waves due to fluctuations in the ionosphere layers, but were later found to be useful for modeling various other wave scattering phenomena, including optical scattering. They have been extensively used to analyze atmospheric light propagation, coherent imaging, speckle metrology, and optical scintillations, amongst others [5].
In Figure 1a, a scheme of a paraxial optical system is illustrated, where there are n sub paraxial systems, M 1 , M 2 . . . M k . . . M n , and n − 1 masks, g 1 , g 2 . . . g k−1 , g k . . . g n−1 , between them, while U in and U out denote the input field and output field, respectively. For example, in atmospheric propagation of light, the sub-systems might be free-space propagation combined with an optical scaling system where the masks are RPSs [4,6]. Figure 1b is a specific case of Figure 1a for two masks and three sub-systems, showing a paraxial system composed of intermittent lenses of focal length f and free-space propagation sections of length L/2 with two masks in between. By choosing the masks as RPSs with properties appropriate for the relevant medium, this figure represents a model [5] for the analysis of the optical memory effect of scattered light in random media [7]. We use Figure 1b as a case study throughout this paper. The basic principle behind the optical memory has its roots in astronomy and was adapted to other fields using well-established astronomic techniques [5,8]. Some other advanced optical memory effect applications might include imaging through the atmosphere, dense fog, or biological tissues [9].
The first order paraxial sub-systems between the masks in Figure 1a might be any combination of beam propagation forms with ideal optical elements between the stages, such as ideal thin lenses, ideal graded-index fiber, etc. Such lenses and fibers are commonly used for scaling [4], imaging [10], and spatial filtering [5] between the RPSs in the system. For example, in Figure 1b, M 1 is a free-space propagation section of length L/2, M 2 is an ideal 2 f focal length system, and M 3 is an ideal 2 f focal length system with an additional is an ideal 2 focal length system, and is an ideal 2 focal length system with an additional free-space propagation section of length 2 ⁄ . Thus, we define the optical system without the masks as the "ideal core" optical system. The model may also contain a combination of RSPs with random absorbing screens. , . . . placed at arbitrary positions. and denote the input field and output field. (b) A private case of (a) for two random phase masks and three sub-systems. By choosing appropriate system properties, the set-up represents a medium with an optical memory effect.
Generally, masks in optical systems have non shift-invariant (SI) behavior; thus, the individual effect of each mask on the output cannot be entirely separated. However, we will show in this paper that the effect of each RPS in the system in Figure 1a can be fully separated under certain conditions. This separation facilitates the analysis of various systems, and helps to understand the relation between the input and the output. Furthermore, it enables separating the influence of each of the system's parameters and enables modeling more complex systems, such as random media with varying behavior along the system path.
This paper is structured as follows. In Section 2, we investigate the general case of two RSPs, and we demonstrate our approach for the specific case study illustrated in Figure 1b. For completeness of the analysis and for the purpose of comparison, we first review the case of a deterministic mask (Section 2.1), and then we develop the main result of the RPS case and demonstrate it by simulation (Section 2.2). In Section 3, we generalize the approach to the multiple RPS case and its applications. Section 4 presents a discussion and conclusions.

Review of Optical Systems with Deterministic Masks
In Figure 2a, a special case of Figure 1a is shown for two masks and three sub-systems; here, the masks are considered to be any general known masks and are not necessarily RPSs. General masks have various optical applications including incoherent holography [11,12], compressive sensing [13], and spectral processing [14].
, and represent sub-paraxial systems before, between, and after masks and , where denotes the coherent input and denotes the coherent output. In Figure 2b, the scheme of the ideal core system, without the masks, is shown; its coherent output is denoted by . The output field in Figure 2c denotes the projection of the mask through the ideal core system from its position to the output plane. Similarly, in Figure 2d is the output projection of the mask . The output field can be calculated by a sequence of generalized convolutions [15] of the output field core system with the projection . . M n with arbitrary multiple masks g 1 , g 2 . . . g k−1 , g k . . . g n−1 placed at arbitrary positions. U in and U out denote the input field and output field. (b) A private case of (a) for two random phase masks and three sub-systems. By choosing appropriate system properties, the set-up represents a medium with an optical memory effect.
Generally, masks in optical systems have non shift-invariant (SI) behavior; thus, the individual effect of each mask on the output cannot be entirely separated. However, we will show in this paper that the effect of each RPS in the system in Figure 1a can be fully separated under certain conditions. This separation facilitates the analysis of various systems, and helps to understand the relation between the input and the output. Furthermore, it enables separating the influence of each of the system's parameters and enables modeling more complex systems, such as random media with varying behavior along the system path.
This paper is structured as follows. In Section 2, we investigate the general case of two RSPs, and we demonstrate our approach for the specific case study illustrated in Figure 1b. For completeness of the analysis and for the purpose of comparison, we first review the case of a deterministic mask (Section 2.1), and then we develop the main result of the RPS case and demonstrate it by simulation (Section 2.2). In Section 3, we generalize the approach to the multiple RPS case and its applications. Section 4 presents a discussion and conclusions.

Review of Optical Systems with Deterministic Masks
In Figure 2a, a special case of Figure 1a is shown for two masks and three subsystems; here, the masks are considered to be any general known masks and are not necessarily RPSs. General masks have various optical applications including incoherent holography [11,12], compressive sensing [13], and spectral processing [14]. M 1 , M 2 and M 3 . represent sub-paraxial systems before, between, and after masks g 1 and g 2 , where U in denotes the coherent input and U out denotes the coherent output. In Figure 2b, the scheme of the ideal core system, without the masks, is shown; its coherent output is denoted by U ideal . The output field U g 1 in Figure 2c denotes the projection of the mask g 1 through the ideal core system from its position to the output plane. Similarly, U g 2 in Figure 2d is the output projection of the mask g 2 . The output field U out can be calculated by a sequence of generalized convolutions [15] of the output field core system U ideal with the projection output of U g 1 and U g 2 , which are defined mathematically in Appendix A. Although the generalized convolution provides an elegant and compact way to express the output field, in this paper, we will carry out our analysis using conventional convolutions (the generalized convolution can be expressed in terms of conventional convolutions [16,17]). The output field U out in the scheme in Figure 2a is obtained, up to a multiplication factor, as [18]: where  [19,20]. For example, a scaled optical FT of mask g(x) by a 2 f system is given, up to a complex factor, x 2 , where b 3 and d 3 originate from the ABCD ray transfer matrix of the paraxial sub-system M 3 for the projection of g 2 onto the output plane ( Figure 2d). Similarly, b 32 and d 32 originate from the ABCD ray transfer matrix of the cascade sub-system [18] of M 2 and M 3 for the projection of g 1 onto the output plane ( Figure 2c). In Equation (1), and henceforth, we adopt a one-dimensional notation for simplicity, and without loss of generality. output of and , which are defined mathematically in Appendix A. Alt generalized convolution provides an elegant and compact way to express the ou in this paper, we will carry out our analysis using conventional convolutions ( alized convolution can be expressed in terms of conventional convolutions [1 output field in the scheme in Figure 2a is obtained, up to a multiplication [18]: is the wavelength, * is the convolution and ℑ is the inverse Fourier Transform (FT) operator. V 1⁄ ℎ( ) is a scalin on a general complex function, with ℎ( ) defined as V 1⁄ ℎ( ) = ℎ (1⁄ ) [19,20]. For example, a scaled optical FT of mask ( ) by a 2 system is give , where and origina ABCD ray transfer matrix of the paraxial sub-system for the projection of output plane (Figure 2d). Similarly, and originate from the ABCD ra matrix of the cascade sub-system [18] of and for the projection of on put plane (Figure 2c). In Equation (1), and henceforth, we adopt a one-dimens tion for simplicity, and without loss of generality. Equation (1) was derived mathematically in [18], but it would be instructiv ute some physical meaning to it. The term V ℑ plays a similar role to coherent impulse response in a 4 system with in the Fourier plane [20]. T Equation (1) was derived mathematically in [18], but it would be instructive to attribute some physical meaning to it. The term V 1 λb −1 [g] plays a similar role to that of the coherent impulse response in a 4 f system with g in the Fourier plane [20]. The term q accounts for generalization to other optical systems than the 4 f example. In a 4 f system, the parameter b is equal to the focal length. For other optical systems, the parameters b and d need to be calculated from the ABCD projection system parameters, as mentioned above. The concept of the forward projections, as shown in Figure 2c practice to project each one towards the image region and the most restricting one is defined as the exit pupil. Then, the scaled FT of the exit pupil can be used to approximate the coherent impulse response. The Fourier scaling parameter is related to the ratio between the aperture and the exit pupil. This conceptual idea is manifested in the field analysis in Equation (1) by the FT of each mask with its own scaling parameter b. The parameter d for each projection is a measure of the shifting of each mask from the Fourier plane, and it approaches zero as the masks are shifted to the Fourier plane. One advantage of the formalism in Equation (1) is the separation of the influence of U ideal from the expression of U out . Unfortunately, only in the particular case of a system with q 32 = q 3,32 = 1 can the effect of each mask be fully separated from any other mask [18], since then Equation (1) reduces to a sequence of pure convolutions and, therefore, the commutative property of the regular convolution operator applies. In this paper, we generalize the effect of each mask and investigate the conditions for such separation for the case of RPSs.

Separating the Effect of Each Random Phase Mask in the Optical System
Consider the case where the masks in Figure 2a are RPSs, and we are interested in the output intensity for a coherent input. The output intensity measurement is the desired result in common optic systems with random media, such as biological and atmospheric light propagation applications [1][2][3][4]21,22]. A practical simulation example of this will be given in the next section. We assume that the RPSs are spatially stationary random processes and that only their statistical properties are known, such as their statistical autocorrelation (AC); therefore, only the average effect of the RPS is measured (Ch.8 in [21]). The commonly used result is the average of many measurements of an ensemble of many RPSs for the same set-up [21]. Consequently, we consider the expected intensity of the system in Figure 2a, modeled by Equation (1) with g 1 and g 2 being RPSs. The derivation is described in Appendix B, and uses the method for developing the Schell theorem (Chap. 5. in [21]), following well-known techniques for calculating the expected value of a statistical process through linear systems (Ch.7 in [23]). For the derivation, we defined the mask g(x) as the multiplication of a deterministic mask p(x) with a RPS t(x), i.e., g 1 (x) = p 1 (x)t 1 (x) and g 2 (x) = p 2 (x)t 2 (x). The deterministic masks p 1 and p 2 may represent, for example, simple apertures or diffractive elements. We start with the case of p 1 = p 2 = 1, and we obtain (Appendix B) where E[•] is the expected value notation, −1 . is the inverse FT operator, and Γ t 1 and Γ t 2 are the statistical AC functions of the first and second RPSs t 1 and t 2 , respectively, where Thus, the average output intensity E U out 2 is obtained by a regular convolution of the intensity's ideal core output |U ideal | 2 with the FT of a scaled AC of each RPS. The parameters b 32 and b 3 reflect the influence of the projection of each RPS on the output, according to its position in the system. This influence is analogous to the concept, as described in Figure 2c,d. These are the same scaling parameters of the FT of g 1 and g 2 for the deterministic masks in Equation (1). As is evident from Equation (2), owing to the convolution operator's commutative property, the average effect of each RPS can be entirely separated from any other RPS and the core output intensity |U ideal | 2 . In Equation (2), and henceforth, we ignore constant multiplicative factors. We continue to the case of p 1 = 1, p 2 = 1 and we obtain up to complex factor (Appendix B), Please note that only the effect of −1 [Γ t 1 (λb 32 x)] * −1 [Γ t 2 (λb 3 x)] is SI, and not the whole system. Even if the output intensity of the core system |U ideal | 2 is SI, then, in contradiction to Equation (2), the influence of the deterministic part p 1 makes the whole core system a non-SI system, i.e., Yet, the effect for each RPS still remains SI (since it is described by a convolution in (3)).
Consider the particular case where the core system is a 4 f system, and the input is an on-axis incoherent point source. This is the case in Figure 1b with L = 0 and with only one mask, t 2 , located in the Fourier plane, and without t 1 (i.e., p 1 · t 1 = 1). For this case, Equation (2) reduces to the well-known impulse case of the one random screen analyzing system with an impulse response of (Ch.8 in [21]) is the average optical transfer function (AOTF), and the scaled inverse FT of 2 is the PSF. Similarly, for the general case in Equation (2), the FT of a scaled AC of each RPS might be considered as its average point spread function (APSF), and their convolution as the total APSF. We continue to the general case of p 1 = 1, p 2 = 1, under the assumption that p 2 (x)t 2 (x) is a stationary process. In such a case we may obtain (Appendix B):

Simulation of RPS Separation
We shall demonstrate the application of Equations (2)-(4) through the system in Figure 3a. With this scheme, a focusing system focuses a beam inside a random medium and the deep tissue focusing plane is imaged by a lens to a sensor, where the whole imaging system focus is found by moving the sensor. Such a focusing system is commonly applied by using an adaptive optics system [7] based on spatial light modulator for modulating the beam amplitude and phase. Consider the practical case [7] where the focusing system focuses inside a medium that might be represented by our case study in Figure 1b [5]. Thus, Figure 3a can be modeled by Figure 3b. This set-up is a special case of a paraxial system with two RPSs; therefore, Equations (2)-(4) express its APSF.
A simulation analysis of Equations (2)-(4) is clearer in the spatial frequency domain rather than the spatial domain, especially for the following graph visualities purpose. We choose the case where the ideal output is an on-axis point source and by taking the FT of Equations (2)-(4) the AOTF is obtained. For the case of p 1 = p 2 = 1 in Equation (2), it is obtained by where ν is the spatial frequency variable. The AOTF for the case of p 1 = 1, p 2 = 1, in Equation (3) is obtained by And, similarly, the AOTF for the case of p 1 = 1, p 2 = 1 in Equation (4) is obtained by  The medium was modeled according to the case study in Fig 1b, with the RPSs modeled as Gaussian. Following [5], the AC of a Gaussian RPS is where ′ is the variance of the derivative of the phase [21]. For our case study, the variances of the first and second RPSs are = (2 / ) / , = (2 / ) /(12 ), respectively [5], where is the wavelength, is the lenses' focal length in Figure 1b, is the width of the diffusive medium, and is the transport mean free path. The transport mean free path is associated with the distance through which light propagates between each scattering event and with the anisotropy factor of the medium [3,7].
2 ⁄ is the distance before and after the 4 system in Figure 1b. The system without the masks is just a relay system due to the 4 system. Consequently, the output field is just the result of freespace propagation of the input field by a distance , as expected for propagation with negligible scattering along the slab. The AOTF accounts for the spatial frequency distribution but does not take into account the intensity loss (by a factor of [7]). The medium was modeled according to the case study in Figure 1b, with the RPSs modeled as Gaussian. Following [5], the AC of a Gaussian RPS is φ is the variance of the derivative of the phase [21]. For our case study, the variances of the first and second RPSs are σ 2 1 = L(2π/λ) 2 /l tr , σ 2 2 = L 3 (2π/λ) 2 / 12l tr f 2 , respectively [5], where λ is the wavelength, f is the lenses' focal length in Figure 1b, L is the width of the diffusive medium, and l tr is the transport mean free path. The transport mean free path is associated with the distance through which light propagates between each scattering event and with the anisotropy factor of the medium [3,7]. L/2 is the distance before and after the 4 f system in Figure 1b. The system without the masks is just a relay system due to the 4 f system. Consequently, the output field is just the result of free-space propagation of the input field by a distance L, as expected for propagation with negligible scattering along the slab. The AOTF accounts for the spatial frequency distribution but does not take into account the intensity loss (by a factor of 9 L 8 λl tr π 4 [7]). Following the experimental conditions reported in [7], we chose L = 260 um, l tr = 14.8 mm, and λ = 632.8 nm. However, it should be emphasized that our case study of the memory effect is valid for any scattering medium or geometry [7]. The focal length f in our model (Figure 1b) is a free parameter that affects the scaling parameter b 3 Hence, we choose f to be the same as for all the other cases of free-space propagation in Figure 1b; this choice (i.e., f = L/2) makes it possible to approach the critical sampling for the simulation of free-space propagation [24] and to enhance the simulation accuracy. Consequently, we chose the transfer function simulation approach and not the impulse response approach to better simulate the free-space propagation (Ch.5 in [24]). The Gaussian RPSs were simulated by convolving an uncorrelated random signal with a Gaussian function in a similar way to the Gaussian Schell-model beam simulations [3,24,25].
The Schell theorem is a special case of Figure 1a with one RPS, where p is a clear aperture. Consider p to be a rectangle aperture; then, the diffraction pattern in the Schell theorem is more pronounced as the ratio between the size of p and the transverse coherence length decreases (Chap. 5. in [21]). This is similar to the effect of p 1 and p 2 in Equation (4). To demonstrate this principle, we chose a rectangle size of p 1 and p 2 of 2.82σ 1 −1 and 0.31σ 2 −1 , respectively. We simulated the expected value of Equation (1) on the ensemble of 300 different RPSs, but with the same statistical properties. Figure 3c shows one of the RPS simulations for t 1 from all 300 simulations of the whole ensemble. Similarly, Figure 3d shows one simulation for t 2 and Figure 3e describes the AOTF simulations after normalization and shows a good match between the simulations and the analytical expressions of Equations (5)- (7). The deterministic OTF of a rectangle aperture has a triangle shape. The fifth graph takes into account the two different size apertures, p 1 and p 2 , in addition to the Gaussian AOTF of t 1 and t 2 . This of course reduces the cutoff spatial frequency.

Generalization for Multiple RPSs
We may generalize the former equations for the system with multiple RPSs shown in Figure 1a. Consider a paraxial system with n sub-systems, M 1 , M 2 . . . M k . . . M n , and with n − 1 RPSs t 1 , t 2 . . . t k−1 , t k . . . t n−1 between the sub-systems. Then, (2) is generalized to where t k−1 is the k − 1 th mask and b n,k originates from the projected ABCD ray transfer matrix of the cascade ideal core sub-systems from the t k−1 position to the output plane.
This generalization is useful for many cases. For example, Figure 4a can describe the system model of two adjacent random media with different transport mean free paths. Such a combination of random media is modeled by the system shown in Figure 4b   properties, which might be analyzed by Equation (8). In the analysis of atmospheric light propagation, this may be relevant when the source or the sensor get closer to or further from each other on the same line of sight. In a biological application, this may be relevant, for example, for evaluating the change in the optical propagation through a multi-layered tissue after removing tissue layers during surgery or due to stretching of the tissue. Although there are now four masks, the overall APSF is still a convolution chain of the effect of each APSF medium, which is one of the advantages of Equation (8). Consequently, for example, if we change the value of l tr of only one medium, the APSF of only this medium is changed. Similarly, suppose we add a third medium before the two media in Figure 4a; in this case, the new overall APSF overall is just a convolution of the overall APSF two of the previous two media with the APSF third of the third medium, i.e., APSF overall = APSF third * APSF two . The average effect of each RPS depends on its position in the system in accordance with the scaling factor of the AC. Thus, if we change the order of the two media or their length in Figure 4a, the effect of each RPS is changed only by scaling.
Another useful case is where the random medium properties change along the system path, as is common in vertical atmospheric propagation or in any other scattering scenario that involves changes in the environmental conditions, such as a temperature or pressure gradient. For such a change, the entire system might be seen as cascading of an infinite number of systems of Figure 1b, with gradual changes in their RPS properties. Consider, for example, N adjacent media as described in Figure 4c. Each medium is represented by the case study model in Figure 1b. The RSPs of each medium are Gaussian, with different l tr according to the medium position. Because each medium is represented by two masks, there are n = 2N masks and n + 1 sub-systems between them. The value of l tr for each two masks of the same medium is identical. Considering a gradually linear increase in l tr by a factor of c tr , we obtain l tr_K = Kc tr l tr , where l tr_K denotes the l tr of the K_th medium. Each medium has two scaling parameters that depend on their position in the whole system. Assuming all the media have the same length ∆L, and since we may choose the free parameter f to be ∆L/2, as in Section 2.2, we obtain: 1_K and σ 2 2_K as the variance of the first RPS g k and the second RPS g k+1 of the same K_th medium, accordingly, where k = K/2 and both variances are a function of the same l tr_K . Thus, using Equation (8), the overall AOTF is obtained by Please note that the third equality in Equation (9) is the general AOTF for any N media, each one with its own l tr_K , length L K , and projection parameters. Only the last line is the private case for a gradual linear decrease in l tr by a factor c tr and with the same length L K = ∆L (Figure 4c). The left summation term in the exponent of the last expression could be considered to be a factor for the average effective l tr of the whole system [26]. Please note that for our case study model (Figure 1b), the dependence of the scattering on the wavelength is assimilated into σ by the l tr value, where l tr originates from Mie theory [7] and is fundamentally a function of λ. Thus, the dependence of the AOTF on λ is also obtained through l tr . Figure 4d shows the AOTF of the particular case where the medium in Figure 3a is replaced by the gradual medium change in Figure 4c. The graphs show the cases of two, three, and four cascading media, i.e., N = 2, 3, 4 and n = 4, 6, 8, where n is the number of the masks. For each case of cascaded N media, the sensor in Figure 3a was moved to the new imaging plane. We chose each medium length in Figure 4c to be the same as for the value of Figure 3a. The transport mean free path decreases from medium to medium by a factor of c lt = 0.25. As it can be seen, the graph shows a good match between the simulations and the analytical expression for the multiple masks in Equation (9). Figure 4a-c are examples for multilayer applications of random media with different properties, which might be analyzed by Equation (8). In the analysis of atmospheric light propagation, this may be relevant when the source or the sensor get closer to or further from each other on the same line of sight. In a biological application, this may be relevant, for example, for evaluating the change in the optical propagation through a multi-layered tissue after removing tissue layers during surgery or due to stretching of the tissue.
In the last examples, we use our case study (Figure 1b) with two RSPs, where the need for multiple RSPs is due to the addition of different media to the system or to a change in the medium properties along the path (Figure 4c) of the same model (Figure 1b). However, multiple RSPs may also be needed for describing other scattering models [1][2][3]6]; thus, the APSF of other multiple RSP models might also be described by Equation (8). For example, it was shown that at least two RPSs are needed to represent a turbulent medium (Ch.9 in [27]). Using only two RSPs enables creating a lab simulator for atmospheric light propagation through two spatial light modulators [4]. However, for the investigation of strong turbulence, a large number of RPSs is needed [6]. Similarly, multiple RSPs are needed to analyze the PSF in the transition regime from ballistic to diffusive light transport [3].
Another valuable application of Equation (8) could be to account for arbitrary shapes of the scattering media. For example, consider the case of a random medium with a known mantle embedded in another random medium (Figure 4e). The multislice method (a.k.a. the beam propagation method) models a three-dimensional object as a series of thin two-dimensional slices separated by homogeneous media. The propagation of the light is analyzed from slice to slice [28]. Therefore, we may approximate the two media as a stack of many RPSs with free-space propagation between them [29]. As the number of slices increases and the distance between them decreases, the approximation is better ( Figure  4f). Each slice can be modeled by a combination of three RPSs. One belongs to the inner medium and the others to the environment. If the p of each RPS is known, as well as its statistical AC, then we may use the model. If the scattering of the environment is much smaller than the cell, then we may approximate each slice as a superposition of only two RPSs, one of the environment and the other of the cell.
One of the main technical challenges of the multislice approach is to determine the slice width to be used in the analysis. The common approach is to compute the result with a specific depth resolution, as in Figure 4f, and then, if required, refining by recomputing it with a higher resolution, as depicted in Figure 4g. The refinement process might need to be repeated until convergence to a stable output. In the conventional way, at each refinement step, the output needs to be calculated. However, using our model for computing U out 2 , the effect of each slice might be separated from any other slice. Thus, only the additional slices would need to be computed, and their effect is added to the previous computation by the convolution operator. Therefore, the model is effective in the sense of computational efficiency.

Discussion and Conclusions
The output complex field amplitude of an optical system is formulated in Equation (1) in terms of convolutions of a maskless system with a scaled FT for each mask, with corresponding quadratic phase multiplications. This formulation facilitates separating the effects on the output field of all the masks from the output core system without masks. It also highlights each mask's influence on the output [18]; however, the effect of each mask on the output cannot be entirely separated from the core system due to the quadratic phase multiplications.
In this paper, we analyze the case where the masks are RPSs. For RPSs, only the statistical properties can be evaluated; thus, we examined the average effect. An elegant and comprehensive approach to analyzing the propagation of the first and second order statistics through the system is by examining the Wigner distribution (WD). For example, in [5], we applied such an approach for the system in Figure 1b. The WD analysis [5] of U out cannot reveal the effect of each individual mask due to the non-SI nature of the system. However, the statistical averaging mechanism of APSF eliminates the quadratic phase terms in Equation (1) and separates the statistical effect of each mask. We have developed analytical expressions, and validated with simulations that this mechanism enables separating the statistical effect of each RPS from that of any other RPS and from the core system.
We conclude by outlining a remarkable difference between systems with deterministic masks and systems with RPS masks. For deterministic masks, the main reason for the SI of U out is the masks' positions. This affects the quadratic phase factors between the convolutions in Equation (1), and, even if U ideal is generally a SI system, U out is a non-SI system. In contrast to deterministic mask systems, for the RPS system, the reason for the non-SI feature of the averaged output E |U out | 2 is the non-SI feature of |U ideal | 2 and not the RPS positions.
The presented results are useful for exploring various optical systems with multiplescattering media, and for cases when the random medium properties are not fixed [30] and vary along the system path. The method applies to any paraxial core system and any cascaded system involving different medium properties. For example, the third equality in (9) is valid for any varying function in l tr along the path in Figure 4c, and additional optical systems might be added before, between, and after the media. The method also applies for other models of scattered light in random medium analysis by RSPs, in addition to the specific model in Figure 1b.
There are many different approaches to physically modeling a random medium [1,3,7,22,29]. This paper has made no new contribution to the physical scattering model at the structural level. Rather, the contribution of this paper is at the system model level by separating the RSPs' effects. We generalize the analysis of any paraxial system with deterministic masks to the case of RPSs. Our model applies to a general framework; it can handle any paraxial core system and any cascaded system of sub-systems with RPSs whose ACs are known. Any paraxial system is defined by the linear canonical transform (LCT) [31,32] and the model is valid for any RPS that is a wide-sense stationary (WSS) process. The generalized LCT convolution and the WSS process have many other application areas beside optics [15,16]. Therefore, the results are valid for any cascaded LCT sub-system with multiplication with WSS processes between the LCT sub-systems.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. The Generalized Convolution
Any paraxial system M can be represented by a corresponding ABCD ray transfer matrix M with the parameters a, b, c, d. Let us denote the sub-systems M 1 , M 2 , and M 3 in Figure 2 with the corresponding matrices, M 1 , M 2 , and M 3 , respectively. The matrices can be cascaded, so we will denote the maskless system in Figure 2b Figure 2d is M 3 , with parameters a 3 , b 3 , c 3 , d 3 .
The ray matrix parameters can be used for the evaluation of the field by virtue of the generalized form of the Kirchhoff-Huygens diffraction integral through a general ABCD system, which is also known as LCT [31][32][33][34][35]. The relation between the output field and the input field of a first order optical system is given, up to complex factor, by with the parameters a, b, c, d of its corresponding ABCD matrix M. Thus, we denote the same subscript for O and M. U out is obtained by cascading the integral three times with the parameters of the corresponding matrices M 1 , M 2 , and M 3 , together with multiplication of the masks, i.e., The integral has the same cascading property as the matrices. Hence, as with (Figure 2b), and, similarly, O 3 (O 2 (g 1 )) = O 32 (g 1 ) = U g 1 (Figure 2c USV Symbol Macro(s) The advantage of the generalized convolution is pronounced in the LCT domain [15,36]. However, O(U · g) can also be expressed with regular convolution with a quadratic phase before and after the convolution [16,17]; then, U out is obtained by Equation (1) [18].