Stray light estimates due to micrometeoroid damage in space optics, application to the LISA telescope

Abstract. The impact on an optical surface by a micrometeoroid gives rise to a specific type of stray light inherent only in the space optical instruments. This causes a double source of light scattering: the impact crater and the ejected contamination. We propose a method of stray light estimation and apply it to the case of the Laser Interferometer Space Antenna telescope. We estimate the backscattering fraction for nominal (4 years) and extended (10 years) mission durations.


Introduction
Scattering due to micrometeoroid damage is a specific type of stray-light, which is inherent only in the space instruments with optics exposed to the space environment. Free flying dust particles in space can hit and damage the optical surface. These will cause an increase of stray-light in the system. The problem of micrometeoroid damage exists since the first space flights. 1 However, only a few papers 2,3 give an estimate of this type of stray-light. This work aims at improving the simulation of scattered stray light that results from the impact of micrometeoroids in space optical instruments.
To estimate light scattering induced by the micrometeoroid damage, we propose a method that consists of four steps: 1. Definition of the environmental conditions (Particulates Environmental Model) of the satellite: estimation of the flux and parameters of the particles that arrive at the critical surfaces. 2. Calculation of the expected damage crater diameter (DCD) and ejected mass due to the micrometeoroid impact. 3. Calculation, with the Peterson model, 4 of the bidirectional scattering distribution function (BSDF) that results from the impact craters. Calculation of the corresponding cleanliness level and slope due to contamination by ejected mass. 4. Optical software (FRED) calculation of the scattered light.
In this paper, we apply this method to the case of the LISA (Laser Interferometer Space Antenna) telescope and consider the sun-orbiting LISA trajectory, 50 Mkm away from the Earth, in the micrometeoroid flux estimates. In the following, each step will be described, and the results will be presented.
1916. 6 The first direct observation of GW was made in 2015 by the ground-based interferometer LIGO 7 in the United States (Nobel Prize 2017). Due to seismic noise, gravity gradient noise, and other effects, ground-based interferometers (LIGO, Virgo) are limited at low frequencies.
Significant improvements can be made by the construction of an underground interferometer (Einstein Telescope) on the stable lithospheric plate. 8 However, even underground detectors will be limited to observing the merger of very compact massive objects, generating a signal at relatively high frequencies. To detect GW at lower frequencies and detect more massive and slower objects, the detector has to be in space and is likely to be a space interferometer mission such as LISA. The LISA frequency band will include GW generated by compact objects captured by supermassive black holes in galactic nuclei, compact binaries in our galaxy and beyond, binary supermassive black holes, and quantum fluctuations in the early universe.
LISA is a GW observatory space mission now in phase A. It is a constellation of three identical satellites, forming an equilateral triangle of arms-length L ¼ 2.5 Mkm. The continuous measurement of the light path length L½1 þ hðtÞ between two test masses reveals the presence of the GW hðtÞ. It is performed via six interferometric heterodyne phase measurements that take place at each end of the triangle arms. In the presence of any light that was not intended to be in the design (stray light), the heterodyne interference can be disturbed. The correct assessment of the stray light, together with the instrument stability, is vital for LISA measurements, which aim at a precision of 10 −21 in the GW hðtÞ strain, that is, in the measurement of the fractional change of the light path length L½1 þ hðtÞ.
To make the phase measurements between beams from distant satellites possible, an optical telescope is used for transmission and receiving of the beams. The schema of the LISA telescope is given in Fig. 1. In the NASA design, 9 it is an afocal Cassegrain telescope, which consists of four mirrors. The mirrors are named in sequential order following the beam from the big entrance aperture: M1 (primary mirror), M2 (secondary mirror), M3, and M4. The current off-axis design provides better performance in terms of diffracted light, in comparison to the on-axis configuration. To provide excellent thermal stability of the telescope, the material of the M1 mirror is chosen to be a Zerodur ðRÞ with a thin silver coating. Nevertheless, the approach developed below is of general applicability and is not specific of this brittle material.
The telescope is a part of the moving optical subassembly (MOSA), which also includes the optical bench and the gravitational reference sensor. From the optical bench, the transmitted beam propagates through the telescope and is sent to the distant satellite. The same telescope is used for collecting the received beam. On the optical bench, part of the transmitted beam is pinched off and used for the interference with the received beam. The phase of the interferometric signal encodes information about the GW hðtÞ.
A thermal shield will surround the mirrors of the telescope, and only mirror M1 will be exposed to the space environment.

Stray Light Problem in the LISA Telescope
In the presence of stray light, the phase measurements can be disturbed. The stray light can be considered as an additional beam if it has the same polarization and optical frequency as the reference beam. Stray light in LISA includes exterior contributions (stars, planets, etc.), Fig. 1 Schematic of the Cassegrain four mirror telescope designed by NASA for the LISA mission. The same telescope is used to expand the transmitted beam and collect the light of the received beam.
diffraction of the beam, multiple reflections in optics (ghosts), and scattering of the nominal beams. The sources of scattering are: 1. Surface microroughness 2. Contamination: particle, molecular, biological on the surface and in the volume 3. Cosmetic defects: digs, scratches either due to the assembling of the optical components or resulting from the impact of micrometeoroids 4. Backscattering in the optical fibers A common way to describe scattering is the usage of the BSDF, which designates how the scattering probability depends on the scattering angles. BSDF is the radiance of the scattering surface normalized by the irradiance of the surface. 10 BSDF is used by optical software (FRED, 11 Zemax, 12 Code V 13 ) to perform stray light simulations. When optical and mechanical geometries are defined and optical properties set, then a BSDF is assigned to each surface. Later, each stray light path is identified. To avoid stray light, the optical designs generally use smooth, clean optics (mirrors, beam dumps), black coating for support structures, unique optical designs, baffles, and stops, etc.

Environmental Conditions: Particulates Environmental Model
The first step of stray light analysis due to micrometeoroid damage is to determine the particulate environment for the satellite. This includes information about the flux, velocity, density, mass, and directivity of the micrometeoroids. In the particular case of LISA, the environment of the satellite is given in the LISA Environment Specification document. 14 Charter 5 of that document contains information about the micrometeoroid distribution. However, below we discuss a general approach to the solution of the problem.
The flux-mass model for meteoroids at one astronomical unit from the Sun has been proposed and presented by Grün et al. 15 (see Fig. 2). It gives the total average meteoroid flux ϕ G ðmÞ (sporadic+ stream average) in terms of the integral flux (i.e., the number of particles per square meter, per year, of mass larger than or equal to a given mass m, impacting a randomly oriented flat plate under a viewing angle of 2π). Except for Earth shielding and gravitational effects (which are negligible at the LISA altitude of 50 Mkm), this flux is omnidirectional. 16 The Grün model accounts for averaging over all directions. Besides this, micrometeoroid streams are not considered. This interplanetary flux is valid for the micrometeoroid mass range of 1 × 10 −18 g to 1 g.
For the impact crater calculation, we use a density value of 2.5 g∕cm 3 for all ranges of mass, as was specified in the LISA Environmental model. 14 However, since the mass density of the micrometeoroids is not a measured quantity, a reasonable assumption should be used instead. Another description of the density is given as a function of steps 16 (see Fig. 2). This function will be used in Sec. 3.4.  Fig. 2 The Grün meteoroid flux-mass model (in black) and the step function of density distribution (in blue). 16 As a first approximation, we use the constant value of meteoroid impact velocity 20 km∕s for all the masses of micrometeoroids. This value is the typical mean velocity for a micrometeoroid impact with a sun-orbiting body, 17 and it is proposed by LISA environmental model. 14 The space debris environment is not considered. The next step is to apply this model to the investigated surface. For this, some parameters of the mission are required: duration of the mission, nature of the critical surfaces (size, material, orientation, shielding), etc. For the calculation of the number of expected impacts, we use the Grün model 14 and the following parameters: • Nominal mission duration is 4 years (extended duration is 10 years). • The primary M1 mirror diameter of the LISA telescope is 0.3 m. We assume that M1 is the only mirror exposed to micrometeoroids.
Because the mechanical structure is unknown at the moment, no correction is made to account for the shielding by the structure surrounding the telescope. The approach presented in this paper corresponds to this worst-case scenario.
Since the flux ϕ G ðmÞ is a cumulative flux, to know the number of expected impacts in a certain mass range (one bin), the difference between neighboring bins should be taken into account. The result of the calculation is given in Fig. 3. The bin size is uniform in the logarithmic scale. The ratio between neighboring bin size is 10 1∕10 . The number of expected impacts is a fractional value (different from integer 2 ), as it is a statistical quantity. In the mass spectrum of the meteoroid (see Fig. 3), we neglect the high mass tail, for which the cumulated flux is lower than 1∕e 2 .

Effect of the Micrometeoroids
The hypervelocity impact of the optical surface causes a double effect in terms of scattering. The direct result of the impact is a microcrater. It causes scattering inside of it and diffraction on the border of the crater. From studies of lunar craters and craters in hardware returned from space, it is known that the shape of the damage crater is approximately circular, independent of micrometeoroid shape or incidence angle. 2 This is because hypervelocity impact is an explosive release of energy in which heat diffuses outward from a point. In this paper, we use a single parameter to describe the impact crater: the DCD. Below, we propose several methods to calculate the DCD.
Another effect of the hypervelocity impact is the contamination of the surface by the ejected material. There is an experimental evidence that ejection takes place due to the strength of the hypervelocity impact. [18][19][20] We cannot write a definite account as to whether contamination goes to the mirror or to the structure around (including other mirrors). The danger with ejected material is that it can contaminate other components in the system, which, for a given contamination level, can generate a higher scattered light contribution to the photodetectors. In principle, the amount of stray light caused by contamination may be of the same order or even larger than that caused by the impact craters. The mechanisms for redeposition go well beyond the scope of the paper, but we can place upper limits by supposing that the impacted mirror cannot not receive more than 100% of the contamination it generates. So, in this paper, we derive an upper limit for the contamination of the impacted mirror, by considering that all ejected matter is deposited back. After DCD estimates, we give expressions to calculate the total mass ejected for single micrometeoroid impact.

Estimate of the damage crater diameter
To calculate the size of the DCD, we use two different models 2,21 with seven sets of parameters in total. The Hörz 21 model is based on the analysis of three laboratory experiments performed by independent investigators to calibrate lunar microcraters. The DCD D (cm) is found to be a function of the mass of the projectile m (g): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 5 3 1 The coefficients C (cm) and Λ are given in Table 1. The DCD as a function of mass for these three Hörz models are given in Fig. 2.
Another model is based on a damage equation, 16 which describes the physics of projectiles impacting a target at high velocities. The DCD D in this case is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 4 5 0 where K 1 is a factor characteristic of the model, d μ (cm) is the micrometeoroid diameter, ρ μ , ρ t (g∕cm 3 ) are the densities of the micrometeoroid particle and target, respectively, v (km∕s) is the impact velocity, α is the impact angle, the crater factor K c is the ratio of the crater radius to the crater depth, and it may be as high as 10 for brittle targets such as Zerodur ðRÞ . 16 The origin of the equation parameters is independent investigations (Gault, Fechtig, McHugh & Richardson, and Cour-Palais), which are summarized in Ref. 16. Typical values of the parameters of Eq. (2) are given in Table 2.
As a conservative approach, the value of the impact angle α is set to 0°. In the LISA telescope, the material for the primary M1 mirror will be Zerodur ðRÞ ceramics with a density of ρ t ¼ 2.53 g∕cm 3 . Zerodur ðRÞ is a brittle material. We assume that the coating on the mirror does not affect the crater formation. We choose K c ¼ 10 as the worst-case scenario.
The results of the different DCD calculations are given in Fig. 4. The difference between different models (up to one order of magnitude) can be explained by the variety of the experimental conditions in the study of hypervelocity impacts, the complexity of the physical phenomena, and the different analytical approaches of the investigations. We assume that some of the used models may overvalue or undervalue DCD for Zerodur ðRÞ material. For this reason,

Estimate of mass ejection
As a micrometeoroid impact is a microexplosion event, some mass will be ejected and can contaminate surfaces including the M1 mirror. Here, we consider this microexplosion process and calculate the total amount of ejected mass M e . To calculate M e , we use the equation derived by Gault 18,19 following the analysis of a range of experimental data: where E i is the projectile kinetic energy in Joule and α is the angle of impact. For brittle materials such as Zerodur ðRÞ , the coefficient K depends on the diameter d μ of the micrometeoroid: To derive a worst-case value, we assume that all the ejected mass will be deposited on the M1 mirror surface. As a result, using Eq. (3), we can calculate the mass ejected for each mass of the micrometeoroid and so to build an appropriate M1 contamination model (see Sec. 3.3.2).

BSDF Calculations
As was mentioned above, the hypervelocity impact of the mirror surface by a micrometeoroid can cause scattering for two reasons: from the impact crater and ejected contamination. Each of them requires a specific analysis.

BSDF due to crater damage
To calculate BSDF(θ) (θ is the scattering angle) from damage craters, we use the Peterson model, which is devoted to calculating the BSDF of the scattering due to digs in optical components. 4 We assume that modeling the scattering from a dig can apply to modeling the scattering from an impact crater of the same diameter. In this model for a crater of a given diameter D, the scattered light is divided into two contributions: • "geometric" scattering or backscattering from surfaces inside the crater, considered as a Lambertian scatterer of diameter D; • diffraction of light that passes around the crater, considered as a circular mask of diameter D.
In the Peterson model, 4 the digs are considered to be circular, and the intensity of diffracted light from a dig (crater) is calculated using the scalar diffraction theory in the Fraunhofer (far-field) limit and Babinet's principle. So, the total BSDF from digs (craters) is calculated as the sum of the geometric and diffraction contributions: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 5 4 2 BSDFðθÞ where N D is the number of digs per unit area, D is the dig (crater) diameter, λ ¼ 1.064 μm is the optical wavelength, and the roll-off angle 4 After integration over the range of micrometeorite masses, the calculated BSDF for different models of DCD is given in Fig. 5. We further proceed with FRED optical software in which the "ABg model" of BSDF is embedded. This model is widely used to describe the scattering due to microroughness 10,22 and involves only three parameters: A, B, and g coresponding to a proportionality factor, a roll-of angle, and a slope, respectively. To implement this Peterson model in the FRED software, we fit the resulting BSDF (see Fig. 5) with the ABg model plus a constant term: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 3 9 7 : The first term in Eq. (5) corresponds to the diffraction of light that passes around the crater, and the second term corresponds to Lambertian scattering of level R inside the crater. The fit of the Peterson BSDF by the ABg model 22 curve was previously used to simulate light scattering by digs in the METIS coronograph. 23 Let us now consider the backscattering from the primary mirror of NASA's model of the LISA telescope. The precalculated BSDF for different models of DCD is shown in Fig. 5. The contribution to the backscattering fraction (BSF) due to the diffractive part varies according to the different models from 1.6% to 8.7% (5.4% on average over the seven models). All the rest is due to Lambertian scattering. So in this particular case, methods based on percentage area coverage 10 give a reasonable estimate of the scattered light due to micrometeoroid damage. The total integrated scatter (TIS) is a ratio of the total scattered power to the incident power. TIS for different models and mission duration is given in Table 3. The values of TIS due to micrometeoroid damage is significant.
The variations of the computed TIS values (Table 3) are large, due to the range of values (up to a factor of 10) for the crater diameter. Presently, there is no reason to prefer or discard one of the models, so we used all available models, from the most pessimistic to the most optimistic. Further work on this topic is required to distinguish the best model to use.

BSDF due to ejected contamination
To calculate the BSDF due to ejected contamination, we assume that the size distribution of particles can be described with the (IEST)CC1246 standard. 24 It describes the number of particles N p (per 0.1 m 2 ) whose diameter is greater than or equal to D p by where S is the slope of particle size distribution, CL is the cleanliness level, and D p is the particle diameter in μm. This is a common way to describe the distribution of the contamination on an optical surface. The model is implemented in FRED optical software and is easy to use as it requires only a few parameters (λ, S, CL, etc.) and relies on the properties of Mie scattering. Here, we will calculate the cleanliness level using the definition of the CL parameter [the largest particle (in microns), which can be found on a surface of 0.1 m 2 ; see Eq. (6)] and using the following considerations: • In the regime of hypervelocity micrometeoroid impacts, the largest ejected particle mass is reported to be proportional to the total mass ejected. 20 For simplification, we assume the worst case when the coefficient of proportionality is equal to one: the biggest ejecta carries most of the ejected mass, up to the total ejected mass (the coefficient is 1). The ratio of 1 is the worst-case scenario. Further numerical and experimental investigation in this field is required. With this assumption, the most massive particle is ejected from the biggest impact micrometeoroid. We apply this assumption only for the largest particle on the surface. The distribution of all the rest ejected particles is assumed to be given by the (IEST) CC1246 standard.
• Ejected mass is mainly target mass (Zerodur ðRÞ ) and has the same density as a target material. Ejected particles are spherical.
• All the ejected mass deposits back onto the surface. Using the flux calculated in Sec. 3.1 and ejected mass in Sec. 3.2.2, we find that after 4 years of exposition, we will have maximum ejected mass equal to 1.14 × 10 −4 g, which corresponds to a diameter of a sphere equal to 441 μm, so CL ¼ 441.
To find the slope S, we use the mass conservation law. The total mass of particles, which is given by this distribution, should be equal to the total mass ejected over the exposition duration: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 6 7 5 where S m is the M1 mirror area, the mass of the micrometeoroid mðxÞ is a function of its equivalent sphere diameter x and M e is the ejected mass given in Sec. 3.2.2, ϕ G ðmÞ is the integrated flux given in Sec. 3.1, Δ is a difference operation between the neighboring bins of the distribution (N p and ϕ G are the cumulative distributions). On the right side of Eq. (7) is the total mass ejected due to all micrometeoroid impacts, and in the left side is the total mass of contaminations, assuming that the distribution of particles on the surface will be given by Eq. (6). When the slope parameter S is chosen correctly, these two masses will be equal. We find that for 4 years, the absolute value of the slope is equal to S ¼ 0.8738, which is quite close to the value 0.926 used in the CC1246D standard. For 10 years of exposition: CL ¼ 743.2 and S ¼ 0.7668.
Up until this point we have ignored the silver coating of the Zerodur ðRÞ mirror. We have assumed that only Zerodur ðRÞ material is ejected. However, we can also assume that the silver coating will play a role in the amount of ejected mass, so in the limiting case we assume that only silver will be ejected. We consider this case only for the mass ejection calculation. In all other cases only uncoated Zerodur ðRÞ material is considered. In this case in Eq. (3), the coefficient is K ¼ 1 for all diameters of the micrometeoroid. The corresponding S and CL coefficients are S ¼ 1.2827 and CL ¼ 216.6 for 4 years of mission duration and S ¼ 1.0998 and CL ¼ 365 for 10 years.

FRED Simulations
The backscattering in the direction of the photodiode of the LISA telescope has been calculated using the FRED simulation software. The scattering model has been applied for the M1 telescope mirror only. The scattering calculation includes two main items: due to impact craters and due to ejected mass contamination.

Scattering due to impact craters
To calculate the scattering due to impact craters in FRED software, we use two embedded BSDF models: Lambertian and ABg. The coefficients of these models were obtained from the fit of total Peterson's BSDF, as was described in Sec. 3.3.1. The calculation results are summarized in Table 4. Despite the high values of the TIS (see Table 3), the values of backscattering are low, as the coupling factor of the M1 mirror is low.
To study the impact of micrometeoroid parameters on the final result, we consider a situation when the density of micrometeoroid follows the step distribution given in Fig. 2 (blue curve), and we simulate the distribution velocities following Taylor's 17 observations, which is a re-evaluation of the Harvard Radio Meteor Project data of about 20,000 meteor observations. The result is summarized in Table 5 (third column). For comparison, the second column contains the values of BSF with the assumptions used in the paper: constant velocity v ¼ 20 km∕s and density ¼ 2.5 g∕cm 3 of the micrometeoroid. No qualitative change is observed. The values are slightly lower, as the considered density of micrometeoroids is lower.
The FRED divot analysis approach was used by NASA 25 to estimate the micrometeoroid damage of the M1 mirror. In their work, the total area occupied by the crater was modeled as a single divot placed on the M1 mirror surface. Depending on the position of the divot on the mirror surface, the BSF obtained in their study is in range from 1.73e-14 to 3.33e-13, which is compatible with the values obtained in this study.

Scattering due to ejected mass contamination
The FRED calculation of contamination was performed with the embedded 1246C 24 standard. The values of CL and S used in the simulation are listed in Sec. 3.3.2. The computed values of backscattering are summarized in Table 6. As Zerodur ðRÞ mirror will be coated with silver, it would certainly cause an effect of contamination population, and so the real value of scattering will be in between two limit cases: mirror material is only Zerodur ðRÞ and mirror material in only silver.
These values are compatible with the highest of the scattering data of Table 4 (6.6e-13 for 4 years and 2.5e-12 for 10 years). So the scattering contribution due to ejected mass is dominant under the assumption that 100% of the ejecta contribute to the M1 contamination and should not be neglected in stray-light estimates caused by the micrometeoroid impacts.

Conclusion
In this paper, we developed a method for the estimation of the stray-light due to the micrometeoroid damage of optical surfaces. It consists of four steps. The first step is the flux calculation   based on the satellite environmental model. The second step is the calculation of the DCD and ejected mass. The third step is the calculation of the corresponding bidirectional reflectance distribution function using the Peterson model and the 1246C standard. The last step is the calculation of the scattering with an optical software, for the detailed application to the considered optical configuration. We have applied the method to the simulation of the scattering of light in the LISA telescope, due to the damage to the primary mirror from the micrometeoroid impact. The results suggest that even under the worst-case assumptions the impact craters and the resulting contamination contribute to an acceptable scattering of light to the LISA detectors. It should be noticed that a contamination due to mass ejection gives a significant contribution. Further work should address a possible contamination due to ejecta toward mirrors other than the primary mirror.
The main sources of uncertainty are in the modeling of the DCD and in the distribution of ejected particles (shape and quantity) on the damaged surface. The modeling presented here should benefit from any future improvements in the experimental data, particularly when optical materials are used as targets for the hypervelocity impact experiments.
The method is straightforward to apply, to modify, and it can be used for any space optical instrument with minor parameter changes. The code is available on GitHub. 26 The model can be used not only for reflective but also for refractive optics as well. The final result of the scattered light calculation depends on the optical design of the telescope, in our case the beam expanding telescope that is required to transmit the emitted beam to the distant spacecraft.