Backscattering Estimation of a Tilted Spherical Cap for Different Kinds of Optical Scattering

: In many optical engineering applications, a spherical cap shaped optical element is widely used such as concave or convex mirrors in reﬂective optics. Such an element can also tilt around the vertex which corresponds to an off-axis optical design. The optical backscattering of such an optical element sometimes could be important. For example, in the space-based gravitational wave detection, the backscattering of such an element could be superimposed with the local oscillator and limits the sensitivity of the spacecraft. The scattered contributions depend on the scattering property of the mirror surfaces and the geometrical arrangement including the radius of curvature, the tilt and the interval between the scattering source and detector plane. Based on random estimation method, this paper starts from the radiometry, combines these variables and calculates the theoretical amount of back scattered light for both diffuse and superpolished surfaces. The results are compared with analytical and ray tracing solution. The conclusions can be used to further improve the optical design of the telescope or extended to other cases where the backscattered light should be controlled.


Introduction
Light scattering can be considered as an unwanted source that causes image degradation, output reduction, polarization jitter or phase disturbance. For an image system, the optical surface scattering is one of the critical sources that results in stray light. the irradiance of scattered light that received by a detector is related to the scatter properties of the surface and the geometry of the optical system [1]. The scattering properties of an optical surface can be related to imperfections of the surface such as roughness [2,3], contaminations [4,5], subsurface damage [6] etc.al. These imperfections are inevitable even for high quality system and limits the final performance of an optical system.
In real optical systems, the geometrical effect also plays an important role and relies highly on the optical configuration. However, compared with intrinsic scattering, the geometrical effect can be controlled and hence the expected stray light performance of an optical system can be achieved through rational arrangement of the optical elements. For example, an off-axis design, together with the field stop and the Lyot stop, could exhibit a lower level of stray light than an on-axis design [7]. Furthermore, Gary L. Peterson [8] derived the in-field scattered light contribution of a single element in the optical paths and James E. Harvey [9] applies the calculation into a two-mirror system. It shows that the contribution of the stray light of a component is affected by the beam radius. These methods assume that the in-field scattering direction is the same with the direction of the specular rays which is called forward scattering. However, in some special arrangements, the scattering path is different from specular rays where scattering is backward. For example, in LISA

Definitions
The scattering from an optical reflective surface can be described by the bidirectional reflectance distribution function (BRDF) which is the differential radiance (L s ) of a certain scattering surface over the differential irradiance E i of the incident surface [17]: BRDF(θ i , θ s , ϕ s ) = dL s (θ i , θ s , ϕ s ) dE i (θ i ) ≈ P s /Ω s P i cos θ s (1) where P s , P i is the scattered power and the incident power. Ω s is the solid angle of the scattered beam. θ i , θ s , ϕ s is the incident angle, scattering angle and azimuth angle. The unit of BRDF is 1/Sr. The relevant geometry is shown in Figure 1. For an isotropic surface, BRDF exhibits rotational symmetry and can be reduced to an in-plane scattering [18]. The radiometric definition of the scattering is straightforward. However, the accurate measurement and modeling can be difficult. The scattering signature is a result of the roughness, the surface structure, contaminations, scratches, the wavelength and the polarization. Furthermore, the digital noise, optical alignment also makes the signature even more complicated. In this paper, the scattering of the spherical cap is assumed to be roughness. There are a lot of existing models to simulate the scattering of roughness scattering such as Harvey model [19], K-correction [20], ABG [21] and so on. Detailed comparisons The radiometric definition of the scattering is straightforward. However, the accurate measurement and modeling can be difficult. The scattering signature is a result of the roughness, the surface structure, contaminations, scratches, the wavelength and the polarization. Furthermore, the digital noise, optical alignment also makes the signature even more complicated. In this paper, the scattering of the spherical cap is assumed to be roughness. There are a lot of existing models to simulate the scattering of roughness scattering such as Harvey model [19], K-correction [20], ABG [21] and so on. Detailed comparisons of these models are not the goal of this paper. The model used in this paper is the two-parameter Harvey model.
The total integrate scattering (TIS) is defined as the ratio of total power scattered by a surface in the reflected or transmitted direction to the incident power. For an opaque surface, TIS is defined as: For high quality isotropic smooth surfaces, the TIS can be used to estimate the rootmean-square roughness which is: On the other hand, when BRDF is independent from the incident angle or scattering angle, BRDF is a constant and the scattering type is Lambertian. Then the ratio of TIS to BRDF is π. These relations are later used to verify the modeling.
As for the radiation transfer, referring to two randomly oriented surfaces in Figure 2, the differential flux d 2 Φ emitted from dA 1 with radiance L and received by dA 2 is determined by: where r is the distance between the two infinitesimal surfaces, θ 1 , θ 2 is the angle of the scattered ray with respect to the surface normal vector → n 1 and → n 2 . When finite surfaces are involved, the total transferred flux can be derived by integration over the total area A 1 and A 2 : BRDF · E in cos(θ 1 ) cos(θ 2 ) r 2 dA 1 dA 2 (5) Optics 2022, 3, FOR PEER REVIEW 4 As can be seen, the total received flux depends highly on the geometrical layout including the shape, distance, orientation and the area of the surfaces. When the scattering source is Lambertian, the radiance term L can be extracted from the integral and the rest divided by the area is often called configuration factor. Some common geometries with different methods can be found in [22]. In the following section, the surface 1 is considered to be a tilted spherical cap and surface 2 is assumed to be a circular planar detector. For simplicity, the view edge effect is not considered which can be found in [23].
When focusing more on the geometrical flux emission characteristics of the emitter itself, the radiance term in Equation (5) can be ignored and the rest is often called Étendue or optical throughput which is defined as [24,25]: As can be seen, the total received flux depends highly on the geometrical layout including the shape, distance, orientation and the area of the surfaces. When the scattering source is Lambertian, the radiance term L can be extracted from the integral and the rest divided by the area is often called configuration factor. Some common geometries with different methods can be found in [22]. In the following section, the surface 1 is considered to be a tilted spherical cap and surface 2 is assumed to be a circular planar detector. For simplicity, the view edge effect is not considered which can be found in [23].
When focusing more on the geometrical flux emission characteristics of the emitter itself, the radiance term in Equation (5) can be ignored and the rest is often called Étendue or optical throughput which is defined as [24,25]: Étendue can be used to determine the theoretical flux transfer of an optical system and is often used in non-imaging or illumination design.

Mathmetical Modeling
The configuration of the system is shown in Figure 3. Relevant symbols can be found in Table 1. The optical axis is along Z axis. The detector is a disk with radius r 1 on XY plane and the polar coordinate of the infinitesimal patch of the detector is ρ2(0 < ρ2 < r 2 ), ϕ 2 (0 < ϕ 2 < 2π). Then the height of the vertex z 0 is determined by: Optics 2022, 3, FOR PEER REVIEW 5  x y z The center of spherical cap before rotation , ,

O2
O2 O2 x y z The center of spherical cap after rotation , , 0 0 0 x y z Cartesian coordinates of the sample before rotation , , Cylinder coordinates of the sample before rotation , , 1 1 1 x y z Cartesian coordinates of the sample after rotation , , 2 2 x y 0 Cartesian coordinates of the detector r Distance between two differential patches. A1, A2 Area of spherical cap and the detector dA1, dA2 The differential area of A1, A2 The radiation source is in the shape of a spherical cap and the center, with a height  The center of spherical cap before rotation x O2 , y O2 , z O2 The center of spherical cap after rotation x 0 , y 0 , z 0 Cartesian coordinates of the sample before rotation ρ 0 , ϕ 0 , z 0 Cylinder coordinates of the sample before rotation x 1 , y 1 , z 1 Cartesian coordinates of the sample after rotation x 2 , y 2 , 0 Cartesian coordinates of the detector r Distance between two differential patches. A 1 , A 2 Area of spherical cap and the detector dA 1 , dA 2 The differential area of A 1 , A 2 The radiation source is in the shape of a spherical cap and the center, with a height of d radius R, is on the Z axis. The spherical center O 1 is (0, 0, d − R). The mirror rotates around X axis with an angle of θ which is common for an off-axis optical layout. After rotation, the spherical center becomes O 2 with the coordinate: For an arbitrary point A 0 , with the polar coordinate ρ 0 (0 < ρ 0 < r 1 ), ϕ 0 (0 < ϕ 0 < 2π), the point after rotation becomes A and the Cartesian coordinates x 1 , y 1 , z 1 are: Then, the distance r is: For an arbitrary point B on the detector, the normal vector is along Z-axis. The normal vector of point A on the spherical cap is along vector → AO 2 . Therefore cos(θ 1 ) and cos(θ 2 ) can be calculated as: After plugging Equations (7)−(11) together with dA 1 , dA 2 into Equation (5), the total flux emitted from the spherical cap, received by the detector can be derived through the integration. In order to solve the high dimensional integration, a simple feasible method is based on random estimation (Monte Carlo method) where N random pairs of patches A i, B i are introduced. The basic idea is to randomly sample detector and the spherical cap and the vectors from A i to B i are assumed to be the scattered rays. The Cartesian coordinate of A i is (x 1i , y 1i , z 1i ) which is a random patch on the spherical cap whereas B i, with the coordinate of (x 2i , y 2i ), is a random patch on the detector. Then total transferred flux become: The Monte-Carlo method has the following advantages. First, the standard deviation is proportional to 1/ √ N regardless the dimension and smoothness of the integrand [26]. Therefore, it is useful to solve high dimension integrals, or the integrand is not continuous. Furthermore, Equation (12) directly samples the positions of the differential patches which maximizes the collecting efficiency compared with sampling direction aimlessly. The Monte-Carlo method also exhibits a great flexibility especially when a visibility function is introduced. For example, a Circ function can be inserted in Equation (12) which means that the scattered rays are selected by an aperture.
In the following parts, Equation (12) is used to calculate the backscattering of the quaternary mirror of the telescope for the Taiji GW detection program. The space-borne instrument can be divided into two parts. The first part is the optical bench contains an interferometer and a Proof Mass. More details can be found in Ref. [10]. The other part is a four mirror off-axis telescope which transmits the laser from interferometer and received the laser from the far-end telescope simultaneously. The telescope delivers a gaussian laser beam (wavelength λ = 1064 nm, transmitted power P tot = 2 W, beam waist w = 1.875 mm) to the far-end telescope with a 400 mm entrance pupil. The received power for the receiving end is around 700 pw. The phase noise φ bsc (t) introduced by the backscattered light is [13,27]: where l 0 is the static optical path length while δl dis (t) is the displacement noise from the scattering surfaces. φ static and δφ dis (t) are the corresponding phase. Then the backscattered field E bsc produces an amplitude and phase modulation to the undisturbed interferometer field E in : where P in is the power inside the interferometer at a reference point for calculation. f r is the fraction of the transmitted power which is backscattered by the telescope and recoupled with the receiving beam. Since P tot , P in are fixed, f r in Equation (14) should be kept small to reduce the impact of backscattering. In order to keep the phase noise at picometer level, the fraction f r is set to at a level of 10 −10 according to the specification of the current Taiji program. Moreover, due to the high coherence of the laser, the impact scattered light can be described by additional phase or amplitude. The additional phase noise can be suppressed by frequency stabilization [28] and then the amplitude noise is dominant [13].
The local optical layout of the tertiary and the quaternary mirror is shown in Figure 4. The scattering from the tertiary mirror is temporarily not considered. The quaternary mirror, which directly faces to the entrance of the optical bench, delivers the scattering from the other surfaces and plays a major part in the backscattering. So, in the design phase, it is particular critical to accurate estimate the scattering contributions of the quaternary mirror. Assuming the telescope is fed by a fundamental mode gaussian beam with total power P tot, Beam waist of w, the irradiance at a certain position of the mirror is: Optics 2022, 3, FOR PEER REVIEW 7 backscattered field bsc E produces an amplitude and phase modulation to the undisturbed interferometer field in E : where in P is the power inside the interferometer at a reference point for calculation. r f is the fraction of the transmitted power which is backscattered by the telescope and recoupled with the receiving beam. Since , tot in P P are fixed, r f in Equation (14) should be kept small to reduce the impact of backscattering. In order to keep the phase noise at picometer level, the fraction r f is set to at a level of 10 10  according to the specification of the current Taiji program. Moreover, due to the high coherence of the laser, the impact scattered light can be described by additional phase or amplitude. The additional phase noise can be suppressed by frequency stabilization [28] and then the amplitude noise is dominant [13].
The local optical layout of the tertiary and the quaternary mirror is shown in Figure  4. The scattering from the tertiary mirror is temporarily not considered. The quaternary mirror, which directly faces to the entrance of the optical bench, delivers the scattering from the other surfaces and plays a major part in the backscattering. So, in the design phase, it is particular critical to accurate estimate the scattering contributions of the quaternary mirror. Assuming the telescope is fed by a fundamental mode gaussian beam with total power Ptot, Beam waist of w, the irradiance at a certain position of the mirror is: Plugging Equation (15) into Equation (12), the total flux of the backscattered light can be estimated by: In Equation (16), BRDFi is a function of wavelength, incident angle i θ and scattering angle s θ . Since the incident gaussian beam is along Z-axis and the considered scattering is backward, the scattering angle is minus and hence the angle of incidence and the scattering angle fulfill: Plugging Equation (15) into Equation (12), the total flux of the backscattered light can be estimated by: In Equation (16), BRDF i is a function of wavelength, incident angle θ i and scattering angle θ s . Since the incident gaussian beam is along Z-axis and the considered scattering is backward, the scattering angle is minus and hence the angle of incidence and the scattering angle fulfill:

Results
In this section, the results of the scattering flux of a spherical cap are presented. In this paper, two kinds of samples are considered, Section 3.1 demonstrates the results for the diffuse sample under different configurations and the results are compared with existing models. Section 3.2 shows the results of the backscattering of a quaternary mirror of the GW telescope which is superpolished. The scattering model here used is the two-parameter Harvey model. In the model, the parameters of b and s are introduced where b means the BRDF value at 0.573 • from the direction of specular direction and s means the slope of BRDF on the log-log plot. Section 3.3 demonstrates the relevant error analysis.

Diffuse Samples
A perfect diffuse surface scatters the incident beam equally in all directions meaning that BRDF is independent from incident angle and the scattering angle. For a diffuse and homogeneous incident beam, the scattered radiance is a constant and can be extracted from the integration. Then Equation (5) can be simplified to: where M is the exitance and F is usually called configuration factor, which is: For such a surface, the geometrical effect is predominant and the configuration factor denotes the fraction of energy emitted or reflected and is collected by the detector. The determination of the configuration factor under different layouts is well-established. Common ways are:

2.
Statistical method, since the scattering is independent from direction, the configuration factor approximates the number of rays received by the detector divided by the number of rays emitted from the scattering source.

3.
Algebraic method, the laws of closeness, reciprocity, distribution and composition can be applied to derive the target configuration factor. 4.
Projection method, under some special layouts, the target surface can be projected on a unit sphere in order to simplify the calculation.
For a spherical cap, by considering the tilt angle, shape, size and the distance, a rigorous analytical integral can be developed however it is complicated to solve. With the numerical methods, an approximate solution can be applied to the problem and exhibits considerable freedom for arbitrary arrangements. In the following, the calculated configuration factors at some special cases are compared with the existing analytical solutions. The results of the diffuse spherical cap are shown in Figure 5 and the number of random patch pairs are 1500. In Figure 5a, R = 330 mm, r 1 = √ 2/2 · R, d ranges from R − R 2 − r 2 1 to R, ρ 2 ranges from 1 to r 1 . When the size of the detector is much smaller than d, the detector can be considered to be infinitesimal and the configuration is calculated though: where sinα = r1/R Detailed derivations with Stokes theorem can be found in Appendix A.
The result is plotted as the red solid line in Figure 5a. Another interesting case is that d becomes close to R − R 2 − r 2 1 . In such a case, the detector plane is on the projection plane of the spherical cap. Since all the radiation emitted from the detector can be collected by spherical cap, the configuration factors can be easily derived with reciprocity relation. The analytical configuration factor for such a case is A 2 /A 1 and when ρ 2 = r 1 , the configuration factor is cos 2 (α/2). The results are also plotted as the blue solid line in Figure 5a. The analytical results are in good agreement with the analytical solutions. For larger r 2 , it shows more relative errors. The reason is straightforward: the standard deviation is inverse proportional to square root of N. Therefore, a larger area needs more sampling pairs which is one of the key parameters for accurate Monte Carlo estimation. More details will be explained in Section 3.3.
where / sin α r1 R  . Detailed derivations with Stokes theorem can be found in Appendix A. The result is plotted as the red solid line in Figure 5a. Another interesting case is that d becomes close to 2 2 1 R R r   . In such a case, the detector plane is on the projection plane of the spherical cap. Since all the radiation emitted from the detector can be collected by spherical cap, the configuration factors can be easily derived with reciprocity relation. The analytical configuration factor for such a case is / 2 1 A A and when 2 1 ρ r  , the configuration factor is cos ( / ) 2 α 2 . The results are also plotted as the blue solid line in Figure 5a.
The analytical results are in good agreement with the analytical solutions. For larger r2, it shows more relative errors. The reason is straightforward: the standard deviation is inverse proportional to square root of N. Therefore, a larger area needs more sampling pairs which is one of the key parameters for accurate Monte Carlo estimation. More details will be explained in Section 3.3.
(a) (b) The other case is the configuration factor changes with the radius of curvature and tilt angle. In the optical layout, r2 = 2.5 mm, d = 118 mm R ranges 330 mm to 5000 mm, θ from 0° to 30°. The results are shown in Figure 5b where the grid is for the calculated results and the blue curve is for the analytical results. Full analytical expressions can be found in Appendix A. The figure shows that as R increasing, the spherical cap becomes close to a disc and the result for this case fulfills the following well-known relation: Relevant quantities can be also found in [29]. On the other hand, it can be seen that the scattering from a spherical cap has small dependence on the orientation which meets the assumption of Lambertian scattering. The other case is the configuration factor changes with the radius of curvature and tilt angle. In the optical layout, r2 = 2.5 mm, d = 118 mm R ranges 330 mm to 5000 mm, θ from 0 • to 30 • . The results are shown in Figure 5b where the grid is for the calculated results and the blue curve is for the analytical results. Full analytical expressions can be found in Appendix A. The figure shows that as R increasing, the spherical cap becomes close to a disc and the result for this case fulfills the following well-known relation: Relevant quantities can be also found in [29]. On the other hand, it can be seen that the scattering from a spherical cap has small dependence on the orientation which meets the assumption of Lambertian scattering.

Backscattering of the Quaternary Mirror of Taiji GW Telescope
In this section, a superpolished spherical cap is considered. Compared with the diffuse sample, the most obvious difference is that the scattering is directional. On the other hand, in practice the profile of the incidence irradiance is also not necessary to be uniform. Therefore, the radiance term in Equation (5) has to be kept in the integral and a quite different conclusion can be drawn. This section deals with such a case and in order to calculate the backscattering of the quaternary mirror of Taiji GW telescope. For simplicity, the total power of the incident gaussian beam is set to be unit which is easy to calculate the fraction f r . The diameter of the mirror is 8 mm and the diameter of the collector is 5 mm. The scattering distributions of the mirror are simulated with Harvey model. Figure 6 is the illustration of BRDF under different angle of incidence where b = 0.001, s = −2. With the considered configuration, TIS is 17.9 ppm which corresponds to σ rms ≈ 0.36nm. In practice, these parameters should be modified according to the real measured data. Then, the flux of the backscattered light can be calculated as a function of scattering distribution, tilt angle and radius. The calculated results are compared with commercial soft ASAP [30] which is famous for the outstanding capability, flexibility, speed and accuracy. The scattering calculation of ASAP is based on ray tracing which additionally generates a set of scattering direction cosines in the importance edges and the energy is not always conserved [31].

Backscattering of the Quaternary Mirror of Taiji GW Telescope
In this section, a superpolished spherical cap is considered. Compared with the diffuse sample, the most obvious difference is that the scattering is directional. On the other hand, in practice the profile of the incidence irradiance is also not necessary to be uniform. Therefore, the radiance term in Equation (5) has to be kept in the integral and a quite different conclusion can be drawn. This section deals with such a case and in order to calculate the backscattering of the quaternary mirror of Taiji GW telescope. For simplicity, the total power of the incident gaussian beam is set to be unit which is easy to calculate the fraction r f . The diameter of the mirror is 8 mm and the diameter of the collector is 5 mm. The scattering distributions of the mirror are simulated with Harvey model. Figure 6 is the illustration of BRDF under different angle of incidence where b = 0.001, s = −2. With the considered configuration, TIS is 17.9 ppm which corresponds to .
In practice, these parameters should be modified according to the real measured data. Then, the flux of the backscattered light can be calculated as a function of scattering distribution, tilt angle and radius. The calculated results are compared with commercial soft ASAP [30] which is famous for the outstanding capability, flexibility, speed and accuracy. The scattering calculation of ASAP is based on ray tracing which additionally generates a set of scattering direction cosines in the importance edges and the energy is not always conserved [31].  Figure 7a where the mirror has tilt angle of 5°. For Monte-Carlo estimation, 10,000 random pairs of patches are chosen and the backscattered flux is plotted the blue curve and the gray curve is the corresponding standard deviation. The results are in good agreement with ASAP calculation (Figure 7a   . The flux increases dramatically as the distance decreases which is because the detector goes into the specular reflected paths which can be seen in Figure 4. On the other hand, for small distance, a larger variance can be witnessed which can be explained by the geometrical effects. As mentioned, BRDF is a strong function of angle of incidence and scattering angle. When the distance is small, the normal vectors of the spherical cap different positions are different meaning that the local patches scatter the incident beam differently. However, as the distance increases and when the area of the spherical cap is much smaller than the distance, the geometrical shape of the scattering source has little effect on scattering. This assumption is commonly used in remote sensing.  Figure 7a where the mirror has tilt angle of 5 • . For Monte-Carlo estimation, 10,000 random pairs of patches are chosen and the backscattered flux is plotted the blue curve and the gray curve is the corresponding standard deviation. The results are in good agreement with ASAP calculation (Figure 7a red curve) for large distance. The backscattering flux of Monte-Carlo methods at d = 118 mm is 3.09 × 10 −10 W whereas the ASAP result is 4.5134 × 10 −10 W. The flux increases dramatically as the distance decreases which is because the detector goes into the specular reflected paths which can be seen in Figure 4. On the other hand, for small distance, a larger variance can be witnessed which can be explained by the geometrical effects. As mentioned, BRDF is a strong function of angle of incidence and scattering angle. When the distance is small, the normal vectors of the spherical cap different positions are different meaning that the local patches scatter the incident beam differently. However, as the distance increases and when the area of the spherical cap is much smaller than the distance, the geometrical shape of the scattering source has little effect on scattering. This assumption is commonly used in remote sensing. In some literature, the geometrical effect is simplified to the solid angle when d is much larger than the size of scattering source. Figure 7b shows the evolvement of the average of incidence angle and the scattering angle from small distance to large distance. The error bars represent the standard deviation of the incident and scattering angles at different positions. For the considered case, the angle of incidence (Blue curve in Figure 7b) is near a constant of 5 • which corresponds to the tilt angle of the spherical cap. The scattering angle (Red curve in Figure 7b) are minus meaning the backward scattering and both the average and the variance for small distance are larger than that of larger distance.

The calculated backscattering flux changes with distance d is shown in
In some literature, the geometrical effect is simplified to the solid angle when d is much larger than the size of scattering source. Figure 7b shows the evolvement of the average of incidence angle and the scattering angle from small distance to large distance. The error bars represent the standard deviation of the incident and scattering angles at different positions. For the considered case, the angle of incidence (Blue curve in Figure 7b) is near a constant of 5° which corresponds to the tilt angle of the spherical cap. The scattering angle (Red curve in Figure 7b) are minus meaning the backward scattering and both the average and the variance for small distance are larger than that of larger distance. Furthermore, the standard deviation of scattering angle has the same tendency of the standard deviation of the final scattering flux which also explains the reason why large distance calculation is stabler and more accurate. So, threshold distance dthreshold is the position from which the geometrical effect becomes stable. It defined as the distance where the standard deviation product of the scattering angle () i σ d and incident angle () s σ d equals to 1/e of the maximum. In detail, the distance dthreshold fulfills: In Equation (22), the standard deviation of angle of incidence and scattering angle are chosen rather than that of the final backscattering which can minimize the impact of the numerical calculation. For the considered configuration, dthreshold is calculated as 34 mm.
Another important parameter for backscattering control is the tilt angle of the spherical cap. As mentioned, the reflection becomes directional and as a result, the backscattering performance of a superpolished surface relies highly on the tilt angle and hence the performance of backscattering can be controlled. This fact is illustrated in Figure 8  Furthermore, the standard deviation of scattering angle has the same tendency of the standard deviation of the final scattering flux which also explains the reason why large distance calculation is stabler and more accurate. So, threshold distance d threshold is the position from which the geometrical effect becomes stable. It defined as the distance where the standard deviation product of the scattering angle σ i (d) and incident angle σ s (d) equals to 1/e of the maximum. In detail, the distance d threshold fulfills: In Equation (22), the standard deviation of angle of incidence and scattering angle are chosen rather than that of the final backscattering which can minimize the impact of the numerical calculation. For the considered configuration, d threshold is calculated as 34 mm.
Another important parameter for backscattering control is the tilt angle of the spherical cap. As mentioned, the reflection becomes directional and as a result, the backscattering performance of a superpolished surface relies highly on the tilt angle and hence the performance of backscattering can be controlled. This fact is illustrated in Figure 8. where the backscattering flux is calculated as a function of tilt angle and the roughness of the mirror. The blue curves are from the Monte-Carlo method and the red curves are for ASAP. The solid lines are for the calculations whose RMS = 0.36 nm, the dashed lines are for the surface with 0.6 nm and the dotted lines are for RMS = 1.42 nm. The backscattering flux increases as the roughness increases. The Monte-Carlo method and the ASAP calculation show the same tendency of decreasing of the total flux with the increasing of the tilt angle. From the graph, one orders drop of backscattering flux can be seen as the tilt angle from 5 • to 10 • . A larger tilt angle usually means a larger scattering angle. As mentioned before, the scattering angle with respect to the specular ray is approximately equals to twice of the tilt angle for a large distance and radius of curvature. Therefore, the received flux can be controlled by properly arranging the tilt angle of the mirror. flux increases as the roughness increases. The Monte-Carlo method and the ASAP calculation show the same tendency of decreasing of the total flux with the increasing of the tilt angle. From the graph, one orders drop of backscattering flux can be seen as the tilt angle from 5° to 10°. A larger tilt angle usually means a larger scattering angle. As mentioned before, the scattering angle with respect to the specular ray is approximately equals to twice of the tilt angle for a large distance and radius of curvature. Therefore, the received flux can be controlled by properly arranging the tilt angle of the mirror. On the other hand, it is also found that when the optical layout fulfills Equation (22), the standard deviation of the incident angle and the scattering angle are small meaning that the ratio of backscattered flux to BRDF is around a constant which equals to: cos( ) cos( ) ( , ) 12 12 Figure 9 shows the ratio of backscattered flux to BRDF under different roughness, As can be seen, the ratio is a constant around 9.5143 4

10
Sr W −  and is independent from the roughness of the surface. Moreover, when the incident irradiance is homogeneous, Equation (23) can be further related to Étendue and the scattered flux can be easily calculated with: On the other hand, it is also found that when the optical layout fulfills Equation (22), the standard deviation of the incident angle and the scattering angle are small meaning that the ratio of backscattered flux to BRDF is around a constant which equals to: E in (x, y) · cos(θ 1 ) · cos(θ 2 ) r 2 dA 1 dA 2 (23) Figure 9 shows the ratio of backscattered flux to BRDF under different roughness, As can be seen, the ratio is a constant around 9.5143 × 10 −4 W · Sr and is independent from the roughness of the surface. Moreover, when the incident irradiance is homogeneous, Equation (23) can be further related to Étendue and the scattered flux can be easily calculated with: Optics 2022, 3, FOR PEER REVIEW 13 Figure 9. The ratio of backscattering flux to BRDF with respect to the tilt angle for the samples under different roughness.

Error Analysis
The quality of the integral can be increased with the following strategies. The first and the simplest one is to increase the number of samplings. As mentioned before, the standard deviation is inverse proportional to N which makes Monte-Carlo to be a

Error Analysis
The quality of the integral can be increased with the following strategies. The first and the simplest one is to increase the number of samplings. As mentioned before, the standard deviation is inverse proportional to √ N which makes Monte-Carlo to be a practical technique for high dimension problems. However, the square root of the number of samplings may lead to low efficiency and in practice there always exists a convergence rate issue even though the N is large. Importance sampling reduces variance by changing the probability density function. In Equation (12), the probability density function is set to be uniform which is A 1 A 2 . Therefore, variance always exist because of the mismatch between the probability density function and the integrated function. An ideal case is to find a probability density function that is proportional to the function and the variance is zero. However, such a distribution is difficult to find. Since the value of the integral must be known beforehand. A common method is to estimate the probability density function according to a certain known term is of the integrand [26]. Since the function is smooth and continuous for the considered optical layout, the improvement is limited and the probability density function is kept uniform. However, the quality of the randomness is critical and directly affects the accuracy of the results. A pseudorandom process appears to be random but is not. It will cause the mismatch between the samplings and the probability density function in Equation (12) and is fatal to the results. More details about such an effect can be found in [32,33]. Modern methods to generate random points include: N-Rooks [34], Multi-jittered [35] which are aiming at acquiring better randomness during the sampling.
The other practical technique is to control the source of variance. Since the sampling patches are independent, by controlling the source of variance can reduce the overall variance. In this paper, controlling the rays with extreme value is an effective method. After verification, the largest source of variance in Equation (12) is the BRDF term especially when near the specular direction. These rays should be less weighted or resampled in the iteration in order to increase the robustness and the accuracy of the algorithm. More details can be found in the source code.

Discussion
In this paper, a method of calculating the direct scattering flux is raised which is not only suitable for diffuse surfaces but also for the case where the incident irradiance and the scattering distribution of the scattering source are not uniform. The calculation starts from radiometry with the method of random estimation. The basic idea is to choose a larger number of random patch pairs between the scattering surface and detector surface and calculate the radiation transfer between these patches. In this way, the collecting efficiency can be maximized. It is also demonstrated that how the integral of flux calculation is simplified under different situations. The result can be applied to heat transfer, optical scattering or illumination calculation. For diffuse surfaces, the configuration factors are calculated and are in good agreement with the analytical solutions under some specific configurations. As for the scattering flux for a superpolished surface, the theoretical best backscattering performance of the quaternary mirror of off-axis telescope in space GW detection is estimated where the incident beam is gaussian. The results show that the distance between the spherical cap and the detector, tilt angle, and the scattering property are the key parameters in scattering control. With the current design, it is recommended that the root mean square roughness smaller than 0.4 nm and the tilt angle larger than 5 • . It is worth mentioning that these conclusions only consider the stray light effects. For the final optical design, more effects should be taken into consideration and therefore, a compromised value should be chosen. For example, a larger tilt angle usually introduces more off-axis aberration such as coma. However, the results could be a reference in the optical design phase. The relevant error analyses are also given in this paper. It shows that the randomness of the sampling and the controlling of the extreme values are crucial for the algorithm. The random estimation method exhibits great flexibility and accuracy for high dimension problems. By following the structure of the paper, the scattering flux of other geometrical layouts could also be calculated. In addition, based on the approach, the problem of radiation transportation could be further solved which is indirect scattering or multi-scattering. Moreover, with random estimation, more information can be integrated such as the visibility or polarization reaction of the surface.

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