Modeling of optical scattering from topographic surface measurements of high-quality mirrors

For simplicity, let us limit the following discussion on smooth surfaces of low-loss and high-reflective mirrors. In other words, we assume the power incident on the surface, $\mathcal{P}_{\mathrm{i}}$, is almost equal to that of the total scattered-out power $\mathcal{P}_{\mathrm{s}}$, which includes specular reflection as well as the other scattered light. Note that, in most of practical cases, it is difficult or rather meaningless to precisely distinguish the specular and scattered components.

Under this condition, the scattering probability density function (PDF) per unit solid angle at a certain point of the illuminated surface is

where $d\Omega_{\mathrm{s}}$ indicates a fraction of solid angle around a scattered angle in concern, which can be expressed by a combination of a polar angle $\theta_{\mathrm{s}}$ and an azimuthal angle $\phi_{\mathrm{s}}$, while $d\mathcal{P}_{\mathrm{s}}$ is the light power scattered within $d\Omega_{\mathrm{s}}$. The angle of incident on the illuminated surface can be also expressed by a similar combination $(\theta_{\mathrm{i}},\phi_{\mathrm{i}})$, but proper choices of the coordinate can set $\phi_{\mathrm{i}}=0$ without losing generality in most of the cases. Note that $\int dP/d\Omega_{\mathrm{s}}\,d\Omega_{\mathrm{s}}=1$ due to this definition or the energy conservation, where the integration is taken over the whole sphere around the concerned point on the illuminated surface. By the way, the traditional BSDF (bi-directional scatter distribution function) corresponds to $(d\mathcal{P}_{\mathrm{s}}/d\Omega_{\mathrm{s}})/(\mathcal{P}_{\mathrm{i}}\cos\theta_{\mathrm{s}})$, which is approximately equal to the scattering PDF when $|\theta_{\mathrm{s}}|\ll 1$.

The traditional Rayleigh-Rice theory relates the topographic height errors of the illuminated surface to the scattering distribution as [13]

where $\lambda$ is a wavelength of the light, $Q$ is a function of $\theta_{\mathrm{s}}$,$\phi_{\mathrm{s}}$, and $\theta_{\mathrm{i}}$ to describe a distribution of the mean light power, and $S_{2}$ is a two-sided 2D (two dimensional) power spectral density (PSD) for a map of the topographic height errors of the surface; see Eq. (6). $S_{2}$ is also a function of the angles, because the spatial frequencies $f_{x}$ and $f_{y}$ are related with the anuglar parameters: $f_{x}=(\sin\theta_{\mathrm{s}}\cos\phi_{\mathrm{s}}-\sin\theta_{\mathrm{i}})/\lambda$, and $f_{y}=(\sin\theta_{\mathrm{s}}\sin\phi_{\mathrm{s}})/\lambda$, respectively.

For the later discussion, here we derive the more reduced forms with some assumptions. Let $\theta_{\mathrm{i}}\simeq 0$. Let us remind the assumption $\mathcal{P}_{\mathrm{s}}\simeq\mathcal{P}_{\mathrm{i}}$. Let most of the scattered light power distribute around the direction of the specular reflection, and so only consider the case of $|\theta_{\mathrm{s}}|\ll 1$. The assumptions are also to set $Q\sim\mathcal{P}_{\mathrm{s}}/\mathcal{P}_{\mathrm{i}}\sim 1$, as $Q$ is a generalized reflectivity in the context of this paper. Then, the scattering PDF becomes $dP/d\Omega_{\mathrm{s}}\simeq 16\pi^{2}S_{2}(f_{x},f_{y})/\lambda^{4}$, where $f_{x}\simeq\theta_{\mathrm{s}}\cos\phi_{\mathrm{s}}/\lambda$ and $f_{y}\simeq\theta_{\mathrm{s}}\sin\phi_{\mathrm{s}}/\lambda$. Furthermore, if $S_{2}$ is circularly symmetric, they reduce to

where

In the following, two widely-used spectral models are introduced; the inverse power law model and the ABC model. Mathematically, the former can be understood as a certain approximation of the latter.

Fig. 1

summarizes our data processing flow. The flow is implemented into a python script. As discussed later, this time the read data are the measured 2D maps of the KAGRA mirrors, which had been characterized in [8]. The map data are saved in ZYGO’s .dat format, and so the prysm library is utilized to read them.

For pre-processing, each map image is cropped to be a 140-mm diameter circle centered at the image’s center. Then, the cropped map is detrended with the poppy library; to put it concretely, components corresponding to the typical low-order Zernike polynomials ($Z_{1}$, $Z_{2}$, $Z_{3}$, and $Z_{4}$ in Table I of [26]) are evaluated, and then subtracted.

Before performing 2D FFT (fast Fourier transform) for the pre-processed map, a 2D Hann window is applied. In the normalized form, it is written as

where $x=\frac{1}{2}+t\cos\phi$ and $y=\frac{1}{2}+t\sin\phi$ for $0\leq t\leq\frac{1}{2}$ and $0\leq\phi<2\pi$, while $w(x,y)=0$ for the other $t$. It takes the maximum 1 at $(x,y)=(\frac{1}{2},\frac{1}{2})$. Note that the domain of the window function should be scaled to the map size in the real coding. The corresponding window is introduced with the scikit-image library, and then multiplied with the map.

After performing 2D FFT for the windowed map, a raw 2D PSD is calculated with

where $Z_{L}^{\diamond}$ is the practical discrete version of $Z_{L}$ in Eq. (5). The window correction factor $c_{w}$ is calculated as $\int_{0}^{1}\int_{0}^{1}w^{2}(x,y)\,dx\,dy$, which corresponds to $2\pi\int_{0}^{1/2}t\,w^{2}(t)\,dt$, or $c_{w}\sim 0.135$. Note that $s_{2}$ is not ensemble averaged yet, so this spectrum has a large estimation error.

One way to obtain $S_{2}^{\diamond}$, which is the practical discrete version of Eq. (6), or the ensemble average of $s_{2}$, is to calculate a circular mean of $s_{2}$ under the assumption that the original mirror map $z(x,y)$ is isotropic. This means the corresponding two-sided 2D PSD should be circularly symmetric. The circular mean can be formally written as

where $f_{x}=f\cos(2\pi k/N)$ and $f_{y}=f\sin(2\pi k/N)$ for $f>0$, and $N$ is a non-negative integer for averaging. Note that $S_{2}^{\diamond}$ is still a two-sided 2D PSD estimate, although $f>0$ is put in the calculation above; in other words, $S_{2}^{\diamond}(f_{x},f_{y})=S_{2}^{\diamond}(f)$.

From $S_{2}^{\diamond}$, relevant parameters in Eq. (11) or (13) are estimated by curve fitting. This way, the two-sided 2D PSD for the measured mirror map is modeled by a simple function with a few parameters. With them, the scattering PDF estimate can be obtained through Eq. (3).

In Fig. 1, there is another branch from the “Raw 2D PSD” node. This path compresses the 2D PSD to a 1D PSD using Eq. (7). Under the assumption that the 2D PSD is circularly symmetric, this corresponds to the Abel transform shown in Eq. (8). Furthermore, the compression automatically results in an averaging effect at the same time. This way, $S_{1}^{\diamond}(f_{x})$, which is a practical discrete version of Eq. (7), is obtained from $s_{2}(f_{x},f_{y})$. Note that $S_{1}^{\diamond}(f_{y})$ can also be obtained by integration on $f_{x}$. In our case, the average $(S_{1}^{\diamond}(f_{x})+S_{1}^{\diamond}(f_{y}))/2$ was adopted for the one-sided 1D PSD estimate. From $S_{1}^{\diamond}$ data, relevant parameters in Eq. (10) or (12) are estimated by curve fitting. These estimated parameters can be compared with the relevant parameters estimated from $S_{2}^{\diamond}$ described above, for consistency check.

This subsection discusses some outcomes when the data processing in Fig. 1 is applied to the measured 2D topographic surface maps of the KAGRA core mirrors as a case study.

Fig. 2 shows the outcome of the “Detrend” node.

There are four core mirrors in the KAGRA interferometer, where their names are ITMX, ITMY, ETMX, and ETMY, respectively⁴⁴4They are acronyms for input test mass x or y or end test mass x or y.; the details for their measurement and characterization are found in [8]. In fact, the method to obtain Fig. 2 is conceptually the same as the one to obtain the relevant figure in [8]. In Fig. 2, each map image is cropped to a 140-mm diameter centered at the mirror center, and the corresponding image data is 351 $\times$ 351 pixels.

Fig. 3

shows the outcome of the “Raw 2D PSD” node, or $s_{2}$ in Eq. (15). For each mirror, the corresponding $s_{2}$ is plot with respect to the spatial frequencies $f_{x}$ and $f_{y}$. The color bar represents the 2D PSDs in terms of the exponent of 10, and the color scale is common for every subplot. The circular symmetry of each $s_{2}$ is apparent.

Outcomes of the “Circ. mean” node, or $S_{2}^{\diamond}$ in Eq.(16) for the four mirrors are plot in Fig. 4

with filled circles. Again note that each 2D PSD is in two-sided form. Then, the worst values among the four data sets at each special frequency are marked with upside-down triangles in the range above $0.03\,\mathrm{mm^{-1}}$.

To perform curve fit, the lower limit of the spatial frequency is set. This threshold is derived as follows. For the KAGRA’s core mirrors, scattered light ($\lambda=1064\,\mathrm{nm}$) to be considered is that spreads more than $110\,\mathrm{mm}/3\,\mathrm{km}\sim 37\,\mathrm{\mu rad}$ in the polar angle [12], where the length of the main Fabry-Pérot arm cavity is 3 km, while each core mirrors at the both ends of the cavity is 110 mm in radius. This corresponds to $\sim 0.035\,\mathrm{mm^{-1}}$ according to Eq. (4). We use 2D PSD values above $\sim 0.03\,\mathrm{mm^{-1}}$ for the curve fit.

The fit curve is drawn with a solid curve in Fig. 4. The fit data are indicated with upside-down triangles, which are the worst-case values among the four mirrors, as already mentioned. This is because it is enough to discuss on the worst case or the upper limit of the noise due to the scattered light for our purpose. Let us remind that our objective is to obtain an estimate of scattering profile that is available for various ray-trace simulation in the future, but not to characterize each mirror’s quality. The real world is quite complicated for this kind of simulation, so certain simplifications and assumptions are unavoidable to obtain usable results efficiently. For example, one of the apparent differences is that rays are described in geometrical optics while the laser light in use is well described in wave optics. Then, the results include various uncertainty. Practically, we further include a safety factor of 10 to 100 (depending on the situation) into the baffle design in the end. Note that, despite these treatments, the currently operated gravitational-wave detectors such as LIGO, Virgo, and KAGRA are still suffered from the stray-light noise through unforeseen coupling paths.

In this way, outcomes of the left “Curve fit” node, or fit parameters can be obtained as $A^{\prime}=2.049\pm 0.141$, $B=15.00\pm 1.22$, and $C=2.198\pm 0.132$ in Eq. (13). Note that $A^{\prime}$ and $1/B$ have the same dimension as of the 2D PSD ($\mathrm{nm^{2}\,mm^{2}}$) and the special frequency ($\mathrm{mm^{-1}}$), respectively, while $C$ is dimensionless. In the actual coding, the curve fit is weighted by the spatial frequency to improve the estimation errors.

The dotted line in Fig. 4 is calculated from the fit curve described above as a power-law limit; $(Bf)^{2}\gg 1$ in Eq. (13). This corresponds to Eq. (11). For the safer upper limit of the scattering profile, the inverse power law estimation may be also useful in some situation.

Outcomes of the “1D PSD” node, or $(S_{1}^{\diamond}(f_{x})+S_{1}^{\diamond}(f_{y}))/2$, are plot in Fig. 5

for ITMX, ITMY, ETMX, and ETMY, respectively. Note that the vertical axis is the one-sided 1D PSD. In the same manner as the 2D PSD case, the worst-case values above the threshold spatial frequency are marked with upside-down triangles. The curve fit is performed for these data with Eq. (12).

In the same manner as the 2D PSD case, outcomes of the right “Curve fit” node are obtained as $A=0.479\pm 0.031$, $B=14.04\pm 0.69$, and $C=2.234\pm 0.043$ for Eq. (12). The curve fit is weighted by the spatial frequency to improve the estimation errors. From the fit result, the relevant $S_{2}(f)$ in Eq. (13) can be also reconstructed. For cross check, this is drawn in Fig. 4 by a dashed line. Compared with the solid line, the overlap with the dashed line is not bad. This fact supports the assumption that the original two-sided 2D PSD is of circular symmetry (also supported by Fig. 3), and means the scattering profile is nearly constant with respect to the azimuthal angles.

As an outcome of “Scattering profile” node, the scattering PDF $dP/d\Omega_{\mathrm{s}}$ can be calculated. To show this, the right vertical axis in Fig. 4 is converted to this scale according to Eq. (3). Also, the scattering angle (polar angle $\theta_{\mathrm{s}}$) that corresponds to the spatial frequency is shown in the upper horizontal axis. Note that the scattering PDF here is of the isotropic surfaces, so should be invariant with the azimuthal angle $\phi_{\mathrm{s}}$ changes.