Numerical Models for Exact Description of in-situ Digital In-Line Holography Experiments with Irregularly-Shaped Arbitrarily-Located Particles

We present the development of a numerical simulator for digital in-line holography applications. In-line holograms of arbitrarily shaped and arbitrarily located objects are calculated using generalized Huygens-Fresnel integrals. The objects are 2D opaque or phase objects. The optical set-up is described by its optical transfer matrix. A wide variety of optical systems, involving windows, spherical or cylindrical lenses, can thus be taken into account. It makes the simulator applicable for design and description of in situ experiments. We discuss future applications of this simulator for detection of nanoparticles in droplets, or calibration of airborne instruments that detect and measure ice crystals in the atmosphere.


Introduction
With the development of electronic array sensors as CCD (Charge Coupled Device), or CMOS (Complementary Metal Oxide Semiconductor), with sufficiently small pixels to fulfill the Shannon conditions, Digital Holography has expanded rapidly [1].This technique now has applications in a large number of domains from investigation of particles in flows [2][3][4][5][6][7], to visualization of cells in biology or medicine, to phase contrast metrology, or to detection of nanoparticles without being exhaustive [8][9][10][11][12][13][14].Research activities in this field have many orientations: some articles treat the elaboration of models to compute and predict holograms.Others emphasize on the development of new image processing techniques for object reconstruction, estimation of ultimate sensitivity of the set-up, or design of specific configurations for in situ applications [15][16][17][18][19][20].In digital holography, the in-line configuration allows the design of a simple set-up with a relatively low number of optical elements.Depending on the application (visualization of bubbles in a cavitation tunnel, visualization of droplets through a transparent engine), the set-up has to be carefully designed.In order to optimize the design, we have developed analytical models that calculate theoretically in-line holograms generated by various kinds of spherical opaque or transparent particles [19,[21][22][23][24].Our analytical models allow to determine precisely reconstruction parameters from parameters of the set-up, through estimation of the coefficients of its optical transfer matrix.Nevertheless, these previous models are limited to the description of relatively simple circular particles, located in the center of the optical axis of the incident laser beam.There is no simulator that can predict the in-line holograms generated by objects that are arbitrarily-shaped and arbitrarily positioned in the field of the incident laser beam.In this paper, a numerical model that can calculate in-line holograms generated by arbitrarily-shaped arbitrarily-located objects is described.It is based on the scalar diffraction theory.Objects can be either opaque, transparent or semi-transparent.They are described by a 2-Dimensional transmission function.The whole set-up is described by optical transfer matrices of all elements.This formalism allows the description of in situ set-ups, that combine beam shaping, protection windows, different kinds of spherical or cylindrical lenses and objectives.Digital reconstruction of opaque and phase objects using the 2D-fractional Fourier transform is then shown to illustrate the potentiality of the simulator.
Section 2 will present the principle of the simulator.Attention will be focused on the description of the optical system by its own optical transfer matrix, and the definition of the objects (irregularly-shaped, arbitrarily-located, opaque or phase objects).The method to evaluate the intensity pattern in the plane of the electronic array sensor will be described in section 3: it is based on the use of generalized Huygens-Fresnel integrals.We will finally introduce in section 4 two main applications of our simulator: detection of inclusions in droplets for detection of nanoparticles, and analysis of ice-particles for calibration of airborne instruments that can detect and measure the size of ice crystals in the atmosphere.

Amplitude Distribution of the Beam in the Plane of the Object
The basic idea of Digital In-line Holography (DIH) is to record the intensity pattern diffracted by an object on an electronic array sensor.The object is then reconstructed numerically using image processing.Figure 1 is the typical set-up of a DIH experiment using an incident Gaussian beam.This specific case is well adapted to our experiments that use in most cases a single-mode laser beam.Part 1 allows an appropriate illumination of the object.Part 2 is the "imaging" system (leading to a "transport" of the diffraction pattern from the object to the sensor).In order to interpret precisely the experimental results, it is necessary to develop a general numerical model that can describe any set-up.We present here the simulator that has been developed, using the scalar theory of diffraction.At the beginning of the propagation, the waist of the incident Gaussian beam is located at the origin (z = 0).z is the longitudinal axis of propagation (horizontal axis on Figure 1).In the next equations, the transverse axes will be noted µ and ν in the plane of the laser source, ξ and η in the plane of the particle, x and y in the plane of the CCD sensor.The electric field amplitude in the incident plane can be written as follows: where ω 0 is the waist of the beam.Then the laser beam propagates through the first optical system from the incident plane to the object plane.This part of the optical set-up is described by an optical transfer matrix M 1 as can be seen in Figure 1.The complex field amplitude of the wave in the plane of the object is noted G 1 .It is obtained from a generalized Fresnel integral [19,[25][26][27].G 1 is given by: where 1 is the optical path between the incident plane and the object plane.A x,y 1 , B x,y 1 , and D x,y 1 are the matrix elements of the matrix M x,y 1 , where M x 1 = M y 1 = M 1 for an axisymmetric system.M x 1 and M y 1 will be different for a cylindrical system.ξ and η are the transverse coordinates in the object plane.Substituting Equation (1) into Equation (2), the complex amplitude G 1 (ξ, η) in the object plane becomes: with , ω 1x,y = λB x,y , where ω 1x,y and R 1x,y are respectively the beam waist radii and the beam curvatures in the object plane.They are deduced from the classical formula of Gaussian beams.It is important to note that this formulation allows the description of cylindrical systems, by separation of the two transfer matrices along both transverse axes : M x 1 and M y 1 .The details of this calculus can be found in the appendix of reference [19].

Definition of the Objects
In order to describe propagation of light through the second part of the system, it is then necessary to evaluate the transmission coefficient of the diffracting element.We have developed a library of opaque, semi-transparent or transparent phase objects.They are described by a transmission coefficient T (ξ, η).The objects are created from the collection of small circular disks.Both opaque and phase objects can be described.Opaque objects are defined by a collection of disks whose transmission coefficient is 0 inside the disk and 1 outside.For transparent objects, the transmission value inside the object is equal to e iφ and equal to 1 outside the object.The procedure is slightly longer.We realize a collection of circular opaque objects.The transmission coefficient of the object is then modified.We can further create an object where φ is not constant.For both opaque and phase objects, we can adjust the characteristic dimensions of the object.Some examples will be presented in the next sections.

Amplitude and Intensity Distributions in the Plane of the CCD Sensor
A second generalized Huygens-Fresnel transform gives the expression of the electric field in the CCD plane: G 2 .It is given by the following relation: where 2 is the optical path between the object plane and the CCD plane.A x,y 2 , B x,y 2 , and D x,y 2 are the coefficients of the matrices M x,y 2 which describe the second part of the optical system (see Figure 1).We have M x 2 = M y 2 = M 2 for an axisymmetric system.T is the transmittance of the object.x and y are the transverse coordinates in the sensor plane.This second integral is evaluated numerically.The intensity distribution recorded by the CCD sensor is then obtained from: where the overline refers to the complex conjugate.

Hologram Analysis by Fractional Fourier Transform
The next step is the reconstruction of the objects from the holograms.Different tools can be used: wavelet transform, Fresnel transform.We use here the 2 Dimensional-FRactional-order Fourier Transform (2D-FRFT) [28][29][30][31].The definition of the 2D-FRFT that we use is: with α p = apπ 2 (p = x, y) and 0 ≤ |α p | ≤ π/2.The kernel of the fractional operator is defined by: with C(α p ) defined by: Coefficient C(α p ) enables conservation of energy.Parameters s p (p = x or y) are defined from the experimental set-up according to [32]: where N p are the numbers of pixels of the image along x or y axes, and the constants δ p are the pixel sizes along both axes (p = x or y).In practice δ x = δ y .To reconstruct the image of the particle, we evaluate the 2D-FRFT of the diffraction pattern.The diffraction pattern contains a term with a linear chirp ϕ.This linear chirp is easily recovered through 2D-FRFT operation.The kernel of the 2D-FRFT operator is indeed composed of a linearly-chirped function, leading to applications in spatial optics when considering diffraction patterns (as here) or in temporal optics when considering chirped pulses [33].The best reconstruction plane is reached when [19]: where ϕ a is the phase contained in 2D-FRFT and ϕ is the phase contained in the diffraction pattern.From this condition, the optimal fractional orders of reconstruction α opt x and α opt y can be obtained.The existence of analytical relations that link optimal order of reconstruction to experimental parameters of a set-up are particularly interesting for different reasons: first, it allows 3D-localization of the particles under study, whatever the set-up is.In addition, it allows localization of defaults, present in the optical system, that create their own diffraction pattern on the sensor (which acts as noise).For example, we could localize the position of a scratch on a plexiglas cylinder.It was deduced from the analysis of the experimental hologram with these analytical relations.It was then confirmed experimentally watching the inner of the pipe once the set-up was uninstalled [23].As a perspective, we think that our simulator combined to such analytical relations will allow to understand the sources of speckle noise that can alter the quality of holograms.

Simulation of Two Objects in Different Longitudinal Planes: Mixing of Opaque and Phase Objects
In this section, two irregular objects are generated in two different longitudinal planes.Each object has its own size and shape as shown in Figure 2. The object in plane 1 is a phase object (left of Figure 2) and the object in plane 2 is an opaque object (right of Figure 2).The numerical setup is presented in Figure 3.We evaluate separately the two diffraction patterns created by object 1 and by object 2 respectively.As both grids do not necessarily have the same pixel size (due to sampling conditions used when realizing both numerical Fresnel integrals), we use an interpolation function to super-impose both diffraction patterns.Noise can be added as well.As in [20], we assume that multiple scattering can be neglected (i.e., the part of the light diffracted by the first object and then by the second object).For the two objects, the system introduces magnification factors between the object plane and the CCD plane.They are given analytically by the following relation: where g 1,2 x,y are the magnification factors between Plane 1 (or Plane 2) and the CCD plane.w 1,2x,y are the beam waist dimensions (along x and y, respectively) of the unperturbed Gaussian beam in Plane 1 (or in Plane 2 respectively).w totx,y are the beam waist dimensions of the unperturbed Gaussian beam in the plane of the CCD sensor.They are given by: A p tot and B p tot (with p = x or y) are the matrix elements of the transfer matrix describing the whole system between the laser source and the CCD sensor: M p tot = M p 2 × M p 1 .In the simulation, the sensor size is 1024 × 1024 pixels and each pixel is 4.4 × 4.4 µm 2 .The magnification factor introduced by the system for the object in Plane 1 is equal to 6.8, while it is equal to 5.3 for the object in plane 2. Figure 4a is the hologram of the two objects recorded by the CCD sensor.Figure 4b,c show the optimal reconstruction of the two objects.For reconstruction of the objects, the 2D-FRFT is applied to the hologram.The optimal fractional orders for the objects in Figure 4b,c are -0.600π/2 and -0.686 π/2, respectively.These two optimal fractional orders of the reconstructed objects give the correct longitudinal positions along z-axis of both objects.These reconstructions show the difference between opaque and phase objects: for the pure phase object in Figure 4b, a bright and a dark contour due to the discontinuity are observed at the edge of the object as mentioned in [13].The intensity inside the reconstructed object can be higher than the continuous background.For the opaque object in Figure 4c, the intensity in the reconstructed object is zero.Moreover, the analysis of the patterns of the reconstructed objects leads to the size and the 2D-shape of the objects.As classically observed with in-line holography, the reconstructed objects are surrounded by fringes due to the twin image.They "scramble" the reconstruction.Nevertheless, the reconstruction remains possible and a quantitative analysis can be carried out from the reconstructed objects.Intensity profiles of these phase and opaque reconstructed objects are plotted in Figure 5a,b, respectively.The size of the phase object can be estimated from the distance between both discontinuities which are indicated by two arrows in Figure 5a.For the opaque object, the width between both discontinuities of the intensity profile represents its size.From the intensity profile of the two objects and their magnification factors (5.3 and 6.8, for the phase and opaque objects, respectively), the real size of the phase object along x-axis is 24.8 µm and the real size of the opaque object along x-axis is 20.1 µm, respectively.Size measurements correspond to the original sizes of the two simulated objects in Figure 2.For different values of the phase shift introduced by the irregularly-shaped phase object, the transverse reconstructed intensity profiles for ϕ = 0.25 π, 0.5 π, and π are presented in Figure 6 (in solid line, dotted line and dashed line respectively).The curves have been plotted on the same figure.A continuous background has been added on these plots to separate them for clarity (it should be suppressed for a quantitative analysis).A reconstructed "phase jump" is clearly observed on the transverse profiles of Figure 6: it is the abrupt transition between points A and B (points plotted in the case ϕ = π).The amplitude of the discontinuity is directly linked to the phase shift introduced by the object.It has been shown that parameter y A −y B yc is a monotonic function versus ϕ for a circular object [13].y A , y B and y c are normalized intensities measured on the three points A, B and C of the transverse profile defined in Figure 6.The present simulator allows to evaluate this parameter for any complex phase object.The method is efficient when 0.2 π < ϕ < 1.6 π.In other cases, the discontinuity is too small to be observed experimentally.Parameter ϕ can be recovered from estimated parameter y A −y B yc .A comparison of experimental results with simulation allows to validate the determination of the phase shift introduced by the object.

Simulation of Phase Objects in a Droplet for Detection of Nanoparticles
We performed recently Digital In-line Holography experiments to reconstruct inclusions in a droplet [21,34].In such experiments, the droplet itself induces a magnification factor which is added to the magnification introduced by the imaging system.Using our formalism, the two sides of the droplet and propagation through the droplet can be taken into account as effective optical elements by introduction of their own transfer matrices in the global optical transfer matrices M 1 and M 2 .Let us consider the experimental set-up of Figure 7.The beam emitted by a fibered-laser diode is collimated using a first lens and then focused in the vicinity of a suspended water droplet with a microscope objective.Calibrated microspheres (diameter 20 µm) are introduced in the droplet.Figure 8 (on the left) shows optimal reconstruction of some of the microspheres.The wavelength of the incident beam is 642 nm.The diameter of the water droplet is 1.5 mm and its refractive index is 1.33.The other parameters are e 1 = 42.8 mm, f 1 = 42.8 mm, e 2 = 135.6 mm, f 2 = 5.5 mm, e 3 = 10 mm, z 1 = 0.825 mm, z 2 = 1.125 mm, and z = 9.7 mm.Let us now present in Figure 8 (on the right) the optimal reconstruction issued from a simulated hologram.In the case of this reconstruction, we have considered 7 circular particles: 4 are opaque objects, 3 are transparent phase objects (phase shift introduced π).The fact that the particles are located in a droplet is not a problem.The different optical elements of the droplet are considered through the appropriate definition of matrices M 1 and M 2 .The different inclusions are clearly reconstructed.The opaque and non-opaque objects are easily identified.The possibility to simulate this experiment is particularly attractive for us.This set-up can indeed be used to detect 50 nm nanoparticles [22]: the droplet is seeded with nanoparticles.The droplet is then heated using second laser (a frequency-doubled Nd:YAG laser 4 ns, 10 mJ pulses in our case).Heating a nanoparticle creates a bubble which is reconstructed using Digital In-Line Holography [22].For more details concerning the description of spherical micro-bubbles in this experiment, please refer to [22].This indirect detection method is particularly efficient.With this new simulator that can describe arbitrarily-shaped phase objects as arbitrarily-shaped bubbles, we should be able to obtain information on morphology of nanoparticles through analysis of the shape of microbubbles in the future [35].Our numerical model allows this.

Library of Ice Crystal Objects for the Calibration of Interferometric Airborne Instruments
We develop an airborne instrument to detect and characterize droplets and ice crystals in the atmosphere [36].The technique used is interferometric particle imaging [37][38][39].It has been chosen because fast algorithms can be developed that analyze directly a global image [36].The instrument records on a CCD sensor the out-of-focus image of the scattering particle.In the case of droplets, the out-of-focus image shows regular interference fringes due to interference of both reflected and refracted beams.In the case of irregular rough particles such as ice crystals, the out-of-focus image shows a speckle-like pattern. of the droplet is deduced from the frequency of the fringes, while size and ellipticity of irregular particles are deduced from size and ellipticity of the speck of light in the speckle pattern [40].In this last case, it can be demonstrated that the 2-Dimensional Fourier transform of the speckle pattern gives the 2-Dimensional autocorrelation of the global shape of the object [41].The development of the airborne instrument requires its calibration in laboratory using a second well-known technique.Digital In-line Holography is perfectly appropriate.But the experimental set-up in laboratory requires experiments in cold chambers with different optics and windows.Our matrix-based numerical model allows to describe such environmental conditions by the appropriate definition of the correct optical transfer matrices.For future studies, we have developed a library of the most-common ice crystals that can be observed in the atmosphere.Some exemples are presented in Figure 9: a dendritic crystal with sector-like ends (top left) a malformed crystal (top right), an ordinary dendritic crystal (bottom left) and a crystal with sectorlike branches (bottom right).Using our simulator, Figure 10 shows the digital in-line hologram generated by a transparent particle whose form is the one of a dendritic crystal with sector-like ends (left) and its optimal reconstruction using 2D-FRFT (right).As the object is not opaque, the exact form of the reconstructed can be affected by different effects as for example an intense central peak, or ambiguous limits in different parts of the reconstruction.We think that the models developed should allow a better interpretation of experiments and artifacts encountered in the reconstructions.With the definition of the object presented in the previous section, any kind of object shapes can be simulated.Objects can be opaque or phase objects.Any optical system can be considered by the choice of the appropriate transfer matrices.

Conclusions
We have developed a numerical simulator for digital in-line holography applications.In-line holograms of arbitrarily shaped and arbitrarily located 2D opaque or phase objects are calculated using generalized Huygens-Fresnel integrals.With this formalism, the optical set-up is described in terms of optical transfer matrices.A wide number of spherical or cylindrical lenses can thus be taken into account, which makes the simulator applicable for design and description of in situ experiments.Noise can be added as well, whose spatial frequency can be arbitrarily defined.We discuss future applications of this simulator for detection of nanoparticles in droplets or calibration of airborne instruments that detect and measure ice crystals in the atmosphere.Digital Holography has applications in a wide range of applied sciences from investigation of particles in flows in fluid mechanics, to visualization of cells in biology or medicine, phase contrast metrology.This numerical simulator should thus become a very powerful tool for analysis and simulation of experimental results in a wide range of optical configurations.

Figure 7 .
Figure 7. Experimental set-up for the reconstruction of droplet's inclusions.

Figure 8 .
(a) Experimental reconstruction of inclusions in a droplet; (b) Simulated reconstruction of opaque and transparent inclusions in a droplet.

Figure 9 .
Figure 9.Some examples from our library of transparent phase objects describing ice crystals.A dendritic crystal with sector-like ends (top left); A malformed crystal (top right); An ordinary dendritic crystal (bottom left) and a crystal with sectorlike branches (bottom right).

Figure 10 .
Figure 10.Transparent phase particle whose form is the one of a dendritic crystal with sector-like ends (top); Example of Digital In-line Hologram generated by this particle (bottom left) and its optimal reconstruction using 2D-FRFT (bottom right).