1. Introduction

When designing an optical system, whatever the application, minimizing the stray light is an important aspect of the work. It is unavoidable that some parasitic light will remain in the system and degrade the performance, so it is essential to identify the origin of this unwanted light in order to reduce its impact. Stray light is often referred to as flare or veiling glare in an imaging system and defined in the ISO 935 and 18844 [1,2] as the light falling on an image which does not emanate from the subject point. The term “unwanted light” summarizes best what most people define as flare [3]. Depending on the application, it can create different issues.

For visible imaging applications flare is a well-known observed phenomenon and can sometimes be used deliberately for artistic reasons in photography, movies, or video games. It can arise from different origins in the system: bright sources outside the field of view of the system that find a path to the sensor, internal reflections on the optics (lenses, filters) or mechanics of the system, scattering due to material imperfections (roughness for example), diffraction on the elements of the system, etc. Flare reduces the image quality as it can reduce the Signal to Noise Ratio (SNR) of the sensor, reduce the dynamic range of High Dynamic Range (HDR) cameras, create saturated areas of the sensor, or create overlapped features in the final image [3–9]. For these reasons, it is seen most of the time as an artifact that should be reduced for better image quality. The two pictures in Fig. 1 illustrate some of the different types of flare on images; what can be seen in the left picture around the bright source is called petal flare and will be explained and reproduced in this work in the final section.

 

Fig. 1. Examples of flare on images (copyright-free images from Ref. [10,11]).

Download Full Size | PDF

For time of flight (TOF) applications in the near infrared range where the distance to the object is evaluated to reconstruct a 3D image of the scene (either with global shutter pixel technology or single photon avalanche diodes), stray light can lead to errors in depth data accuracy [12–15]. The degradation can be quite important, for example when a weak signal from a distant object is affected by strong scattering/reflectivity signals from close objects as explained in Fig. 2. As an example, Mure-Dubois [12] reports an error in the range of 100mm to 400mm at 4m distance.

 

Fig. 2. Light scattering in a TOF camera (Ref. [12,13]).

Download Full Size | PDF

Even though a lot of work can be done for these different applications to reduce these stray light issues by post-processing of the data or calibration of the system [7,12–14,16–20], it is critical to identify the origin of the problem in the optical module to reduce the amount of stray light perturbing the system.

Analytical calculations in the past [21–26], and ray-tracing software today [19,27–29] are powerful methods to model light propagation in optical systems and estimate the degradation of performance due to stray light. These studies require a knowledge of the surface properties of the elements of the system. The principle of the analytic calculations done in the past to model stray light in optical modules is to consider the diffraction and scattering of the light on the optical elements of the system. For that purpose, the Bidirectional Scattering Distribution Function (BSDF) is widely used to describe the scattering due to the surface roughness of optical elements [30]. Peterson [22] has for example used the Harvey model [31] to describe and calculate the scattering of rough surfaces of optical elements in an imaging system. More recently Caron, Mariani, and Ritt have extended these calculations for their work and the readers could refer to these interesting reviews and works for more details [25,26,28,29,32]. Of course, any BSDF model could be defined on the surfaces and nowadays most ray-tracing software tools available in the market have different scatter models that can be applied directly on the surfaces of the optical elements and the housing of the system (ABg model, Lambertian, Mie scattering, imported BSDF data, etc.). When only the reflectance of a surface is required for the analysis (as studied here in our work), the bidirectional reflectance distribution function (BRDF) rather than the full BSDF, which includes both reflectance and transmittance, can be used.

Nevertheless, in all these studies, the reflectivity of the detector itself in the parasitic light of the system is either not considered or a simple coating model on the surface of the detector is used to simulate a given amount of specular reflectivity. But considering the wavelength of the applications (visible for photography, or near infrared range for TOF and LIDAR applications as some examples) and the size of the pixels of the detector in the market nowadays (between 1µm to 4µm for standard CMOS imagers and between 5µm to 10µm for TOF pixels or SPAD) the detector will act like a diffraction grating and reflect power in different directions corresponding to the grating orders. Analytical models of a grating could be applied to the detector surface to obtain the correct diffracted order directions depending on the wavelength, the pixel size, and the incident angle, but they cannot predict what fraction of the incident power will be reflected to each diffracted order. To be able to estimate the correct diffuse reflectivity of the detector for each incidence angle of the source (i.e. the BRDFs of the pixels) an electromagnetic simulation of pixels solving Maxwell equations is required. The finite-difference time-domain (FDTD) method is now widely used in the optical community for the simulation of imaging or TOF sensors [33–37], as well as for the simulation of surface reflectivity [38–41]. Thus, this solver allows us to obtain accurate BRDFs of our pixels, and subsequently use them as the surface property of our detector in the ray-tracing simulation for stray light analysis of the system. The coupling between two simulation worlds (electromagnetic and ray-tracing) with vastly different scales is an emergent field in the scientific community for different applications [42–50]. In this work we present a demonstration of such an FDTD to ray-tracing simulation workflow that couples the light propagation at different spatial scales (µm and mm) in order to analyse stray light in an optical module. We believe this methodology could be used to complement and improve the stray light model in the aforementioned references in this field.

The work is organized as follows: in Section 2, we describe the principle of the simulation flow, explaining the constraints considered and the general methodology; in Section 3, we detail the BRDF calculation by electromagnetic simulation and demonstrate this methodology by comparing with the experimental BRDF characterization of some pixels; finally, in Section 4, we apply this coupling methodology to reproduce some stray light effects observed in the characterization of real modules.

2. Micro to macro scale coupling simulation methodology

The methodology we present here allows us to simulate the BRDF of our sensor pixels by electromagnetic simulation and use it as a surface property in our ray-tracing simulation for image sensor or TOF applications. But it is important to note that this methodology could be extended to every application where a BRDF needs to be simulated with electromagnetic simulation and used afterwards for ray-tracing simulations. In general, the electromagnetic simulation evaluates the coherent scattering properties of the surface, from which the incoherent scattering properties required for the BRDF format can be generated. In this work we use Ansys Lumerical FDTD software [51] for the electromagnetic simulation, and Zemax OpticStudio [52] for the ray-tracing simulation, but this simulation flow is compatible with any other software.

A schematic of an optical module for visible imaging is represented in Fig. 3. In this case it is composed of a housing, several lenses, an infrared filter, and a CMOS sensor. The light reaching each pixel is defined by two main parameters: the angular distribution which corresponds to the f-number (f/#) of the objective (ratio between its focal length and its entrance pupil diameter) and the chief ray angle (CRA) which is the ray passing through the center of the aperture stop and the center of the pixel of interest. Figure 3(a) shows what happens for pixels at the sensor center with the insert showing a zoomed view of two Backside Imager (BSI) technology pixels (blue and green pixels) with schematic geometrical rays and electromagnetic simulation of a green wavelength. BSI technology, in comparison with Frontside Imager (FSI) technology, has better optical performance because metallic interconnections are not in the light path [37], and is now the standard in the CMOS image sensor industry. Figure 3(b) shows what happens for pixels at the sensor edge with the insert also showing a zoomed view of two BSI technology pixels (blue and green pixels) with schematic geometrical rays and electromagnetic simulation of a green wavelength. The pixels form a locally periodic pattern on the image sensor surface, and the pitch refers to the minimum distance between pixels. For color sensitive detectors, the pixels are arranged into a 2 by 2 locally periodic structure with red, green and blue (RGB) color filters, called the Bayer pattern, with a period that is twice the pitch, as shown in the top view of Fig. 3.

 

Fig. 3. Representation of an image sensor module. Figure 3(a) shows rays for an on-axis field and Fig. 3(b) shows rays for an off-axis field with inserts representing light propagation inside two pixels (blue and green) of a 1.1µm pitch sensor (on the left is a schematic of geometrical rays, and on the right the electric field intensity from an FDTD simulation of a green wavelength). Between the inserts is a top view of the RGB sensor with the periodic Bayer RGB pattern (2 by 2 pixels) simulated by FDTD.

Download Full Size | PDF

As represented in Fig. 3(b), to re-center the off-axis incoming light on the photodiode, color filters (if present) and microlenses are shifted in most commercial sensors today. Each pixel on the matrix has a different design because the radial shift of each filter/microlens is optimized depending on the CRA corresponding to the position on the sensor and, as a result, the optical efficiency and the reflectivity will also be different. Thus, the pixel BRDFs must be calculated for several incidence angles (defined by the illumination cone and the CRA), for different pixel designs (defined by the position on the matrix), and for the wavelengths of interest. Of course, on the detector, the local variation of the pixel design and illumination is small and simulating only a few different areas of interest is sufficient to estimate the local BRDFs over the entire sensor plane. We therefore simulate different small groups of pixels at different locations on the sensor (i.e., corresponding to different fields in the scene) that receive the same uniform illumination; the spatial extent of each group is small compared to the spatial variation of the source. The reader could refer to our previous work that describes the CMOS simulation methodology for more details [33]. At the pixel level, light is distributed uniformly: spatially over the pixel area and angularly inside a cone defined by the exit pupil of the objective and the pixel. With electromagnetic simulations, it is possible to simulate these groups of pixels for different incidence angles (with a plane wave source) to get the corresponding BRDF as will be explained in more detail in Section 3. Finally, in the ray-tracing software, these BRDFs are imported and used as surface properties at the corresponding location on the detector plane. This electromagnetism to ray-tracing simulation flow is summarized in Fig. 4.

 

Fig. 4. Representation of the electromagnetism to ray-tracing simulation flow for stray light analysis of an optical module considering the reflectivity of the detector. As noted on the right, each incident angle (and indeed polarization) must be simulated separately because FDTD is a fully coherent electromagnetic solver.

Download Full Size | PDF

3. BRDF with electromagnetic simulation

The BRDF is the ratio of reflected radiance to the irradiance of a surface for a given incident angle, reflected angle, and wavelength. To fully characterize the BRDF using FDTD (or similar electromagnetic solver) requires many individual simulations, one for each angle of incidence and polarization (S or P). Multiple plane waves at different incident angles cannot be solved in the same simulation because FDTD is fully coherent and interference between the different sources would result. However, the frequency-dependent response can be extracted from a single FDTD simulation, with some care, because it is a time domain method so many wavelengths can be solved in a single simulation. To calculate the BRDF we must consider the reflected light in the far field, where it behaves as a plane wave with well-defined angle and polarization. This is achieved by performing a near to far field projection of the FDTD near field results. In the far field, quantities such as the Poynting vector can be easily determined from the electric field intensity which allows us to calculate radiance quantities required for the determination of the BRDF. The irradiance can be determined from the plane wave source conditions used in the FDTD simulation. Unpolarized results can be easily calculated from the full polarization information available from the FDTD simulations. In the remainder of this section, we explain the full details of this BRDF calculation with Ansys Lumerical FDTD.

The pixel structure is assumed to be locally periodic, which is a good assumption given that the shifts in the microlenses and color filters occur slowly from the center to the edge of the image sensor, on length scales which are much larger than the wavelength. The simulation can therefore be reduced to the smallest, periodically replicated group of pixels, which is often four when a Bayer pattern is used for visible applications (as shown in the top view of the sensor in Fig. 3). In this locally periodic approximation, the Bayer pattern forms a two-dimensional grating with well-defined primitive lattice vectors ${\vec{a}_1}$ and ${\vec{a}_2}$. The excitation source is a plane wave with a given incident angle and polarization. For boundary conditions, we use perfectly matched layer (PML) absorbing boundaries above and below the pixels to absorb both the reflected light and the light transmitted to silicon substrate carrier. On the lateral boundaries, either Bloch-periodic boundaries or split field boundaries, also known as the broadband fixed angle source technique (BFAST), can be used. The Bloch-periodic boundaries result in a wavelength dependent incident angle, while the BFAST method results in a substantial increase in simulation time, especially at steeper incident angles. In this work, we chose Bloch-periodic boundaries because the wavelength dependence of the source angle can be accounted for during the analysis. During the FDTD simulation, the transient electromagnetic fields are Fourier transformed to the frequency domain and normalized to the source pulse used to excite the system [53]. By calculating these fields on a plane above the image sensor surface, we can calculate the Poynting vector in the frequency domain and therefore the total reflected power, normalized to the source power, at any desired wavelength. Typical sizes of the FDTD grid are on the order of λ_min/10 where λ_min is the minimum wavelength (in each material) over the desired bandwidth and standard methods can be used to ensure convergence.

In order to calculate the reflected power to each angle (or each diffracted order), we must perform a near to far field projection using well known Green’s function methods [53]. In general, calculating the far field from the near field involves evaluating integrals of the tangential electric and magnetic fields over a near field surface for each position in the far field. When the far field position approaches an infinite distance compared to the dimensions of the near field surface (the Fraunhofer diffraction limit), the far field position can be defined by a wave vector $\vec{k} = {k_0}\vec{u}$ and a distance $R = |{{{\vec{r}}_{far}} - {{\vec{r}}_{near}}} |$ where k₀ is the wavenumber in the far field medium, $\vec{u}$ is the unit vector given by $\vec{u} = ({{{\vec{r}}_{far}} - {{\vec{r}}_{near}}} )/R\; $ and ${\vec{r}_{near}}$ is the center of the near field surface. The integrals over the near field surface can then be reduced to a series of Fourier transforms of the near electromagnetic fields.

Here, we know the reflected near fields on a surface corresponding to one period of a periodic system. We can therefore calculate the contribution of a single unit cell to the far field, $\vec{E}_0^{far}({\vec{k}} )$ from the near field of a single period, $\vec{E}_0^{near}({\vec{r}} )$. From there, we can calculate the reflected far field of a finite number of periods simply by summing the contributions of each period with the appropriate phase

From Eq. (1), we can calculate the normalized power (and polarization, if desired) of reflected light at any angle from the surface, as a function of the incident angle and polarization. This can be done either in the limit of an infinite or a finite number of periods. This information on the reflective properties of the image sensor surface can be saved for subsequent analysis in ray-tracing. We note that the FDTD simulation provides the fully coherent results, including phase and polarization, which can be post-processed into incoherent results as required when passing the information to a ray tracer. There are several methods which can be used to store the reflective properties for ray-tracing. For example, the grating order efficiencies could be used directly to determine the fraction of rays reflected to each order, while the direction could be determined from Eqs. (2) and (3). While this method is preferable in principle, support for wavelength and incident angle dependent diffraction of 2D gratings is not always available in ray-tracing tools. An alternative is to use the tabular BRDF data format that most ray-tracing software can handle easily. When using a tabular BRDF data format, the diffracted order efficiencies cannot be employed directly because the reflectance of a grating as a function of output angle is composed of Dirac delta functions that cannot be easily interpolated. Instead, we want to treat the ideal plane waves as finite sized beams and thereby create a BRDF that is smooth enough to allow for interpolation. One method to achieve this is to apply a weighting function in Eq. (1) of:

This is approximately equivalent to the angular reflection that will occur for a finite sized incident beam of radius R₀. Here, we are neither so concerned with the precise value of R₀ nor are we concerned with the specific choice of function used, but rather with smoothing the Dirac delta functions present in a BRDF of an ideal grating into functions that can be managed through the tabular BRDF formalism available in ray-tracing software. Choosing a larger value of R₀ gives narrower diffracted beams, which requires sampling the BRDF at more angles. A smaller value of R₀ is less accurate but allows sampling the BRDF at fewer angles. More quantitatively, the above weighting function will result in a far field composed of multiple beams of the form

In order to demonstrate this electromagnetic BRDF simulation methodology of pixels, we performed some BRDF characterizations by measurement. The reader could refer to a Guarnera publication [54] for a complete review of BRDF representation and acquisition. The characterizations were performed with a mini-Diff equipment from LightTec [55] at 940nm at normal incidence on two different process technologies with different pixel pitch (4.6µm pixel and 2.2µm pixel without color filters). The illumination spot size in characterization was around 1mm diameter. We performed the measurements on samples without a radial shift of the microlenses (i.e., all pixels were identical) in order to be equivalent to the simulation conditions. Figures 5(a) and 5(b) show, respectively, the results from characterization and simulation on the 2.2µm pixel pitch. We see that both results qualitatively agree in terms of energy and direction angles corresponding to the ones predicted by the grating equation (Eqs. (2) and (3)) for the 940nm wavelength at normal incidence in the air. A log scale is used to emphasize the signal in the diffracted orders. We also notice that the weak 2^nd orders are not detected in characterization, due to the limited detector responsivity of the equipment for these conditions [55]. Figures 5(c) and 5(d) show, respectively, the results from characterization and simulation on the 4.6µm pixel pitch. Once again, the simulation results qualitatively match the characterized ones with angles corresponding to the grating equation. The small differences observed could be explained by the simulation model. Even though we try to model the pixels as accurately as possible, there are always some differences with the final process realized on a silicon wafer (refractive indices of the materials, layer thicknesses, microlens shape, isolation trenches profile etc…). Furthermore, the simulated beams have an angular width that is determined by our choice of R₀, while the experimental results are limited by the angular resolution of the measurement device.

 

Fig. 5. BRDF characterization and simulation comparison for different pixel pitch at 940nm at normal incidence. Figures 5(a) and 5(c) show characterization BRDF result for 2.2µm and 4.6µm pixel respectively. Figures 5(b) and 5(d) show FDTD simulation BRDF result for 2.2µm and 4.6µm pixel respectively. The white rings correspond to the polar angles in 10 degree increments up to 90 degrees.

Download Full Size | PDF

We also indicate in Table 1 the absolute values of reflectance from simulation and experimental characterization at 940nm. For reflectance measurements on scattering samples, it is often interesting to combine two different tools: a BRDF equipment like the mini-diff equipment we use here (or a more sophisticated goniometer if available) that gives the angular information but captures less accurately the total diffuse reflectance, and an integrating sphere equipment that accurately measures the total reflectance value but without angular information. For the latter we performed the measurement with a Lambda 1050 from Perkin Elmer [56] to compare the total reflectance with simulation. The values indicated in the Table 1 for the measurement correspond to the mean value measured over 5 samples with uncertainty reported in parenthesis. The results in Table 1 show a good agreement between simulation and characterization (considering the uncertainty of the measurements), for both specular and total integrated reflectance. One should notice from the Table 1 that the reflectance can be quite significant for some pixel technologies, especially at 940nm on BSI technology due to the low absorption of silicon and therefore the light reflected from the metallic interconnections under the silicon substrate. As already discussed, it can create important stray light issues in the optical module that have to be identified and reduced.

Table 1. Image sensor reflectance results comparison

View Table

As it is difficult to accurately characterize by measurement local BRDFs on sensors with a radial shift of microlenses (it would require spot sizes of a few microns in diameter with good angular detection accuracy and wavelength scanning of the source), the FDTD simulation methodology is an efficient way to get the BRDFs of pixels at different incidence angles and wavelengths in order to predict potential issues in the final optical module.

4. Application to stray light analysis

In this section, we use our micro to macro scale optical simulation coupling flow to reproduce and explain stray light issues observed on module. The module is a standard one used for imaging applications in the visible, composed of several lenses and an infrared filter to cut off the signal above 650nm (see the module representation in Fig. 4). We use it in combination with two different sensors: one with a 1.4µm pixel pitch, and another with a 1.75µm pixel pitch. Note that due to the Bayer pattern, the period of the grating simulated by FDTD is twice the pitch. We have observed with laboratory experiments important red ghost images on the sensors when illuminating the module with a bright light source. The ghost patterns observed were composed of several secondary images whose positions in the image depended on the sensor used (i.e. different pixel pitch), which were clear indications that the problem may have arisen from sensor diffraction reflectivity. To confirm this, we performed electromagnetic BRDF simulations of these two pixel pitches at the sensor center and imported the data into the ray-tracing software in order to find the element(s) on the module causing the ghosts. In Fig. 6 we show the BRDF for these two pixel pitches at 630nm wavelength for normal incidence, as well as the incoherent sum of the BRDFs over a sampling of incident angles over the full aperture of the module (F/2.8) at the sensor center.

 

Fig. 6. BRDF for the two pixel dimensions simulated at 630nm. Figures 6(a) and 6(c) show the BRDF at normal incidence for the 1.75µm and 1.4µm pixel pitch respectively. Figures 6(b) and 6(d) show the incoherent sum of the BRDFs corresponding to a sampling of incident angles over the full aperture (F/2.8).

Download Full Size | PDF

The incoherent sums of the BRDFs in Fig. 6(b) and 6(d), showing the sampling of multiple incident angles over the aperture, illustrate the angular spread and power of each diffracted order due to the f/# of the module at this wavelength. The use of an incoherent sum has been justified in previous work [33], and also corresponds to what will be observed when the BRDF is used in a ray tracer where each ray is treated incoherently. We can observe that the sampling of incident angles in the FDTD simulations should be much higher to obtain results corresponding to a full diffuse source. Indeed, if enough incident angles were simulated, we would expect to see smooth circles with a diameter equal to the full aperture centered at each diffracted order in Figs. 6(b) and 6(d). Therefore, Figs. 6(b) and 6(d) underscore some of the challenges of working with a tabular BRDF data format to represent diffraction from gratings because, for ideal accuracy, a large number of incident angles, output angles and wavelengths would be required to correctly represent the scattering behavior. Nonetheless, the relatively small number of incident angles used are sufficient to obtain qualitative agreement with experiment, as we will see in Fig. 8, despite the fact that artifacts of the limited angular sampling are visible in the simulated results. Improvements to avoid the use of a tabular BRDF data format, such as the use of plugins to better handle the diffraction order efficiencies of the periodic image sensor surface, while analytically accounting for the directions of the diffracted orders, will be the subject of future work.

Once the FDTD electromagnetic simulations were performed, the BRDF data was imported into the ray-tracing software using a tabular text file format and applied as a surface property at the center sensor plane as explained in Section 2. The results of the final ray-tracing simulation at 630nm on the two module configurations are shown in Fig. 7 as well as the captured images with a bright light source. The measured images show red ghost patterns that are well reproduced with the simulation flow. The size, spacing, and directions of the ghosts are comparable between measurement and simulation in the two module configurations (yellow diagonal lines and yellow dot squares have been added to the images to more easily compare dimensions between the patterns).

 

Fig. 7. Comparison images with ghost patterns between measurement and simulation. Figure 7(a) shows the image captured with the module using the 1.4µm pixel sensor. Figure 7(b) is the corresponding ghost ray simulation at 630nm. Figure 7(c) shows the image captured with module using the 1.75µm pixel. Figure 7(d) is the corresponding ghost ray simulation at 630nm.

Download Full Size | PDF

To better understand the origin of the ghosts, a detailed stray light analysis in ray-tracing software is necessary. The analysis shows that the ghost pattern comes from light that is reflected to various diffracted orders of the periodic sensor that is again reflected on the IR filter and comes back as ghosts captured by the sensor. In Fig. 8 we have represented some simulated rays to illustrate the optical path leading to ghost images on the sensor. The ghost pattern is thus highly dependent on the wavelength, the pixel pitch, the incidence angle on the sensor, the infrared filter to sensor distance and the thickness of the infrared filter. This demonstrates that electromagnetic simulation of the sensor coupled with ray-tracing simulation of the module is the only way to reproduce this kind of behavior.

 

Fig. 8. Ray-tracing simulation at 630nm of the module with BRDF electromagnetic simulation data applied at the sensor surface. The zoomed view on the right shows the rays generating the ghosts on the sensor, due to small reflections on the sensor (that behaves as a weak reflection grating), that are again reflected by the IR filter back towards the sensor.

Download Full Size | PDF

For this specific issue, a variety of solutions could be explored, such as optimizing the sensor design to reduce reflectivity, changing the IR filter to sensor distance to defocus the ghosts, or optimizing antireflective coatings on the IR filter.

5. Conclusion and discussion

In this paper, we present a new methodology for analyzing stray light in optical systems that originates from image sensor reflections. This methodology is based on a reflection analysis by electromagnetic simulation of the sensor’s locally periodic pixels, coupled to a ray-tracing simulation of the module, and could be used for many different applications including imaging and time of flight modules. This methodology could also be applied at any wavelength and any surface for which a BRDF can be calculated with electromagnetic simulation. We have shown that we are able to reproduce by simulation ghost patterns seen on real images. Thus, this methodology could be a useful additional tool for engineers to detect and solve stray light issues during the design of modules.

The method presented is not limited to a tabular BRDF file format, which was chosen here because it is supported by most ray-tracing software. Indeed, a tabular BRDF file format is limiting for representing reflections from periodic gratings because it requires an artificial broadening of the diffracted orders calculated by electromagnetic simulation, as well as high angular sampling. Furthermore, if more than one wavelength is included, it requires high spectral sampling due to the shift in grating order angles with wavelength. Therefore, a more direct approach in ray-tracing to handling the grating order efficiency information extracted from electromagnetic simulation is desirable and will be the subject of future work.