Prof Lixin Guo

Prof Bing Wei

Prof Lu Bai

Prof Yiping Han

Prof. Ruike Yang

Prof Min Zhang

Prof. Hongfu Guo

Profs Renxian Li

Jiajie Wang

Fabrice R A Onofri

Qingwei Duan

Marseille Université on my experimental and theoretical work. In addition, I would like to acknowledge the patience from Claude Rozé and Loïc Méès for being the Monitoring Committee members (CSI) of this thesis and for assessing the progress and quality of the work. Many thanks to Dr. Ruiping Yang and Dr. Saïd Idlahcen of the Coria

List of Figures

The shapes of natural raindrops with increasing radii (from [6]). . . . . . . 6.5 Parameters of the tightly focused elliptical Gaussian beam in Fig. 6.11. . . Table 6.6 Parameters of the elliptical Gaussian beam in Figs. 6.13(a) and 6.14. . . .

XV

Annotation of Symbols

Chapter 1 Introduction

This chapter introduces the research background and significance of the thesis from the aspect of particle scattering. It reviews the research history and current status of related fields, briefly describes the analytical and numerical methods for light scattering, and puts more effort on the introduction of two approximation methods: the geometrical optics approximation (GOA) and the vectorial complex ray model (VCRM), where the incident wave is represented by a collection of discrete rays. Finally, the research content and method, the content of each chapter, the overall framework, the results and the innovations of the thesis are introduced.

Research background and significance

Light scattering by particles is a field of active research with high relevance for science and engineering. It promotes our understanding of the interaction between light and particles and concerns an increasing number of modern technologies, ranging from the measurement and the manipulation of particles to optical communications in turbid media.

Particles are abundant in natural and artificial environments. In many fields of science and engineering, such as combustion diagnosis, atmospheric optics, remote sensing, scattering by interplanetary dust grains, bio-optical imaging, colloidal chemistry, materials science, coating technology and particle sizing technology, one often has to deal with particles ranging from nanoscale to macroscopic.

Since the advent of laser, optical measurement technologies have been favored for its non-contact, accurate and fast merits. The techniques based on elastic light scattering for particle analysis (such as phase Doppler, extinction method, rainbow technology and small particle imaging technology) have found wide applications. The basic principle of these techniques is based on the fact that the intensity distribution, polarization characteristics and spectral characteristics of scattered light are all related to the particles [1,2]. According to the characteristics of the scattered light, a lot of information about the particles can be retrieved, such as the size, shape, speed of motion, temperature, density and compositions.

Efficient theoretical and numerical tools to predict the scattering characteristics of light by various particles are thus essential.

Many models and methods have been developed for analyzing the light scattering by particles. Each method generally has its own range of applicability determined primarily by the particle size relative to the wavelength of incident light. The scattering by those particles that are very small compared to the wavelength can be calculated by the Rayleigh approximation [1,2]. The particles of size comparable to the wavelength lie in the range commonly called the resonance region, where the separation of variables techniques [1,2] and numerical methods [3,4] are usually utilized. And those particles of size much larger than the wavelength can be addressed by the approximation methods, in which the incident wave is represented as a collection of discrete rays, for example, the GOA method [1].

Among the particles of various forms in nature, one group of them are of smooth surface due to surface tension, such as the raindrops [5][6][7][8], droplets in industry [9][10][11][START_REF] Li | Measurement of transient evaporation of an ethanol droplet stream with phase rainbow refractometry and high-speed microscopic shadowgraphy[END_REF][START_REF] Li | Simultaneous measurement of refractive index, diameter and colloid concentration of a droplet using rainbow refractometry[END_REF], underwater bubbles [START_REF] Langley | Critical-angle scattering of laser light from bubbles in water: measurements, models, and application to sizing of bubbles[END_REF][START_REF] Hansen | Mie scattering as a technique for the sizing of air bubbles[END_REF][START_REF] Arnott | Unfolded optical glory of spheroids: backscattering of laser light from freely rising spheroidal air bubbles in water[END_REF][START_REF] Zhang | Influence of bubbles on scattering of light in the ocean[END_REF][START_REF] Mougin | Path instability of a rising bubble[END_REF][START_REF] Lindken | A novel PIV technique for measurements in multiphase flows and its application to two-phase bubbly flows[END_REF][START_REF] Fujiwara | Bubble deformation and flow structure measured by double shadow images and PIV/LIF[END_REF][START_REF] Li | Influence of characteristics of micro-bubble clouds on backscatter lidar signal[END_REF][START_REF] Onofri | Near-criticalangle scattering for the characterization of clouds of bubbles: particular effects[END_REF][START_REF] Cubaud | Bubble dispenser in microfluidic devices[END_REF][START_REF] Yang | Gas content effect on bubble motion[END_REF][START_REF] Oubei | Performance evaluation of underwater wireless optical communications links in the presence of different air bubble populations[END_REF][START_REF] Al-Lashi | A novel high-resolution optical instrument for imaging oceanic bubbles[END_REF][START_REF] Chen | Effects of air bubbles on underwater optical wireless communication[END_REF][START_REF] Wu | Phase critical angle scattering for measurement of transient nanoscale growth rate of a micron-sized bubble[END_REF][START_REF] Selmke | Bubble optics[END_REF], liquid jets [START_REF] Clanet | Transition from dripping to jetting[END_REF][START_REF] Han | Characterization of initial disturbances in a liquid jet by rainbow sizing[END_REF][START_REF] Song | An optimization scheme for the measurement of liquid jet parameters with rainbow refractometry based on Debye theory[END_REF][START_REF] Song | Measurement of temperature gradient in a heated liquid cylinder using rainbow refractometry assisted with infrared thermometry[END_REF][START_REF] Amini | Effect of gravity on capillary instability of liquid jets[END_REF][START_REF] Radev | Review on the instability and optics of capillary jets and glass fibres: a fruitful collaboration between institute of mechanics and IUSTI[END_REF][START_REF] Moallemi | Breakup of capillary jets with different disturbances[END_REF] and so on. As the effects of gravity and air forces become non-ignorable, their geometric shapes are deformed from a sphere or an infinite cylinder, indicating that their scattering characteristics cannot be resolved by the separation of variables techniques. Besides, their size parameters1 usually range from hundreds to thousands or even larger, making the calculation far beyond the capabilities of the existing numerical methods. The developments of alternative approaches allowing fast and accurate (at least asymptotical) calculations are needed for carrying out in-situ characterization of these particles.

For a nonspherical particle of size much larger than the wavelength (size parameter is generally over 100), the approximate methods based on the assumption that the light wave can be represented by a bundle of discrete rays provide feasible ways for calculating the scattering field. Many researchers have contributed to the development of GOA for dealing with spherical or spheroidal particles/bubbles [1,[START_REF] Lock | Ray scattering by an arbitrarily oriented spheroid. I. Diffraction and specular reflection[END_REF][START_REF] Lock | Ray scattering by an arbitrarily oriented spheroid. II. Transmission and crosspolarization effects[END_REF][START_REF] Xu | Extension of geometrical-optics approximation to on-axis Gaussian beam scattering. I. By a spherical particle[END_REF][START_REF] Xu | Extension of geometrical-optics approximation to on-axis Gaussian beam scattering. II. By a spheroidal particle with end-on incidence[END_REF][START_REF] Yu | Geometrical optics approximation of light scattering by large air bubbles[END_REF][START_REF] He | Light scattering by a spheroidal bubble with geometrical optics approximation[END_REF][START_REF] Yu | Model for computing optical caustic partitions for the primary rainbow from tilted spheriodal drops[END_REF], faceted particles [START_REF] Yang | Light scattering by hexagonal ice crystals: comparison of finite-difference time domain and geometric optics models[END_REF][START_REF] Yang | Geometric-optics-integral-equation method for light scattering by nonspherical ice crystals[END_REF] and natural raindrops [6]. In a GOA model, each ray is characterized by its direction, amplitude, polarization and phase. However, for a particle of curved surface, the GOA encounters difficulties or even obstacles in calculating the shift in ray intensity due to the divergence/convergence of wave (divergence factor) and the shift in phase due to focal lines, except for a spherical or spheroidal particle. For this reason, the applications of GOA are under restrictions.

In this context, Ren et al. [START_REF] Ren | Vectorial complex ray model and application to twodimensional scattering of plane wave by a spheroidal particle[END_REF] proposed the vectorial complex ray model (VCRM), a new model for describing the interaction of light rays with a large particle. The distinctive features of VCRM are that: 1) the curvature of wavefront is integrated as an intrinsic property of a light ray, by which one can directly calculate the divergence factor and the phase shift due to focal lines; 2) the ray directions and the Fresnel coefficients are calculated by the wave vectors and their components; 3) the interference of all scattered rays in the whole region is taken into account.

Since the light scattering by a nonspherical particle of size much larger than the wavelength is still a challenging problem in the field of light scattering, this thesis is devoted to the extension of VCRM, with aims to solve the light scattering by infinite cylinders of various forms, and more importantly, to solve the scattered intensity of plane wave and shaped beam by large nonspherical particles in 3D space.

Research status of light scattering by particles

This section summarizes the research status of various calculation methods for the light scattering by particles.

Analytical methods

The search for an exact analytical solution to the scattering field has been traditionally reduced to solving the vector Helmholtz equation for the time-harmonic electric field using the separation of variables techniques in one of the few coordinate systems in which this equation is separable. For this reason, the separation of variables technique are limited to particles of simple forms such as spheres, spheroids, circular or elliptical cylinders.

Lorenz [START_REF] Lorenz | Lysbevaegelsen i og uden for en af plane lysbølger belyst kugle[END_REF] and Mie [START_REF] Mie | Beiträge zur optik trüber medien, speziell kolloidaler metallösungen[END_REF], independently, derived the solution for an isotropic homogeneous sphere in 1890 and 1908. The method that they proposed is now the well-known Lorenz-Mie theory (LMT). The LMT is a rigorous solution of the Maxwell equations and contains all the effects that contribute to the scattering. It was expounded later in depth by Van de Hulst [1] and by Bohren and Huffman [2]. Aden and Kerker provided a solution for the scattering by a coated sphere in 1951 [START_REF] Aden | Scattering of electromagnetic waves from two concentric spheres[END_REF]. Wait in 1955 presented a solution for the scattering by a homogeneous, isotropic, infinite circular cylinder [START_REF] Wait | Scattering of a plane wave from a circular dielectric cylinder at oblique incidence[END_REF]; and it was extended to the scattering by an isotropic, infinite elliptical cylinder by Kim and Yeh in 1991 [START_REF] Kim | Scattering of an obliquely incident wave by a multilayered elliptical lossy dielectric cylinder[END_REF]. Oguchi in 1973 and Asano and Yamamoto in 1975 derived a general solution for the scattering by homogeneous, isotropic spheroidal particles [START_REF] Oguchi | Scattering properties of oblate raindrops and cross polarization of radio waves due to rain-Calculations at 19.3 and 34.8 GHz[END_REF][START_REF] Asano | Light scattering by a spheroidal particle[END_REF].

To give a clear interpretation to the various physical processes in scattering, the Debye series expansion (DSE) method [START_REF] Debye | Das elektromagnetische feld um einen zylinder und die theorie des regenbogens[END_REF] was developed for a homogeneous [START_REF] Hovenac | Assessing the contributions of surface waves and complex rays to far-field Mie scattering by use of the Debye series[END_REF] and coated [START_REF] Lock | Rainbow scattering by a coated sphere[END_REF] sphere respectively in 1992 and 1994. In 2006, Li et al. [START_REF] Li | Debye series for light scattering by a multilayered sphere[END_REF] derived the formula for the DSE in the light scattering by a multilayered sphere and introduced an efficient algorithm permitting stable calculation for a large multilayered sphere. Wu and Li in 2008 [START_REF] Wu | Debye series of scattering by a multi-layered cylinder in an off-axis 2D Gaussian beam[END_REF] pre-sented a simplified but rigorous iterative formula for the DSE in the light scattering by an infinite circular cylinder of multiple layers. In 2010, Shen and Wang [START_REF] Shen | Calculation of Debye series expansion of light scattering[END_REF] provided a stable, reliable and robust algorithm for the DSE calculation of the scattering of plane wave by a uniform sphere in a wide range of sizes and refractive indices. Xu and Lock contributed to the development of DSE for the light scattering by a spheroid [START_REF] Xu | Debye series for light scattering by a spheroid[END_REF], and by a homogeneous or coated nonspherical particle in combination with the extended boundary condition method [START_REF] Xu | Debye series for light scattering by a nonspherical particle[END_REF][START_REF] Xu | Debye series for light scattering by a coated nonspherical particle[END_REF].

With the advent of lasers and their growing applications, the LMT met one of its fundamental limitations, i.e., the assumption that the incident wave must be a plane wave. To address this problem of LMT, in 1988, G. Gouesbet and G. Gréhan et al. presented a theoretical description of the scattering of a Gaussian beam by a homogeneous and isotropic particle of spherical shape, which is now known as the generalized Lorenz-Mie theory (GLMT) [START_REF] Gouesbet | Light scattering from a sphere arbitrarily located in a Gaussian beam, using a Bromwich formulation[END_REF][START_REF] Gouesbet | Generalized Lorenz-Mie Theories[END_REF]. An equivalent approach for the scattering by a dielectric sphere located arbitrarily within a Gaussian beam was presented by Barton et al. [START_REF] Barton | Internal and near-surface electromagnetic fields for a spherical particle irradiated by a focused laser beam[END_REF][START_REF] Barton | Internal fields of a spherical particle illuminated by a tightly focused laser beam: Focal point positioning effects at resonance[END_REF]. In the framework of GLMT, [START_REF] Ren | Radiation pressure forces exerted on a particle arbitrarily located in a Gaussian beam by using the generalized Lorenz-Mie theory, and associated resonance effects[END_REF]1996 derived the radiation pressure (pushing force) and the reverse radiation pressure (pulling force) exerted by Gaussian beam on spherical particles [START_REF] Ren | Radiation pressure forces exerted on a particle arbitrarily located in a Gaussian beam by using the generalized Lorenz-Mie theory, and associated resonance effects[END_REF][START_REF] Ren | Prediction of reverse radiation pressure by generalized Lorenz-References Mie theory[END_REF].

Onofri et al. [START_REF] Onofri | Electromagnetic scattering from a multilayered sphere located in an arbitrary beam[END_REF] presented a solution in the framework of GLMT for the scattering of an arbitrary shaped beam by a multilayered sphere. Ren and Méès et al. [START_REF] Ren | Scattering of a Gaussian beam by an infinite cylinder in the framework of generalized Lorenz-Mie theory: formulation and numerical results[END_REF][START_REF] Méès | Scattering of a Gaussian beam by an infinite cylinder with arbitrary location and arbitrary orientation: numerical results[END_REF] radially inhomogeneous droplet illuminated by an arbitrary shaped beam in a wide range of size parameters [START_REF] Wang | Internal field distribution of a radially inhomogeneous droplet illuminated by an arbitrary shaped beam[END_REF]. More comprehensive reviews of the GLMT can be found in [START_REF] Lock | Generalized Lorenz-Mie theory and applications[END_REF][START_REF] Gouesbet | Generalized Lorenz-Mie theories and description of electromagnetic arbitrary shaped beams: Localized approximations and localized beam models, a review[END_REF].

It is unlikely that these analytical solutions will be significantly extended in the future.

Indeed, the solution for the simplest finite nonspherical particles, spheroids, is already so complex that it behaves like a numerical solution as Mishchenko pointed out in [3].

Numerical methods

Numerical methods provide exact solutions to the scattering of electromagnetic/light wave by nonspherical particles. Most of them fall into two broad categories [3]: the differential equation methods which compute the scattering field by solving the vector wave equation in the frequency or in the time domain; and the integral equation methods which are based on the volume or surface integral counterparts of Maxwell's equations, with exceptions of those hybrid techniques or methods. Here a brief introduction is made to those commonly used numerical methods for the scattering of light by nonspherical particles.

In the discrete dipole approximation (DDA), a scatterer is divided into small cubical subvolumes ("dipoles"), and the interaction between the dipoles are approximated based on the integral equation for the electric field [START_REF] Purcell | Scattering and absorption of light by nonspherical dielectric grains[END_REF][START_REF] Draine | The discrete-dipole approximation and its application to interstellar graphite grains[END_REF][START_REF] Yurkin | The discrete dipole approximation: An overview and recent developments[END_REF][START_REF] Yurkin | The discrete-dipole-approximation code ADDA: Capabilities and known limitations[END_REF]. It can be applied to simulating the light scattering by finite 3D objects of arbitrary geometry.

The finite-difference time domain (FDTD) method solves the Maxwell equations in the time domain by using the finite-difference analog. It was originally developed by Yee in 1966 [START_REF] Yee | Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media[END_REF]. After that, the FDTD method has been extensively applied to solving various electromagnetic problems [START_REF] Taflove | Numerical solution of steady-state electromagnetic scattering problems using the time-dependent Maxwell's equations[END_REF][START_REF] Umashankar | A novel method to analyze electromagnetic scattering of complex objects[END_REF][START_REF] Britt | Solution of electromagnetic scattering problems using time domain techniques[END_REF][START_REF] Furse | Improvements to the finite-difference time-domain method for calculating the radar cross section of a perfectly conducting target[END_REF][START_REF] Wang | FDTD study on wave propagation in layered structures with biaxial anisotropic metamaterials[END_REF][START_REF] Li | Investigation of composite electromagnetic scattering from shiplike target on the randomly rough sea surface using FDTD method[END_REF][START_REF] Cao | Scattering of targets over layered half space using a semi-analytic method in References conjunction with FDTD algorithm[END_REF][START_REF] He | New hybrid FDTD algorithm for electromagnetic problem analysis[END_REF]. The most relevant literatures are [START_REF] Yang | Light scattering by hexagonal ice crystals: comparison of finite-difference time domain and geometric optics models[END_REF][START_REF] Yang | Finite-difference time domain method for light scattering by small ice crystals in three-dimensional space[END_REF], where Yang and Liou et al. proposed the solutions for the light scattering by ice crystal particles.

The T-matrix method may be the most accurate and efficient numerical method for solving the scattering of electromagnetic wave by a nonspherical particle, and can achieve the widest range of size parameter (up to 200 or 300 dependent on the particle shape [START_REF] Bi | Efficient implementation of the invariant imbedding T-matrix method and the separation of variables method applied to large nonspherical inhomogeneous particles[END_REF][START_REF] Bi | A numerical combination of extended boundary condition method and invariant imbedding method applied to light scattering by large spheroids and cylinders[END_REF][START_REF] Doicu | An overview of the methods for deriving recurrence relations for T-matrix calculation[END_REF][START_REF] Zhai | Improvements in the computational efficiency and convergence of the Invariant Imbedding T-matrix method for spheroids and hexagonal prisms[END_REF]).

As one renowned method to compute the T-matrix in the T-matrix formulation of light scattering, the extended boundary condition method (EBCM), or Waterman's T-matrix method, was initially proposed by [START_REF] Waterman | Matrix formulation of electromagnetic scattering[END_REF]1971 [100, 101]. It was later developed by Barber and Hill [START_REF] Barber | Light scattering by particles: computational methods[END_REF], Mishchenko [START_REF] Mishchenko | T-matrix computations of light scattering by large spheroidal particles[END_REF][START_REF] Mishchenko | T-matrix computations of light scattering by nonspherical particles: A review[END_REF][START_REF] Mishchenko | Scattering, absorption, and emission of light by small particles[END_REF] and others. The standard EBCM [START_REF] Mishchenko | Scattering, absorption, and emission of light by small particles[END_REF] is very efficient but will encounter a loss of precision when the particle size becomes larger, the maximum size being affected by the particle aspect ratio. On the other hand, the invariant imbedding method (IIM) was originally proposed by Johnson in 1988 [START_REF] Johnson | Invariant imbedding T matrix approach to electromagnetic scattering[END_REF]. Bi, Sun and Yang et al. [START_REF] Bi | Accurate simulation of the optical properties of atmospheric ice crystals with the invariant imbedding T-matrix method[END_REF][START_REF] Sun | Invariant Imbedding T-matrix Method for Light Scattering by Nonspherical and Inhomogeneous Particles[END_REF] and others have made significant contributions to the development of

IIM. It is based on an electromagnetic volume integral equation and obtains the T-matrix by

growing the scattering volume incrementally in a shell-by-shell manner. The IIM is applicable to particles of relatively large size parameters and extreme aspect ratios. But it is not as efficient as the EBCM due to the large number of differential shells required to discretize the particle volume. Bi and Yang et al. [START_REF] Bi | Efficient implementation of the invariant imbedding T-matrix method and the separation of variables method applied to large nonspherical inhomogeneous particles[END_REF] improved the computational efficiency of IIM by combining the IIM with the LMT and applied to the spheroids and cylinders of size parameters beyond the convergence limit of EBCM. They also proposed a numerical combination of IIM with EBCM [START_REF] Bi | A numerical combination of extended boundary condition method and invariant imbedding method applied to light scattering by large spheroids and cylinders[END_REF], which not only enhanced the efficiency of IIM but also extended the convergence limit of EBCM.

However, a numerical method is severely limited by the particle size. Generally speak-ing, when the size parameter of particle is over 100, the computation using a numerical method becomes rather time-and memory-consuming.

GOA

The geometric optics approximation (GOA) (otherwise known as the ray tracing or ray optics approximation) is an approximate method for the light scattering by large particles, where the incident wave is described by a bundle of light rays. Van de Hulst in his renowned book [1] presented a relatively systematic description of GOA for the scattering of plane wave by a spherical particle, where the direction, amplitude and phase for each scattered ray of perpendicular or parallel polarization were derived. Besides, many researchers have extended and applied the GOA to various scattering problems including:

• The scattering by particles of various forms. Glantschnig and Chen [START_REF] Glantschnig | Light scattering from water droplets in the geometrical optics approximation[END_REF] simplified the calculation for the superposition of the externally reflected rays and the refracted rays with the diffraction and obtained a closed formula for the scattering in the forward angular range (0 • -60 • ) by a water droplet. Adler and Lock et al. [START_REF] Adler | Rainbow scattering by a cylinder with a nearly elliptical cross section[END_REF][START_REF] Lock | Amplification of high-order rainbows of a cylinder with an elliptical cross section[END_REF] examined the scattering by an infinite cylinder of deformed cross section in terms of ray tracing, rainbow angles and the rainbow intensity due to Fresnel coefficients. Hovenac [START_REF] Hovenac | Calculation of far-field scattering from nonspherical particles using a geometrical optics approach[END_REF] developed an algorithm for predicting the far-field scattering from a particle being symmetric about the optical axis. Lock [START_REF] Lock | Ray scattering by an arbitrarily oriented spheroid. I. Diffraction and specular reflection[END_REF][START_REF] Lock | Ray scattering by an arbitrarily oriented spheroid. II. Transmission and crosspolarization effects[END_REF] addressed the reflection, transmission and diffraction of the light rays scattered by an arbitrarily oriented spheroid. Sadeghi et al. [6] extended the GOA to the 3D scattering field near the rainbows produced by natural large raindrops (neither spheres nor spheroids). Yang and Liou [START_REF] Yang | Geometric-optics-integral-equation method for light scattering by nonspherical ice crystals[END_REF] proposed a new geometrical optics model which used the ray-tracing technique to solve the near field on the ice crystal surface and then transformed the near field to the far field based on the electromagnetic equivalence theorem. Bi et al. [START_REF] Bi | Assessment of the accuracy of the conventional ray-tracing technique: Implications in remote sensing and radiative transfer involving ice clouds[END_REF] assessed the uncertainties with the conventional geometrical optics in remote sensing and radiative transfer simulations.

• The scattering by bubbles. In 2008, Yu et al. [START_REF] Yu | Geometrical optics approximation of light scattering by large air bubbles[END_REF] studied the scattering of plane wave by a spherical bubble based on GOA, where the total reflection effect was taken into account to improve the calculation accuracy. In 2012, He et al. [START_REF] He | Light scattering by a spheroidal bubble with geometrical optics approximation[END_REF] studied the scattering of plane wave by a spheroidal bubble with end-on incidence, where the effects of the size and the aspect ratio of bubble and the width of incident beam on the scattering patterns were analyzed. In 2016, Sentis and Onofri et al. [START_REF] Sentis | Scattering of light by large bubbles: Coupling of geometrical and physical optics approximations[END_REF] improved the GOA by combining it with the physical optics approximation (POA) in modeling the scattering properties of large spherical bubbles, where the interference between higher-order rays, the Goos-Hänchen shift, the tunneling phase and the weak caustic associated with the critical angle were taken into account.

• The scattering by coated particles. Lock et al. [START_REF] Lock | Rainbow scattering by a coated sphere[END_REF] presented an intuitive interpretation to the first-order rainbow from a coated sphere by using ray theory in 1994. In 2004, Xu et al. [START_REF] Xu | Geometrical-optics approximation of forward scattering by coated particles[END_REF] provided an algorithm by means of GOA for calculating the scattered intensity by a coated spherical particle in the forward angular range (0 • -60 • ). In 2014, Zhai et al. [START_REF] Zhai | Coated sphere scattering by geometric optics approximation[END_REF] examined the scattering processes of a coated sphere with the GOA method. They parameterized the light rays interacting with a coated sphere and simplified the calculation for those terms of degeneracy paths and repeated paths.

• Ray tracing in inhomogeneous media. The path of a light ray in inhomogeneous media is curved. In the 1960s, Montagnino [START_REF] Montagnino | Ray tracing in inhomogeneous media[END_REF] and Marchand [START_REF] Marchand | Ray tracing in gradient-index media[END_REF] contributed to the development for the methods for tracing the ray paths in inhomogeneous media. In 1982, Sharma et al. [START_REF] Sharma | Tracing rays through graded-index media: a new method[END_REF] proposed an efficient method for tracing the light rays passing through gradient-index media. Later, Sharma achieved a method for computing the optical path in gradient-index media, which proved to be more accurate and faster after a comparison with other methods [START_REF] Sharma | Computing optical path length in gradient-index media: a fast and accurate method[END_REF]. As for the spherical particles of inhomogeneous media, Li et al. [START_REF] Li | Geometrical-optics approximation of forward scattering by gradientindex spheres[END_REF] presented a solution in 2007 to the scattered intensity of plane wave by a gradient-index sphere in the forward angular range (0 • -60 • ). In 2008 and 2015, Lock et al. [START_REF] Lock | Scattering of an electromagnetic plane wave by a Luneburg lens. I. ray theory[END_REF][START_REF] Lock | Scattering of a plane electromagnetic wave by a generalized Luneburg sphere-Part 1: Ray scattering[END_REF] analyzed the scattering of plane wave by a Luneburg lens.

• The scattering by chiral particles. In 2015, the rainbow angles of a chiral sphere were calculated using ray tracing by Wu et al. [START_REF] Wu | Geometric optics interpretation for rainbow scattering of a chiral sphere[END_REF]. In 2019, the scattering pattern of plane wave by a chiral sphere in the forward direction (0 • -90 • ) has been solved with GOA by Lu et al. [START_REF] Lu | Geometrical optics approximation for forward light scattering by a large chiral sphere[END_REF].

• The scattering of shaped beam. In 2006, Xu et al. extended the GOA for calculating the scattered intensity of circular Gaussian beam by a sphere [START_REF] Xu | Extension of geometrical-optics approximation to on-axis Gaussian beam scattering. I. By a spherical particle[END_REF] and by a spheroid with end-on incidence [START_REF] Xu | Extension of geometrical-optics approximation to on-axis Gaussian beam scattering. II. By a spheroidal particle with end-on incidence[END_REF], where the Gaussian beam was described by a bundle of light rays of different propagation directions, amplitudes and phases. To calculate the optical forces of an arbitrary shaped beam, Shao et al. in 2019 decomposed an arbitrary light beam into plane waves. The intensity of each plane wave was dependent on the Fourier angular spectrum; and each plane wave was further represented by a bundle of light rays parallel in propagation direction [START_REF] Shao | Calculation of optical forces for arbitrary light beams using the Fourier ray method[END_REF].

• The scattering by absorbing particles. In 2005, Chang [START_REF] Chang | Ray tracing in absorbing media[END_REF] detailed a ray tracing method to characterize the plane wave propagating in a lossy media. Yu et al. in 2009 presented a method based on GOA to address the light scattering by an absorbing sphere, where the effective refractive index and the effective refractive angle were introduced, and the formulas for the phase shifts due to reflections and refractions were derived [START_REF] Yu | Geometrical optics approximation for light scattering by absorbing spherical particles[END_REF]. In 2018, Lindqvist et al. [START_REF] Lindqvist | Ray optics for absorbing particles with application to ice crystals at near-infrared wavelengths[END_REF] derived a ray-optics solution which took into account the inhomogeneous nature of the represented wave inside an absorbing particle and applied it to the light scattering by ice particles in the near-infrared wavelengths.

• The coupling with physical optics. In 2015, Huang et al. combined the GOA and the concept of divergence factor with the physical optics (PO) and the physical theory of diffraction for the scattering by perfectly electrical conducting targets [START_REF] Huang | GO/PO and PTD with virtual divergence factor for fast analysis of scattering from concave complex targets[END_REF]. Sentis et al.

in 2016 coupled the GOA with the POA and developed an improved version of GOA, which allows one to predict the scattering pattern by large bubbles with high accuracy [START_REF] Sentis | Scattering of light by large bubbles: Coupling of geometrical and physical optics approximations[END_REF].

• The application to optical forces. Ashkin in 1992 calculated the forces of single-beam gradient radiation pressure on micron-sized dielectric spheres in the ray optics regime [START_REF] Ashkin | Forces of a single-beam gradient laser trap on a dielectric sphere in the ray optics regime[END_REF]. Zhou et al. in 2012 calculated the optical forces on a triaxial ellipsoid by vectorial ray tracing [START_REF] Zhou | Calculation of optical forces on an ellipsoid using vectorial ray tracing method[END_REF]. Later Shao et al. in 2019 combined the Fourier optics and ray optics for calculating the optical force of an arbitrary shaped beam exerted on a spherical particle [START_REF] Shao | Calculation of optical forces for arbitrary light beams using the Fourier ray method[END_REF].

• The application to inelastic scattering. The theoretical treatment of inelastic scattering by wave theory is rather complicated and, most of the calculations are limited to a spherical particle of homogeneous medium or with inclusions of size smaller than the wavelength.

Being flexible and efficient in computation, the ray optics was utilized to calculate the inelastic scattering by a large spherical particle [START_REF] Zhang | Hybrid inelastic-scattering models for particle thermometry: polarized emissions[END_REF][START_REF] Velesco | Geometrical optics calculation of inelastic scattering on large particles[END_REF], by a particle of any shape where the particle surface was described by triangles [START_REF] Weigel | Inelastic scattering by particles of arbitrary shape[END_REF], and by a spherical particle with inclusions [START_REF] Weigel | Inelastic scattering on particles with inclusions[END_REF].

A GOA method provides clear insight into the reflection and refraction processes and has advantage in the cases where the exact numerical methods are hard to achieve or no rigorous theory exists, as concluded by Xu et al. [START_REF] Xu | Extension of geometrical-optics approximation to on-axis Gaussian beam scattering. I. By a spherical particle[END_REF]. However, a careful reader may notice that most of the GOA methods available today have addressed only the scattered intensity by spherical or spheroidal particles, or faceted ice crystals. This is mainly due to the fact that although the GOA is flexible in principle, it usually encounters difficulties or even obstacles in accounting for the divergence factor and the phase shift due to focal lines, which restricts the precision and application of GOA.

VCRM

To resolve these difficulties encountered in GOA, Ren et al. in 2011 proposed the vectorial complex ray model (VCRM) [START_REF] Ren | Vectorial complex ray model and application to twodimensional scattering of plane wave by a spheroidal particle[END_REF]. In VCRM, the curvature of wavefront is integrated as an intrinsic property of a light ray, by which the divergence factor and the phase shift due to focal lines can be calculated directly. Furthermore, the ray directions and the Fresnel coefficients are calculated by the wave vectors and their components, which avoids the tedious calculation when a nonspherical particle is involved. Although it is a high-frequency approximation model, the VCRM provides one feasible way for calculating the scattered intensity of plane wave or shaped beam by a large particle of arbitrarily smooth surface2 .

The validity of VCRM has been examined numerically by comparing with the multilevel fast multipole algorithm (MLFMA) for the light scattering by a large ellipsoidal particle [START_REF] Yang | Comparison of scattering diagrams of large non-spherical particles calculated by VCRM and MLFMA[END_REF]. Furthermore, Onofri et al. experimentally examined the VCRM for the light scattering by a large oblate droplet trapped in an acoustic field [START_REF] Onofri | Experimental validation of the vectorial complex ray model on the inter-caustics scattering of oblate droplets[END_REF]. However, the VCRM is still in its early stage of development. The existing numerical implementations of VCRM address only the scattering in a few simple configurations. One of them is that all light rays should propagate in a two-dimensional (2D) plane, for example, in one of the symmetric planes of a spheroidal particle.

In the past few years, several researchers have been trying to find a solution to the 3D scattered intensity by large nonspherical particles in the framework of VCRM. Recently, in Yang's doctoral thesis which was defended in the December of 2019 [START_REF] Yang | Numerical simulation of light scattering of a pendent droplet by statistic vectorial complex ray model[END_REF], an alternative implementation of VCRM using statistic approach has been achieved for the 3D scattered intensity by a pendent droplet. But, the statistic approach requires a huge number of photons, at least 10 8 to 10 10 to obtain an acceptable result for a pendent droplet of 2 mm. This statistic approach, though effective, is so time-and memory-consuming that it violates the efficient merit of VCRM.

Main content and frame arrangement

Because of the flexible and efficient merits of VCRM, this thesis is devoted to the extension of VCRM to the light scattering by an infinite cylinder of arbitrarily smooth cross section, and more importantly, to the scattered intensity of plane wave and shaped beam by a large nonspherical particle in 3D space. This thesis focuses mainly on the following issues:

1. The extension of VCRM to the 2D scattering of plane wave by an infinite cylinder of any smooth cross section.

2. The calculation method based on VCRM for the 3D scattered intensity by a large particle of any smooth surface.

3. The application of the proposed method to calculating and analyzing the scattered intensity of plane wave by a real liquid jet of complex surface in 3D space.

4. The examination of the proposed method by experiments.

5. The extension to the 3D scattered intensity of shaped beam, an elliptical Gaussian beam for example.

The scattering by inhomogeneous particles is not involved here, because the variation of wavefront in inhomogeneous media remains unknown.

The thesis includes seven Chapters. The content of each Chapter is outlined as follows:

Chapter 1 introduces the research background and the significance of this thesis, reviews the research history of related fields, points out the issues that this thesis aims to solve, introduces the main content and the frame arrangement, and outlines the research results and innovations.

Chapter 2 discusses in depth the classical GOA method for the scattering of plane wave by circular cylinders and spheres.

Chapter 3 introduces the interaction of light rays with a cylinder in the framework of VCRM and, extends the VCRM to the scattering of plane wave by an infinite cylinder of any smooth cross section. The proposed algorithm is applied to solving the full diagram of the scattered intensity by a cylinder whose cross section takes the shape of a natural raindrop.

The effects of shape deformation, refractive index and the direction of incident wave on the scattering patterns are investigated and quantitatively analyzed.

Chapter 4 proposes the calculation method for the 3D scattering of plane wave by a large particle of any smooth surface. In the framework of VCRM, the ray tracing, divergence factor, phase shifts due to focal lines and optical path, and cross polarization in 3D scattering are addressed by an elegant way using vectorial rays and wave-front curvature. To account for the superposition of scattered ray in 3D space, a triangulation-based interpolation method is proposed, thus breaking through the bottle-neck problem for VCRM in calculating the 3D scattered intensity.

Chapter 5 realizes the calculation for the 3D scattered intensity of plane wave by a real liquid jet. The geometric model of a continuous water jet near the nozzle is established through image edge detection and data fitting. By applying the calculation method proposed in Chapter 4, the 3D scattering field by the jet is successfully simulated. Emphasis is put on the 3D intensity distribution near the first-order and second-order rainbows. The effect of the stream-wise curvature of jet surface on the scattering characteristics is analyzed. The difference from the scattering field by an infinite cylinder is discussed. An experimental setup is also established for measuring the 3D scattered intensity by a liquid jet, and the result by simulation is examined by the experiment.

Chapter 6 proposes a ray description method in the framework of VCRM for incident elliptical Gaussian beams. After the validation of the proposed model by comparing with the GLMT for spherical particles, the 3D scattered intensity of elliptical Gaussian beam by a real liquid jet which has a complex geometry has been successfully calculated. The scattering characteristics of the elliptical Gaussian beams of different divergence angles are investigated. The spatial characteristics of the scattering field of a tightly focused elliptical Gaussian beam is analyzed.

Chapter 7 concludes the work of the current thesis and gives perspectives for the future.

The innovations

This thesis addresses the 2D and 3D scattered intensity by a large particle of arbitrarily smooth surface in the framework of VCRM. The innovations are as follows:

1. Based on the VCRM, an algorithm is proposed for solving the light scattering by an infi-nite cylinder of arbitrarily smooth surface in the high-frequency limit. For the scattering by an infinite cylinder, the rigorous LMT is limited to the cylinders which have circular or elliptical cross sections; while the GOA method, except for a circular or elliptical cylinder, is hard to calculate the divergence factor and the phase shift due to focal lines. The algorithm proposed in this thesis allows one to obtain the direction, amplitude, phase, polarization and divergence factor for each scattered ray, and the final interference intensity of the scattered rays from an infinite cylinder of arbitrarily smooth surface. It provides an effective tool for calculating and analyzing the scattering patterns by a cylinder of cross section ranging from simple to complex.

2. The bottle-neck problem of VCRM in calculating the 3D scattered intensity has been solved. Inspired by the idea of triangulation in the area of Computer Graphics, the direction of each 3D scattered ray is represented as a vertex in the (φ, ψ) space, where φ and ψ are respectively the azimuth angle and the elevation angle. These vertexes which represent the directions of scattered rays are firstly meshed by triangles. Then, inside each generated triangle, interpolation is carried out for the amplitudes and phases at the grid points which are enclosed by the triangle, according to the amplitudes and phases of the three vertexes (three scattered rays). The overlapping triangles at a grid point (φ i , ψ j ) account for the interference effect in the scattering direction of azimuth angle φ i and elevation angle ψ j . By introducing this triangulation-based interpolation method, the calculation for the superposition of the scattered rays in 3D space can be achieved. Besides, the ray tracing, divergence factor, phase shifts due to focal lines and optical path, and cross polarization encountered in 3D scattering are addressed in the framework of VCRM. It provides an effective and efficient approach to the scattered intensity in 3D space by a large smooth particle of any shape.

3. The scattered intensity of plane wave in 3D space by a real liquid jet has been successfully simulated and experimentally validated. The study of the light scattering by a liquid jet is motivated primarily by the needs for developing optical means to characterize the jet size and refractive index (temperature). However, the scattering pattern of a real jet can hardly be calculated using the existing analytical theories or numerical methods because of the jet's complex geometry and large size (compared to wavelength). In this thesis, the 3D far-field scattered intensity of plane wave by a real liquid jet has been solved.

It is found that due to the stream-wise curvature of jet surface, the scattered rays of different orders by a real jet are naturally separated in the 3D space, leading to scattering patterns that have rarely been observed. Besides, an experiment has also been carried out.

The agreement between the simulated 3D scattered intensity and the experimental result validates the proposed algorithm which is based on VCRM while allows to calculate 3D scattering field.

4. Many researchers have contributed to the development and application of the GLMT for the scattering of a Gaussian beam by particles. However, the GLMT is applicable only to particles of simple forms so that the separation of variable method could be carried out.

In this thesis, the ray model for incident elliptical or circular Gaussian beam is proposed in the framework of VCRM, thus making it possible to calculate the scattered intensity of laser beam by a large particle of any smooth surface. As an example, the 3D far-field scattered intensity of elliptical Gaussian beam by a real liquid jet has been achieved.

The scattering characteristics near the first-order and second-order rainbows for incident elliptical Gaussian beams of different divergence angles are investigated in 3D space.

Chapter 2 Light scattering by a circular cylinder or a sphere based on GOA

Overview

The scattering of plane wave by an infinite circular cylinder under normal incidence or by a spherical particle is one of the basic problems in the filed of light scattering [1]. When the transversal size of cylinder or the radius of sphere is much larger than the wavelength of light, the geometrical optics approximation (GOA) method provides an asymptotic solution to the scattering field. Although it is not as exact as the LMT or DSE, the calculation by GOA is much simpler and more efficient. In this chapter, the GOA method for calculating the scattered intensity by a circular cylinder under normal incidence and by a sphere are discussed. It provides the fundamental conceptions of the interaction of light rays with a particle, thus facilitating the understanding of the VCRM and its extensions in the following chapters.

Light scattering by an infinite circular cylinder

Calculation method

Consider an infinite circular cylinder illuminated by a plane wave as shown in Fig. 2.1.

The Cartesian coordinate system Oxyz is set such that the z axis is along the cylinder axis and the x axis coincides with the direction of the incident plane wave. The transversal radius of the cylinder is denoted as a. The refractive indices of the surrounding medium and the cylinder are 1 and m, respectively. The scattering angle θ denotes the direction of a scattered ray and is counted in the xy plane from the x axis.

The incident plane wave is regarded as a bundle of light rays of same direction, equal amplitude and equal phase. For an incident light ray, it is subjected to continual reflections and refractions by the cylinder surface. Thus, theoretically speaking, there is an infinite series of emergent rays. An externally-reflected ray is of order p = 0, a transmitted ray without internal reflection is of p = 1, while an emergent ray which undergoes p -1 times of internal reflections is a p-order scattered ray as illustrated in Fig. 2.2. Such a definition for the order of a scattered ray gives certain convenience in later calculation.

For the convenience of calculation, the complementary angles of the incident angle α and the refracted angle β are introduced as τ and τ ′ , respectively. τ = 0 indicates the grazing incidence while τ = π/2 indicates the central incidence. For one of the emergent rays that originate from the same incident ray, its angle with the cylinder surface is always τ. For any light ray inside the cylinder which originates from this incident ray, its angle with the cylinder surface is a constant of value τ ′ as shown in Fig. 2.2.

From Fig. 2.2, one can derive that the total deviation angle of a p-order emergent ray from the incident direction (x axis) is

θ ′ = 2τ -2pτ ′ (2-1)
In practice, the scattered light is observed in [0, 2π], or in the region of [0, π] considering the symmetry of the scattering field. The scattering angle θ in the interval [0, π] is related to the total deviation angle by

θ ′ = c2π + qθ (2-2)
where c is an integer equal to the times that the emergent ray has crossed the x axis. q is equal to +1 or -1 ensuring that θ is in the interval [0, π].

At the first reflection, the ratio of the reflected amplitude to the incident amplitude is calculated according to the Fresnel reflection coefficients:

r ⊥ = sin τ -m sin τ ′ sin τ + m sin τ ′ (2-3) r ∥ = m sin τ -sin τ ′ m sin τ + sin τ ′ (2-4)
where the subscript ⊥ indicates the polarization perpendicular to the scattering plane xy, while ∥ indicates the polarization parallel to the scattering plane xy. These two polarizations are respectively along ê⊥ and ê∥ as shown in Fig. 2.1.

At a reflection inside the circular cylinder, τ and τ ′ are reversed compared to the first reflection. Thus, the reflection coefficients at an internal reflection are -r ⊥ and -r ∥ . According to the conservation of energy, the refraction coefficient t X and the reflection coefficient r X are related by

t X 2 m sin τ ′ sin τ + r X 2 = 1 (2-5)
where X =⊥ or ∥. Thus, the refraction coefficient (always positive) can be expressed as

t X = ( 1 -r X 2 ) 1/2 ( sin τ m sin τ ′ ) 1/2 (2-6)
For a scattered ray of order p ≥ 1, it has undergone p -1 times of internal reflections and two times of refractions: one refraction from air to cylinder and one refraction from cylinder to air. The refraction coefficients at the two refractions are denoted as t X and t ′ X , respectively. According to Eq. (2-6) and the fact that τ and τ ′ are reversed at the two refractions, we have

t X t ′ X = ( 1 -r X 2 ) (2-7)
Then, the variation of the amplitude of a p-order scattered ray due to the reflection and refraction coefficients is obtained as

ε X,p =        r X if p = 0 (1 -r X 2 )(-r X ) p-1 if p ≥ 1 (2-8)
where the term (1r X 2 ) corresponds to the two refractions, while (-r X ) p-1 the p -1 times of internal reflection. On the other hand, a reflected/refracted ray pencil might be converging or diverging and its intensity will be higher or lower accordingly. Consider a ray pencil of width dl in the xy plane, thickness dz along the axis of cylinder and uniform light intensity I 0 . It illuminates the surface of cylinder in an area of dsdz where the complementary angle of incident angle varies from τ to τ + dτ. The incident energy S inc is I 0 dldz. By simple geometric principles, one can deduce that dl is equal to sinτds where ds = adτ as shown in Fig. 2.3. Thus, the light energy incident in the area dsdz can also be expressed as

S inc = I 0 a sin τdτdz (2-9)
The incident energy is then subjected to continual divisions because of refractions and reflections. For a p-order scattered ray pencil of polarization X, the light energy is reduced to

S X,p = S inc ε 2 X,p (2-10)
Suppose the scattering angle of this scattered ray pencil ranges from θ to θ + dθ as shown in Fig. 2.3. The transversal area of the scattered ray pencil, at a distance r far from the emergent point, is denoted as ds ′ dz. In fact, the transversal area of the incident ray pencil is so small that the width dl tends to be zero. Thus, the area of the scattered ray pencil on the cylinder can be regarded as zero, i.e. s 0 → 0 (s 0 is illustrated in Fig. 2.3). Then, the arc length of the scattered ray pencil ds ′ is approximated to be

ds ′ = rdθ (2-11)
Then, the light intensity of the scattered ray pencil observed at r and of polarization X can be obtained by I X,p = S X,p ds ′ dz , that is

I X,p = I 0 a sin τdτdzε 2 X,p rdθdz = a r I 0 ε 2 X,p D (2-12)
The coefficient D is called the divergence factor of a light ray and is defined for an infinite circular cylinder as

D = sin τ dθ/dτ (2-13)
From Eqs. (2-1) and (2-2), we have dθ/dτ = ±dθ ′ /dτ. Considering that the divergence factor D should be positive, dθ/dτ is calculated to be

dθ dτ = dθ ′ dτ = 2 |m sin τ ′ -p sin τ| m sin τ ′ (2-14)
By omitting the term I 0 in Eq. (2-12) and then multiplying by kr, the dimensionless scattered intensity in the far field is expressed by

I X,p = kaDε 2 X,p (2-15)
where k is the wave number of the scattered ray (k = 2π/λ). The dimensionless term ka is the size parameter which depicts the size of particle relative to the wavelength of light.

For an incident beam of coherent length much larger than the size of particle, the inter-ference between different scattered rays should be taken into account. Besides the amplitude discussed in the preceding part, the calculation for the interference field needs also the information about the phase for each scattered ray.

For a particle illuminated by plane wave, there are three factors that affect the phase of a scattered ray: reflection, optical path and focal lines. The time factor is chosen as exp(iωt).

(1) The phase shift due to reflections ϕ X,r . According to the Fresnel formulas, the refraction coefficients are always positive, which means the refraction does not alter the phase; while the reflection coefficients may be negative, indicating the reflection may change the sign of the amplitude and thereby introduce a phase shift of π (-1 = e iπ ).

On the other hand, when total reflection occurs, all the energy of light is reflected and the Fresnel reflection coefficients r ⊥ and r ∥ are equal to exp(iδ ⊥ ) and exp(iδ ∥ ), respectively, where the phase shifts δ ⊥ and δ ∥ are given by

               δ ⊥ = 2 arctan ( √ sin 2 α-M 2 t cos α ) δ ∥ = 2 arctan ( √ sin 2 α-M 2 t M 2 t cos α ) (2-16)
where M t is the refractive index of the refracted medium relative to the incident medium.

A scattered ray might have undergone one or more times of reflection. Since these reflection coefficients, being positive or negative, real or imaginary, are already included in the factor ε X,p defined in Eq. (2-8), the phase shift of a p-order scattered ray due to reflections, ϕ X,r , can be retrieved directly by the argument (phase) of ε X,p :

ε X,p = ∥ε X,p ∥ exp(iϕ X,r ) (2-17) (2)
The phase shift due to optical path ϕ p,OP . The phase shift ϕ p,OP of a scattered ray is caused by its optical path, usually compared to a reference ray. The reference ray, free of the refraction by the particle, arrives at the center of particle in the the direction of the incident ray and then emerges in the same direction as the scattered ray.

As shown in Fig. 2.4, the externally-reflected ray (p = 0) has a shorter optical path than the reference ray, thus it has a positive phase when compared to the reference ray: where σ equals a sin τ. The scattered rays of p ≥ 1 have longer optical paths than the reference rays. For the scattered ray of p = 1, the phase shift due to optical path is

ϕ 0,OP = 2π λ 2σ (2-18) W ' W W a ' W s s O 0 p 1 p 2 p V V V V Fig 2.
-2π
λ ms, where the negative sign indicates that the phase lags for a longer optical path. s is the geometric length between two successive interactions, equal to 2a sin τ ′ . Then, its phase relative to the reference ray is calculated by

ϕ 1,OP = - 2π λ ms -(- 2π λ 2d) = 2π λ (2a sin τ -2ma sin τ ′ ) (2-19)
Then, for the scattered ray of order p, the phase shift relative to the reference ray is deduced as

ϕ p,OP = 2π λ 2a(sin τ -pm sin τ ′ ) (2-20) (3
lines is p -(1 -s)/2
, where p is the order of the scattered ray and s takes the sign of dθ ′ /dτ. Thus, the phase shift due to focal lines in the case of an infinite circular cylinder is

ϕ p,FL = π 2 [ p -(1 -s)/2 ] (2-21)
The total phase of a scattered ray of polarization X (X =⊥ or ∥) and order p is the summation of these three kinds of phase shifts:

ϕ X,p = ϕ X,r + ϕ p,OP + ϕ p,FL (2-22)
It should be noted that the phase shifts due to optical path and focal lines, ϕ p,OP and ϕ p,FL , are independent of polarization.

Finally, according to Eqs. (2-15) and (2-22), the complex amplitude of a p-order scattered ray of polarization X can be obtained as

E X,p = √ π 2 kaDε X,p exp(iϕ X,p ) (2-23)
The coefficient √ π/2 is added so that the amplitude calculated by GOA is consistent with that by LMT [START_REF] Jiang | Theoretical study of light scattering by an elliptical cylinder[END_REF].

Since this thesis is devoted to the scattering of light, an approximate model is used for the diffraction component [2,[START_REF] Jiang | Scattering from an elliptical cylinder by using the vectorial complex ray model[END_REF][START_REF] Jiang | Theoretical study of light scattering by an elliptical cylinder[END_REF]:

E d = 1 + cos θ 2 sin(ka sin θ) sin θ exp ( i 3π 2
) (2-24)
The scattering field at the scattering angle θ is then the summation of the complex amplitudes of those rays scattered in this direction and the corresponding diffraction component:

E X (θ) = ∞ ∑ p=0 E X,p (θ) + E d (θ) (2-25)
In practice, the upper limit of p is usually set as 7 because the scattered rays of a higher order p have much lower intensity. Then, the scattered intensity in the scattering direction θ is computed by

I X (θ) = E X (θ) E * X (θ) (2-26)

Results and discussions

Since the GOA method is valid only for a particle of size much larger than the wavelength, it is necessary to examine its applicable domain before its application. In Fig. 2.6, the scattered intensity calculated by the GOA method described in the preceding section is compared with LMT for three circular cylinders of same refractive index but different size parameters. Here, the perpendicular polarization is concerned. From the comparison, it may conclude that:

• The prediction for light scattering based on the GOA method is satisfying only when the size parameter 2πa/λ is over 100.

• The larger the size parameter, the better the agreement of GOA with LMT will be.

• When 2πa/λ > 100, the discrepancy of GOA from LMT exists mainly near the caustics (the rainbow angles for example) and in the region where surface wave is important.

Although the deficiency in predicting the light intensity near caustics and the effect of surface waves, the GOA method still captures most part of the scattering diagram for a particle of size parameter over 100.

The GOA method allows to distinguish the contributions from different orders of scattered rays and thus provides a clear interpretation to the scattering pattern. Considering that the first-order rainbow will be repeatedly mentioned in the following chapters, a detail discussion on the formation mechanism and the intensity distribution of the first-order rainbow from a circular cylinder is given here. and refractive index [11,[START_REF] Wu | Phase critical angle scattering for measurement of transient nanoscale growth rate of a micron-sized bubble[END_REF][START_REF] Han | Study of refractometry of rainbow and applications to the measurement of instability and temperature gradient of a liquid jet[END_REF][START_REF] Wu | Phase rainbow refractometry for accurate droplet variation characterization[END_REF].

It should be noted that at the rainbow angle, the scattered intensity is incorrectly predicted as infinity, while beyond the rainbow angle the scattered intensity of p = 2 rays drops abruptly to zero. This abrupt change of intensity arises from the intrinsic defect of ray optics. To tackle with this problem, one feasible way is by taking into account the diffraction effect in the region where light intensity shows discontinuity. To improve this defect of ray optics needs much more effort and is beyond the scope of current thesis.

Light scattering by a spherical particle

Compared to the calculation method for the scattering by an infinite circular cylinder as given in Sec. 2.2.1, the divergence factor D and the numbers of focal lines are different in the scattering by a spherical particle. This is because an infinite circular cylinder has a curvature equal to zero along the cylinder axis, while a sphere does not. This section covers the solution of GOA to the scattering of plane wave by a spherical particle. For concise, those calculation processes which are identical to the ones for a circular cylinder will not be repeated here.

Calculation method

Because of the spherical symmetry, the scattered light in one meridional plane is exactly the same as the scattered light in another meridional plane. Thus, one needs only the calculation for the scattered light in one of the meridional planes. In one meridional plane of the spherical particle, the ray paths are same with the paths in the cross section of a circular cylinder, if the circular cylinder and the sphere have identical radius and refractive index.

For a p-order ray of polarization X (X =⊥ or ∥) scattered by a sphere of radius a, its intensity at a distance r from the sphere is given by [1]

I X,p = a 2 r 2 Dε 2 X,p I 0 (2-27)
where the divergence factor D is defined by

D = sin τ cos τ sin θdθ/dτ = m sin τ sin τ ′ cos τ 2 sin θ|m sin τ ′ -p sin τ| (2-28)
These symbols in Eqs. and (2-28) have been explained in the preceding Sec. 2.2.1.

By omitting the term I 0 in Eq. (2-27) and then multiplying by (kr) 2 , the dimensionless scattered intensity in the far field by the sphere is expressed as

I X,p = (ka) 2 Dε 2 X,p (2-29)
where k is the wave number of the scattered ray (k = 2π/λ). The amplitude formula is then obtained as

E X,p = ka √ Dε X,p (2-30)
Besides the focal lines of type A as for an infinite circular cylinder, the focal lines of type B in Fig. 2.5 also exist for a spherical particle. In fact, the focal lines of type A are caused by the convergence of wave in the meridional plane, while the focal lines of type B arise from the convergence of wave in the direction perpendicular to the meridional plane.

For a sphere, the number of the type A focal lines for a p-order scattered ray is given by p -(1s)/2, where s takes the sign of dθ ′ /dτ as has been discussed for a circular cylinder.

While the number of the type B focal lines is given by -2c + (1q)/2 [1], where c and q have been discussed in Eq. (2-2).

Thus, the total phase shift due to focal lines in the case of a spherical particle is

ϕ p,FL = π 2 [ p -2c + 1 2 s - 1 2 q ] (2-31)
The phase shifts due to reflection and optical path are the same as those for a circular cylinder and will not be repeated here. Then, by substituting the value of ϕ p,FL into Eq.

(2-22), the total phase ϕ X,p of a light ray scattered by a spherical particle can be obtained.

The complex amplitude of a p-order scattered ray of polarization X in the scattering by a spherical particle is then expressed as

E X,p = ka √ Dε X,p exp(iϕ X,p ) (2-32)
Since it is not within the scope of this thesis to discuss the diffraction component, the forward diffraction is simply calculated by applying Babinet's principle. The amplitude of the forward diffraction of plane wave by a sphere is given by [START_REF] Lock | Ray scattering by an arbitrarily oriented spheroid. I. Diffraction and specular reflection[END_REF][START_REF] Ren | Vectorial complex ray model and application to twodimensional scattering of plane wave by a spheroidal particle[END_REF]]

E d = (ka) 2 J 1 (kaθ) kaθ (2-33)
where J 1 () is the Bessel function of the first kind of order 1. A more sophisticated model for the forward diffraction can be found in [START_REF] Verrier | Digital in-line holography in thick optical systems: application to visualization in pipes[END_REF][START_REF] Shen | Diffraction of a plane wave by an infinitely long circular cylinder or a sphere: solution from Mie theory[END_REF].

Results and discussions

In Fig. 2.9, the scattered intensity calculated by the GOA is compared with LMT for three homogeneous spherical particles of different radii. The relative refractive index m is 1.3322 and the wavelength of incident plane wave λ is 0.6328 µm. Here the perpendicular polarization is concerned. Since the comparison between GOA and LMT has been made for circular cylinders of size parameters ranging from 50 to 300 in Fig. 2.6, the comparison for the spherical particles will focus on larger size parameters.

One can see that the agreement of GOA with the rigorous LMT is generally satisfying.

And, as the radius a is increased from 50 µm to 300 µm, the agreement is getting better. One notable discrepancy is near the first-order and second-order rainbow angles (129.3 • < θ < 137.8 • ). The incorrect prediction for the light intensity near rainbow caustics arises from the intrinsic defect of the ray model of light as have been discussed in Sec. 2.2.2, rather than the calculation method itself. Another notable discrepancy is caused by the fact that the GOA method does not include the surface wave effects [START_REF] Hovenac | Assessing the contributions of surface waves and complex rays to far-field Mie scattering by use of the Debye series[END_REF]. For example, the discrepancy in the backward scattering area where θ > 165 • is due to the surface wave of order p = 2. Fig. 2.10 analyzes the contributions of the surface waves of different orders, calculated with the Mieplot software [START_REF] Laven | MiePlot: A computer program for scattering of light from a sphere using Mie theory & the Debye series[END_REF], in the scattering by spherical particles of different radii.

One can see from Fig. 2.10(a) that for a spherical particle of radius a = 50 µm, the surface wave of order p = 2 has considerable intensity for θ > 165 • , so comes the notable discrepancy of GOA from LMT in this region as shown in Fig. 2.9(a). Besides, the surface waves of orders p = 1, 3 and 5 are also worthy of consideration in the scattering regions where 55 • < θ < 110 • , which leads to some perceivable discrepancies in these regions as shown the amplitudes of the scattered rays are increased with the size parameter ka. Thus, the effects of surface waves on a spherical particle of a larger size are weakened. As shown in Fig. 2.9, the scattered intensity for 55 • < θ < 110 • and θ > 165 • calculated by GOA agrees better with the LMT for a larger spherical particle. In the measurement of the size and the refractive index (temperature) of a spherical droplet/raindrop by the scattered light, the commonly-utilized information is the scattering pattern of the first-order rainbow. The zoomed view of Fig. 2.9(c) for the first-order rainbow is presented in Fig. 2.11. Except the discrepancy in the immediate vicinity of the rainbow angle (caustic), the agreement of GOA with LMT is very good in the other parts of the firstorder rainbow. And the rainbow fringes in those parts (θ > 138 • ) are fairly adequate in the inverse calculation for the size and refractive index. One distinct advantage of GOA over the rigorous LMT is the computational efficiency for large particles, which is critical for carrying out real-time measurements.

Summary

In this Chapter, the GOA for calculating the light scattering by infinite circular cylinders and spheres are discussed in detail. The calculation results are compared with the rigorous LMT. The comparisons show that the GOA method is valid when the size parameter of particle is over 100. And, the scattered intensity calculated by GOA agrees better with the LMT for a larger spherical particle.

However, since the light rays in the framework of GOA do not reflect the divergence or convergence characteristic of the represented wave, the GOA method encounters difficulties or even obstacles in accounting for the divergence factor and the phase shift due to focal lines when the geometry of particle/target becomes complex. In the following parts of this thesis, the VCRM, which resolves these difficulties encountered in GOA, will be introduced and extended.

Chapter 3 Light scattering by a cylinder of arbitrarily smooth surface based on VCRM

The study of the light scattering by an infinite cylinder provides an easier way to understand the interaction of light wave with particles. With regard to the practical uses, an efficient method to predict the light scattering by a cylinder is essential in the optical measurement techniques for the size and refractive index of optical fibers [START_REF] Adler | Rainbow scattering by a cylinder with a nearly elliptical cross section[END_REF][START_REF] Marcuse | Light scattering from unclad fibers: ray theory[END_REF][START_REF] Świrniak | Approximate solution for optical measurements of the diameter and refractive index of a small and transparent fiber[END_REF][START_REF] Świrniak | Numerical analysis of primary rainbows from a homogeneous cylinder and an optical fiber for incident low-coherent light[END_REF] and liquid cylinders [START_REF] Song | An optimization scheme for the measurement of liquid jet parameters with rainbow refractometry based on Debye theory[END_REF][START_REF] Song | Measurement of temperature gradient in a heated liquid cylinder using rainbow refractometry assisted with infrared thermometry[END_REF][START_REF] Han | Characterization of initial disturbances in a liquid jet by rainbow sizing[END_REF], and in the rendering of biological fibers such as hair and fur [START_REF] Bustard | Investigation into the scattering of light by human hair[END_REF][START_REF] Marschner | Light scattering from human hair fibers[END_REF][START_REF] Yan | An efficient and practical near and far field fur reflectance model[END_REF][START_REF] Khungurn | Azimuthal scattering from elliptical hair fibers[END_REF]. However, the available electromagnetic approaches for the light scattering by infinite cylinders are usually limited to those cylinders of circular or elliptical cross sections. To the author's best knowledge, there is still no satisfying solution for calculating the intensity distribution of the scattering field by a cylinder of other cross section.

In this chapter, based on the vectorial complex ray model (VCRM), a calculation method is proposed to calculate the light scattering, in the high-frequency limit, by an infinite cylinder of arbitrarily smooth cross section.

The difference of VCRM from GOA

For a particle of size much larger than the wavelength, i.e. in the high-frequency limit, the incident wave can be approximated by a bundle of rays, and each of them consists of propagation direction, polarization, amplitude and phase. The directions of the refracted and reflected rays can be calculated according to the Snell's law. The amplitudes are calculated by the Fresnel refraction and reflection coefficients. The phase lags when a light ray experiences a longer optical path; and it is also well known that the phase of reflected ray may have a shift compared to the phase of the incident ray. Besides the aforesaid ones, the amplitude and phase of a light ray are also affected by:

• The divergence/convergence of the represented wavefront. It is understandable by considering the fact that the light intensity of a convergent or divergent beam is varied at different axial positions, a result to meet the rule of energy conservation.

• The phase due to focal lines. A convergent beam has a phase shift after passing the focal lines, or the focal point (two crossing focal lines). This anomalous shift of phase is attributed to the confinement of wave either in spatial or in temporal, or even in spectral domain [1,[START_REF] Feng | Physical origin of the Gouy phase shift[END_REF][START_REF] Visser | The origin of the Gouy phase anomaly and its generalization to astigmatic wavefields[END_REF][START_REF] Boyd | Intuitive explanation of the phase anomaly of focused light beams[END_REF][START_REF] Andresen | Spectral analog of the Gouy phase shift[END_REF].

As have been discussed in the chapter 2, the divergence of wavefront (divergence factor) and the phase shift due to focal lines can be obtained in analytical forms based on GOA when a plane wave is refracted and reflected by a sphere or by a circular cylinder. But, more particles (or objects in a general sense) in nature are lack of the symmetry like a sphere or a circular cylinder. It then becomes rather tedious or even impossible for the conventional GOA method to compute the divergence of beam and the numbers of focal lines, which limits the precision of GOA in dealing with nonspherical particles.

In this context, Ren et al. [START_REF] Ren | Vectorial complex ray model and application to twodimensional scattering of plane wave by a spheroidal particle[END_REF] proposed the vectorial complex ray model (VCRM). Its advantages over the GOA method are that

• In VCRM, the curvature of wavefront is integrated as an intrinsic property of a light ray, by which one can calculate the divergence factor and the phase shift due to focal lines directly. It is because of this that the VCRM removes the biggest obstacle for GOA in calculating the light scattering by large nonspherical particles.

• The ray directions and the Fresnel coefficients are calculated by the wave vectors and their components, which avoids the tedious calculation in GOA when nonspherical particle is involved.

The VCRM is receiving more and more attention for its capability to solve the scattering of wave by large nonspherical particles with sufficient precision and high efficiency.

Calculation method

As to the 2D light scattering by an infinite cylinder, the available numerical implementation of VCRM addressed only the scattering by infinite cylinders of circular or elliptical cross sections [START_REF] Jiang | Scattering from an elliptical cylinder by using the vectorial complex ray model[END_REF]. In this section, a calculation method based on VCRM is proposed to solve the light scattering by an infinite cylinder of arbitrarily smooth cross section. The study is limited to the case when the cylinder is under normal incidence, i.e. all light rays stay in the transversal plane of the cylinder.

Consider a homogeneous cylinder (Fig. 3.1) of refractive index m and geometry given in the Cartesian coordinate system Oxyz by:

f (x, y) = 0 (3-1)
The z axis coincides with the main axis of the cylinder (perpendicular to the paper plane).

The cylinder has a smooth surface, in other words, the function f (x, y) is differentiable to the second order or more. The reason for this requirement will be explained later. The wave numbers in the surrounding medium and in the cylinder are k = 2π/λ and k m = mk, respectively, λ being the wavelength of light in free space.

Incident rays Cylinder

Ray tracing

In VCRM, the direction of a light ray is represented by the wave vector as shown in Fig. 3.2. Consider a light ray of wave vector ⃗ k i incident on the cylinder. The unit vectors normal to and tangent to the interface of cylinder are n and τ, respectively. The Snell's law indicates that the tangent components of the refracted wave vector ⃗ k t , the reflected ⃗ k r and the incident ⃗ k i are equal:

k (t) τ = k (r) τ = k (i) τ (3-2)
where k (i) τ = ⃗ k i • τ is the inner product of ⃗ k i with τ. Then, the normal components of the refracted wave and the reflected wave, k (t) n and k (r) n , are determined respectively by

k (t) n = ± √ k 2 t -k (t) τ 2 and k (r) n = -k (i) n (3-3)
where k t is the wave number of the refracted wave and the ± takes the sign of

k (i) n with k (i) n = ⃗ k i • n. The unit normal vector n is given by n = f ′ x ( f ′ x 2 + f ′ y 2 ) 1/2 • x + f ′ y ( f ′ x 2 + f ′ y 2 ) 1/2 • ŷ = n x • x + n y • ŷ (3-4)
where f ′ x and f ′ y are the derivatives of the surface function f (x, y) = 0 with respect to x and y, respectively. Besides, the unit tangent vector τ is given by τ

= -f ′ y ( f ′ x 2 + f ′ y 2 ) 1/2 • x + f ′ x ( f ′ x 2 + f ′ y 2 ) 1/2 • ŷ = τ x • x + τ y • ŷ (3-5)
After obtaining the normal and tangent components of the refracted and reflected rays, a simple transformation leads to the their expressions in the Cartesian coordinate system.

n Ŵ i k k r k k t k k Fig 3.
For example, the expression of the refracted wave vector in the Oxyz system is obtained by

⃗ k t = [ k (t) x k (t) y ] = [ k (t) τ τ x + k (t) n n x k (t) τ τ y + k (t) n n y ] (3-6)
where n x , n y and τ x , τ y are the x-and y-components of n and τ, respectively, which have already been given in Eqs. (3)(4) and (3)(4)(5). For a cylinder under normal incidence, the z-component of wave vector is zero.

To ease the calculation for the coordinates of the next interaction, a distance factor η is introduced to bridge the coordinates of two successive interaction points (x p , y p ) and

(x p+1 , y p+1 ):        x p+1 = x p + ηk x y p+1 = y p + ηk y (3-7)
where (k x , k y ) is the wave vector from the current interaction point (x p , y p ) to the next point (x p+1 , y p+1 ). By the fact that both (x p , y p ) and (x p+1 , y p+1 ) satisfy the same function f (x, y) = 0, we have

f (x p + ηk x , y p + ηk y ) -f (x p , y p ) = 0 (3-8)
This is an equation with η being the sole unknown. For an object of simple shape, this equation can be solved analytically. Otherwise, it is to be solved numerically by, for example, the Newton downhill method. It is worth noting that several solutions may be found for η if the cylinder has a complex geometry. One can prove this by drawing a line across a given closed curve, and counting the points of intersection as the geometry of the closed curve becomes complex. But, only the smallest positive value is the correct solution for η because a ray, unlike a line, has its starting point and positive direction. By substituting the solved η into Eqs. (3-7), the coordinates of the next interaction point can be obtained.

To designate the order of a scattered ray, the traditional notation (commonly used in the GOA and DSE) is adopted: the externally-reflected rays are of order p = 0, the transmitted rays without internal reflection are of p = 1, while the emergent rays which have been through p -1 times of internal reflections are the p-order scattered ray as illustrated in Fig. 3.3. The exit point and the scattering direction of a p-order ray is denoted as W p and ⃗ k p , respectively.

0 k k 0 W 1 W 2 W 1 p 2 p 1 k k 2 k k Incident ray 0 p Fig 3.
3 Notations for the order, the exit point and the direction of a scattered ray.

Amplitude

In the scattering by a dielectric cylinder, the amplitude of a light wave is affected by: 1) the Fresnel reflection/refraction coefficients and, 2) the divergence/convergence of wavefront. The Fresnel reflection and refraction coefficients in conventional geometric optics or ray optics are calculated with the incident and refracted angles [START_REF] Born | Principles of Optics[END_REF]. However, the acquisition of these angles is a rather daunting work when a nonspherical particle is involved, especially for the scattering in 3D space. For this reason, in VCRM the normal components of wave vectors are used instead of the incident and refracted angles to calculate the Fresnel coefficients:

r ⊥ = k (i) n -k (t) n k (i) n + k (t) n ; r ∥ = M 2 t k (i) n -k (t) n M 2 t k (i) n + k (t) n t ⊥ = 2k (i) n k (i) n + k (t) n ; t ∥ = 2M t k (i) n M 2 t k (i) n + k (t) n (3-9) k (i)
n and k (t) n are the normal components of the incident and the refracted wave vectors, respectively; and they should adopt their absolute values in Eq. (3)(4)(5)(6)(7)(8)(9) considering that the replaced incident and refracted angles range in [0, π/2]. M t is the refractive index of the refracted medium relative to the incident medium. For a refraction from air into the particle, M t = m; while for a refraction from the particle into air, M t = 1/m. The refraction and reflection coefficients at the interaction point W j ( j = 0, 1, 2, ...) are denoted as t X, j and r X, j , respectively. For a p-order scattered ray, the variation of amplitude caused by refraction and reflection coefficients is

ε X,p =            r X,0 , p = 0 t X,0 t X,p ∏ p-1 j=1 r X, j , p ≥ 1 (3-10)
Eq. (3-10) indicates that, for example, a scattered ray of order p = 2 has been through one refraction (t X,0 ) at W 0 , one internal reflection (r X,1 ) at W 1 and an eventual refraction (t X,2 ) into the surroundings at W 2 .

On the other hand, at the refraction or reflection by a particle/target of curved surface, the wave-front shape is altered. For example, an incident wave of planar wavefront is no longer a plane wave if it is refracted or reflected by a curved interface. The convergence or divergence of wavefront leads to the variation in light intensity, hence the amplitude.

For an infinite cylinder illuminated by a plane wave with normal incidence, the scattered rays have cylindrical wavefronts of varied curvatures. One principal direction of the wavefront is along the cylinder axis (z-axis) and the corresponding curvature is zero. The other one is in the xy plane; and at each interaction, the curvature radii 1 of the refracted and the incident wavefronts, R t and R i , are bridged by the curvature radius of cylinder ρ through [START_REF] Ren | Vectorial complex ray model and application to twodimensional scattering of plane wave by a spheroidal particle[END_REF][START_REF] Jiang | Scattering from an elliptical cylinder by using the vectorial complex ray model[END_REF][START_REF] Deschamps | Ray techniques in electromagnetics[END_REF]]

k (t) n 2 k t R t = k (i) n 2 k i R i + k (t) n -k (i) n ρ (3-11) Replacing k (t) n and k t respectively by k (r) n (k (r) n = -k (i) n
) and k i returns the curvature radius of the reflected wavefront R r :

k (i) n 2 k i R r = k (i) n 2 k i R i + -2k (i) n ρ (3-12)
For an infinite cylinder (or a 2D particle/target in a general sense) of cross section given by function f (x, y) = 0, the curvature radius ρ can be calculated through

ρ = ( f ′ x 2 + f ′ y 2 ) 3/2 f ′ y 2 f ′′ xx -2 f ′ x f ′ y f ′′ xy + f ′ x 2 f ′′ yy (3-13)
where f ′′ xx , f ′′ xy and f ′′ yy are the second-order derivatives. A careful reader may have noticed 1 Curvature radius is the reciprocal of curvature.

that the curvature radius ρ may be negative or positive. In fact, the sign of ρ determines whether the local surface of particle is convex or concave along the normal vector defined in Eq. (3-4):

• A positive curvature radius indicates the particle at this point has a convex surface along the normal vector;

• A negative curvature radius indicates a concave surface along the normal vector.

The calculation for the curvature radius ρ requires the second-order derivatives of the cylinder's surface function. Thus, the method based on VCRM for the scattering by cylinders is currently restricted to cylinders of smooth surface.

Besides the sudden change at refraction or reflection, a converging or diverging wavefront is subjected to gradual variation in homogeneous medium. For a wavefront of curvature radius R, after a propagation distance d, its curvature radius becomes

R ′ = R + d (3-14)
Here a positive curvature is attributed to a diverging wavefront. the incident ray (for a ray of central incidence, the impact factor is 0; while for the grazing ray, it equals 1). The incident ray is then refracted by the circular cylinder at W 0 and W 1 .

We can see that the wave-front shape experiences sudden changes precisely at the refraction.

Besides, the wavefront undergoes incremental change in the part from W 0 to W 1 and in the part from W 1 to the far field.

For the cylindrical wavefronts illustrated in Fig. 3.4, suppose the light intensity at the wavefront of curvature radius R is I, while the light intensity at the wavefront of curvature radius R ′ is I ′ . According to the conservation of energy, the energy fluxes passing through these two wavefronts are equal, indicating I ′ |R ′ | = I|R|. Thus, we have

I ′ = R R ′ I (3-15)
From Eqs. (3-11)-(3-13), the curvature radii of the refracted and the reflected wavefronts at each interaction with a cylinder of any smooth surface can be calculated. The curvature radius of the refracted wavefront at the interaction point W j ( j = 0, 1, 2, 3...) is denoted as R (t) j , while the curvature radius of the reflected wavefront at W j as R (r) j . For an externally reflected ray (p = 0) of curvature radius R (r) 0 , the light intensity at the reflection point is calculated as I X,0 r X,0 2 , where I X,0 is the intensity of the incident plane wave for polarization X (X =⊥ or ∥). After a distance of r from the exit point W 0 , its intensity can be calculated according to Eq. (3-15) as

I X,p=0,r = R (r) 0 R (r) 0 + r I X,0 r X,0 2 (3-16)
For a transmitted ray (p = 1), the incident wavefront have been through two refractions at W 0 and W 1 , at which the curvature radii of the refracted wavefronts are R (t) 0 and R (t) 1 , respectively. Besides the sudden changes at refractions, the wavefront is subjected to incremental change from W 0 to W 1 and from W 1 to the observation point. Thus, the light intensity of this p = 1 ray at the observation point is given by

I X,p=1,r = R (t) 0 R (t) 0 + s 0 R (t) 1 R (t) 1 + r I X,0 t 2 X,0 t 2 X,1 (3-17)
where s 0 is the geometric length from W 0 to W 1 , and r is the geometric length from the exit point W 1 to the observation point.

For a scattered ray of p = 2, the incident wavefront has been through a refraction at W 0 , an internal reflection at W 1 and a refraction at W 2 (W 2 is the exit point). The curvature radii of the refracted wavefront at W 0 is denoted as R (t) 0 , the reflected wavefront at W 1 as R (r) 1 and the refracted wavefront at W 2 as R (t) 2 . Besides these sudden changes, the wavefront is subjected to incremental change from W 0 to W 1 , from W 1 to W 2 and from W 2 to the observation point.

Thus, the light intensity of this p = 2 ray at the observation point is given by

I X,p=2,r = R (t) 0 R (t) 0 + s 0 R (r) 1 R (r) 1 + s 1 R (t) 2 R (t) 2 + r I X,0 t 2 X,0 r 2 X,1 t 2 X,2 (3-18)
where s 1 is the geometric length from W 1 to W 2 .

For a scattered ray of p ≥ 2, the incident wavefront has been through a refraction at W 0 , p -1 times of internal reflections respectively at W 1 , W 2 , ..., W p-1 and an eventual refraction at the exit point W p . By taking into account both the sudden changes at the reflections and refractions and the incremental changes along the ray paths, the light intensity of this p-order ray at the observation point is given by

I X,p,r = R (t) 0 R (t) 0 + s 0 R (r) 1 R (r) 1 + s 1 R (r) 2 R (r) 2 + s 2 ... R (r) p-1 R (r) p-1 + s p-1 R (t) p R (t) p + r I X,0 t 2 X,0 r 2 X,1 r 2 X,2 ...r 2 X,p-1 t 2 X,p (3-19)
where s j is the geometric length from W j to W j+1 ( j = 0, 1, ..., p -1), and r is the geometric length from the exit point W p to the observation point.

The whole effect of wave-front variation on the light intensity of a p-order scattered ray, calculated at a distance r from the exit point, is depicted by the divergence factor D p,r , with

D p,r =                          R (r) 0 R (r) 0 + r , p = 0 R (t) 0 R (t) 0 + s 0         p-1 ∏ j=1 R (r) j R (r) j + s j         R (t) p R (t) p + r , p ≥ 1 (3-20)
Thus, the p-order scattered intensity given in Eq. (3-19) can be expressed as

I X,p,r = D p,r I X,0 ε 2 X,p (3-21)
where ε X,p is the variation of amplitude due to refraction and reflection which has been discussed in Eq. (3)(4)(5)(6)(7)(8)(9)(10). The subscript r indicates the light intensity is calculated at a distance r from the exit point.

In the far-field case, i.e. r → ∞, the divergence factor for a p-order scattered ray is calculated as

D p =                    R (r) 0 , p = 0 R (t) 0 R (t) 0 + s 0         p-1 ∏ j=1 R (r) j R (r) j + s j         R (t) p , p ≥ 1 (3-22)
The dimensionless intensity of a light ray scattered by an infinite cylinder is then calculated in the far field as

I X,p = π 2 kD p ε 2 X,p (3-23)
A constant π/2 is multiplied so that the dimensionless scattered intensity calculated by GOA and VCRM is consistent with that by LMT [START_REF] Jiang | Theoretical study of light scattering by an elliptical cylinder[END_REF].

Then, the far-field amplitude for a p-order light ray of polarization X scattered by an infinite cylinder is obtained as

E X,p = ε X,p √ kD p π/2 (3-24)
When the refractive index of the cylinder has an imaginary part m i , an attenuation factor exp(-km i ∑ p-1 j=0 s j ) is to be multiplied, where ∑ p-1 j=0 s j is the total geometric path passed inside the cylinder.

Phase

When the cylinder is illuminated by light rays from a coherent light source, the phase of each scattered ray should be calculated correctly to account for the interference effect.

Here, the time factor is chosen as exp(iωt), which is consistent with that in [1].

(1) Phase differences of the incident rays When a shaped beam is concerned, the phases of the incident rays relative to a given plane (beam waist for example) are different. The calculation for the phase differences within the incident rays will be addressed in Chapter 6. In this chapter, the incident plane wave is under consideration, and the phase differences among the incident rays are zero.

(2) Phase shift due to optical path ϕ p,OP Those scattered rays emerging in the same direction will interfere with each other in the far field or on the focal plane of a focal lens. Since they may have undergone different optical paths, the corresponding phase shift ϕ p,OP differs between them. The optical path of a scattered ray is calculated by setting a reference ray. The reference ray arrives at the center of the cylinder in the direction of the incident ray and then emerges out in the same direction as the concerned scattered ray, as if there is no particle. Fig. 3.6 shows the schematic diagram in computing the optical path of a p = 1 ray relative to its reference ray. Defined in this way, the reference ray is consistent for all those scattered rays propagating in the same direction. In the scattering of plane wave by a spherical particle or a circular cylinder, the optical path can be calculated analytically as did in Chapter 2. For a cylinder of any cross section, an elegant formula is proposed here to calculate the phase shift due to optical path:

ϕ p,OP = 2π λ         σ i + σ p -m p-1 ∑ j=0 s j         (3-25)
where

σ i = ---→ W 0 O• ki and σ p = ---→ OW p
• kp as shown in Fig. 3.6. The ki and kp are the normalized wave vectors of the incident ray ⃗ k i at W 0 and the emergent ray ⃗ k p at W p , respectively. s j is the geometric length from the interaction point W j to W j+1 , and the sum ∑ p-1 j=0 s j represents the total geometric length experienced inside the cylinder.

(3) Phase shifts at reflection According to the Fresnel formulas shown in Eq. (3-9), r X may be negative, indicating the reflection may introduce a phase shift of π. Besides, when total reflection occurs, the normal component of the refracted wave vector in Eq. (3-3) becomes imaginary, i.e.

k (t) n = -i √ k τ 2 -k 2 t .
Consequently, r X is a complex number of unitary modulus and argument δ X

given by

               δ ⊥ = 2 arctan ( √ k i 2 -k (i) n 2 -k t 2 /k (i) n ) δ ∥ = 2 arctan ( √ k i 2 -k (i) n 2 -k t 2 /k (i) n /M 2 t ) (3-26)
where k i and k t are the wave numbers of the incident wave and the refracted wave, respectively. The explanation to k (i) n and M t can be found in Eq. (3)(4)(5)(6)(7)(8)(9).

The total phase shift of a p-order scattered ray due to reflection is then the argument (phase) of ε X,p .

ϕ X,r = arg(ε X,p ) (3-27)
(4) Phase shift due to focal lines A converging wave experiences a phase shift as it passes through its focus, a phenomenon known as the Gouy shift [START_REF] Gouy | Sur une propriété nouvelle des ondes lumineuses[END_REF][START_REF] Feng | Physical origin of the Gouy phase shift[END_REF][START_REF] Visser | The origin of the Gouy phase anomaly and its generalization to astigmatic wavefields[END_REF]. More explicitly, the phase advances by π/2 at the passage of each focal line or by π at a focal point (two crossing focal lines). In fact, the type A focal lines defined in [1] can be equivalently determined according to the convergence of wavefront in the meridional plane of sphere, while the type B focal lines according to the convergence of wavefront in the direction perpendicular to the meridional plane. In GOA, it is a daunting task to address the focal lines, especially when a nonspherical particle is involved. But in VCRM, the curvature of wavefront is integrated as an intrinsic property of a light ray, which allows us to determine directly whether and where a focal line is encountered.

The curvature radius of a light ray refracted or reflected by an infinite cylinder of arbitrarily smooth surface is calculated by Eq. (3)(4)(5)(6)(7)(8)(9)(10)(11). If the curvature radius is changed from negative to positive between two successive interaction points or between the exit point and infinity, it indicates that the wavefront has undergone a conversion from converging to diverging. Then, we can infer that a focal line has been encountered. Fig. 3.7 illustrates the focal line of a cylindrical wavefront, at which the wavefront is transformed from converging to diverging.

By counting the total number N of the focal lines that a light ray has encountered in its interaction with the cylinder, the corresponding phase shift is therefore

ϕ p,FL = Nπ/2 (3-28)
Finally, the complex amplitude for a p-order scattered ray of polarization X is obtained as ẼX,p = E X,p exp [ i

( ϕ X,r + ϕ p,OP + ϕ p,FL )] (3-29)
The diffraction component is calculated according to Eq. (2-24) with replacement of the radius of circular cylinder by the half width of the illuminated part.

Description of the scatterer

The method described in Sec. 3.2 is applicable to the scattering of plane wave by an infinite cylinder of arbitrarily smooth cross section. Here, it is applied to solving the light scattering by a composite elliptical cylinder (CEC) [START_REF] Adler | Rainbow scattering by a cylinder with a nearly elliptical cross section[END_REF][START_REF] Lock | Amplification of high-order rainbows of a cylinder with an elliptical cross section[END_REF]. In the Cartesian coordinate system as shown in Fig. 3.8, the upper part (y ≥ 0) of a CEC is defined by:

x 2 /a 2 + y 2 /b 2 1 = 1 (3-30)
while the lower part (y < 0) by

x 2 /a 2 + y 2 /b 2 2 = 1 (3-31)
where b 1 and b 2 are the semi-axes along y axis of the upper and the lower parts, respectively.

a is the common semi-axis in x direction. The direction of the incident wave makes an angle θ 0 with respect to the x axis, while the scattering angle θ is measured from the x axis. The angles are positive for counterclockwise rotations, while negative for clockwise rotations.

At (x, y), the curvature radius of the cylinder is given by

ρ = a 2 b 2 ( x 2 /a 4 + y 2 /b 4 ) 3/2 (3-32)
where b equals b 1 when y ≥ 0 or b 2 when y < 0.

The reason for choosing such a model is that, by changing the geometric parameters a, b 1 or b 2 , a CEC may take a variety of cross sections ranging from simple to highly-deformed.

A circular cylinder and an elliptical cylinder are among the special cases. Moreover, such a shape, when extended in the three-dimensional form, is a better model for a real raindrop/droplet than a sphere or a spheroid [5][6][7][8].

Suppose the wave vector originating from the interaction point W p is (k x , k y ), its next interaction with the cylinder occurs at W p+1 , the coordinates of which can be obtained by introducing the distance factor η, as have been given in Eq. (3-7). Two cases are involved according to the locations of the two successive interactions.

( (2). The two successive interactions occur on different parts: from upper to the lower part; or from lower to the upper part. In the former situation, η is given as

η = ( -B + √ B 2 -4AC ) / (2A) (3-34) with A = k 2 x a 2 + k 2 y b 2 2 B = 2 ( k x x p a 2 + k y y p b 2
2

) C = ( 1 b 2 2 - 1 b 2 1 ) y 2 p (3-35)
In the latter situation, the parameter η is calculated in the same way except exchanging b 1 and b 2 .

Then, the position of next interaction is determined by substituting the value of η into Eq. (3-7).

Results and discussion

Based on the method described above, a code has been written for calculating the scattered intensity of plane wave by CEC. The incident plane wave (λ = 0.6328 µm) is simulated by 4000 discrete parallel rays of equal phase and equal amplitude. These light rays propagate in the xy plane and make an incident angle θ 0 with the x axis. After the incidence, the light rays are subjected to continual refractions and reflections by the cylinder. Here, only the scattered rays of orders p = 0, 1, 2,...7 participate in the calculation considering that a higher-order ray has much lower intensity.

Since there is no other method available for calculating the scattered intensity of light by a cylinder of deformed transversal shape and transversal size much larger than the wavelength, the proposed method and the code are validated by comparing the results with existing theories for two special cases.

Firstly, the CEC is degenerated into a circular cylinder by setting b 1 = b 2 = a. The calculated scattering diagram is compared with the results of LMT [1,2], DSE [START_REF] Li | Debye series of normally incident plane-wave scattering by an infinite multilayered cylinder[END_REF] and GOA [1] as shown in Fig. 3.9. For VCRM, GOA and DSE, the scattered light of orders p ≤ 7 are taken into account. The relative refractive index m is 1.333. The incident angle θ 0 equals 0. The electric vector is polarized along the main axis of the cylinder (z axis), namely, the perpendicular polarization is concerned.

One can see that, when the CEC takes a circular cross section, the calculated result of VCRM is identical to that of the GOA. Moreover, the result of VCRM agrees well with those of the LMT and the DSE. The notable discrepancy is found in the immediate vicinity of caustics, for example, at the first-and the second-order rainbow angles. Another discrepancy, perceivable in the regions where 80 • < θ < 110 • and where θ → 180 • , is caused by surface waves. These discrepancies arise from the limitation of ray model [6,[START_REF] Hovenac | Assessing the contributions of surface waves and complex rays to far-field Mie scattering by use of the Debye series[END_REF][START_REF] Onofri | Experimental validation of the vectorial complex ray model on the inter-caustics scattering of oblate droplets[END_REF][START_REF] Laven | Supernumerary arcs of rainbows: Young's theory of interference[END_REF] rather than the calculation method itself. but it is worth considering the reasons for these differences. Unlike an infinite cylinder, the curvature of the elongated ellipsoid along z axis is never zero. This non-zero curvature causes convergence/divergence of wavefront and consequently introduces a shift in phase.

Besides, the formulas of divergence factor (Eq. (3)(4)(5)(6)(7)(8)(9)(10)(11)[START_REF] Li | Measurement of transient evaporation of an ethanol droplet stream with phase rainbow refractometry and high-speed microscopic shadowgraphy[END_REF][START_REF] Li | Simultaneous measurement of refractive index, diameter and colloid concentration of a droplet using rainbow refractometry[END_REF][START_REF] Langley | Critical-angle scattering of laser light from bubbles in water: measurements, models, and application to sizing of bubbles[END_REF][START_REF] Hansen | Mie scattering as a technique for the sizing of air bubbles[END_REF][START_REF] Arnott | Unfolded optical glory of spheroids: backscattering of laser light from freely rising spheroidal air bubbles in water[END_REF][START_REF] Zhang | Influence of bubbles on scattering of light in the ocean[END_REF][START_REF] Mougin | Path instability of a rising bubble[END_REF][START_REF] Lindken | A novel PIV technique for measurements in multiphase flows and its application to two-phase bubbly flows[END_REF][START_REF] Fujiwara | Bubble deformation and flow structure measured by double shadow images and PIV/LIF[END_REF][START_REF] Li | Influence of characteristics of micro-bubble clouds on backscatter lidar signal[END_REF][START_REF] Onofri | Near-criticalangle scattering for the characterization of clouds of bubbles: particular effects[END_REF]) and amplitude (Eq. (3-24)) for a cylinder are different from those for an ellipsoid (see Eqs. ( 2) and (3) in [START_REF] Ren | Scattering and transversal divergence of an ellipsoidal particle by using vectorial complex ray model[END_REF]), thus affecting the magnitude.

The scattering characteristics of light by CEC which has deformed cross section will be given and discussed in the following parts of this chapter. By a deformed cylinder, we mean a cylinder with cross section being neither circular nor elliptical.

Scattering patterns for diverse deformations

The cylinder is deformed such that its cross section approximates the shape of a natural raindrop [5]. The relation between such a deformed cylinder and a raindrop is much like that between a circular cylinder and a sphere. The size of a raindrop is described by the radius of the equivalent volume sphere r 0 , which is related to a, b 1 and b 2 through [8]

r 3 0 = a 2 ( b 1 2 + b 2 2 ) = a 3 ( ξ 1 2 + ξ 2 2 ) (3-36)
The axis ratios ξ 1 and ξ 2 are defined respectively by

ξ 1 = b 1 a , ξ 2 = b 2 a (3-37)
A natural raindrop of size larger than 0.4 mm cannot be treated as a spherical particle [6], and the shape deformation increases as the raindrop continues to grow. The two axis ratios ξ 1 and ξ 2 as functions of r 0 are given by [8]:

                   ξ 1 + ξ 2 = 2 ( 1 + f 1 • r ( f 2 + f 3 •r 0 + f 4 •r 3 0 ) 0 ) -1 ξ 1 -ξ 2 = f 5 • ( 2 ξ 1 + ξ 2 -1 ) • r 0 • exp ( -f 6 • √ r 0 ) (3-38)
with the coefficients f 1 = 0.08001, f 2 = 2.414, f 3 = -0.2911, f 4 = 0.009831, f 5 = 13.44 and f 6 = 2.508. The unit of the equivalent radius r 0 is millimeter. Fig. 3.11 shows the calculated ξ 1 and ξ 2 for different raindrop sizes. For a droplet of small size, the two axis ratios ξ 1 and ξ 2 tend to be 1, indicating the droplet approximates a spherical particle. As the size of raindrop becomes larger, the deformation from a spherical particle becomes more and more evident.

In this part, the effect of the transversal deformation of cylinder on its scattering field is investigated. Fig. 3.12 shows the calculated results for the intensity distribution of the scattering field by five CECs of diverse deformations. Here, the perpendicular polarization is concerned. The direction of the incident wave deviates from the x axis by -20 • , i.e. the incident angle θ 0 equals -20 • . The relative refractive index m is 1.333; and the wavelength of the incident wave λ is 0.6328 µm. These five CECs are deformed such that their cross sections approximate the shapes of a raindrop at different sizes. The equivalent radius of the raindrop r 0 ranges from 0.5 mm to 1.5 mm, while the corresponding geometric parameters a, b 1 and b 2 of the CECs are tabulated in Table 3.1. The axis ratios b 1 /a and b 2 /a are shown in the inset of Fig. 3.12. As r 0 ranges from 0.5 mm to 1.5 mm, the size parameter 2πr 0 /λ varies from approximately 5000 to 14900, far beyond the capabilities of the existing electromagnetic approaches. From Fig. 3.12, one can see that the scattering field has a distinct variation as the cylinder is deformed differently. One remarkable phenomenon is the shift of the first-order rainbow in the downward scattering region (θ < 0 • ), whose scattering angle θ is changed by 16.4 • as the equivalent radius r 0 is increased from 0.5 mm to 1.25 mm. The sensitivity indicates that the first-order rainbows from two CECs of diverse deformations have distinctly-different scattering angles. This adds to the evidence that the twinned rainbows in nature arise from the scattering by two groups of raindrops which have different sizes and at least one of them is nonspherical [6]. Besides, it is found that even a single deformed CEC could bring forth twinned second-order rainbows as shown in Fig. 3.12 when r 0 = 1.50 mm, in which case the cross section of CEC has a considerable deformation (b 1 /a = 0.931, b 2 /a = 0.762).

Scattering patterns for varied refractive indices

Besides the deformation of the transversal shape, the refractive index of CEC is another factor worthy of consideration. Here, the medium of CEC is assumed to be pure water. The refractive indices of pure water for λ = 0.6328 µm at different temperatures are tabulated [START_REF] Harvey | Revised formulation for the refractive index of water and steam as a function of wavelength, temperature and density[END_REF] in Table 3.2. The intensity distribution of the scattering field by a CEC at these temperatures are presented in Fig. 3.13. The perpendicular polarization is concerned here. The cross section of the CEC takes the shape of a raindrop when r 0 = 0.75 mm (a = 0.760 mm, b 1 /a = 0.984 and b 2 /a = 0.935). The incident angle θ 0 is equal to -20 • (same with that in Fig. 3.12). From Fig. 3.13, one can see that as the temperature of water is increased from 5 • C to 30 • C, the scattering field is subjected only to a slight change. This is due to the fact that the refractive index m is altered only by 0.0019 as given in Table 3.2. Without doubt, even a small variation in m has an impact on the scattering pattern because the refractive index

Scattering patterns for different incident angles

Due to the asymmetric cross section of a deformed CEC, the scattered light is distributed differently if the direction of incident wave is changed. Thus, the scattering patterns for different incident angles are also needed for a full understanding of the scattering characteristics of light by a cylinder of deformed cross section.

For a CEC of fixed cross section and refractive index, the variation of its scattering field with the direction of incident wave is shown in Fig. 3.15. The CEC used here is of the same cross section as the one in Fig. 3.13, while the refractive index m is fixed at 1.333. The incident angle θ 0 ranges from -20 • to 20 • in steps of 10 • , which simulate the cases when the sun has an inclination of ±20 • relative to the horizontal plane of the concerned raindrops.

Since the scattering angle θ, defined in Fig. 3.8, is measured from the preset x axis, the scattering diagram has a holistic movement with the incident angle θ 0 as shown in Fig. 3.15.

But, unlike a circular cylinder, the scattering field by a CEC of deformed cross section is altered, if we go into detail, for different incident angles. For example, in the downward scattering regions (θ < 0), the fifth-order rainbow has rather low intensity when θ 0 = -20 • , while it is submerged by the second-order rainbow when θ 0 = -10 • , but it becomes relatively clear when θ 0 = 10 • and 20 • . One can see that the angular width between the first-order and the second-order rainbows (a region known as the Alexander's dark band) for the CEC is much wider than that for the circular cylinder. It indicates that a natural raindrop has a wider Alexander's dark band than a spherical particle of equivalent radius. This is in accordance with the prediction made by Sadeghi et al. [6]. Besides, as shown in Fig. 3.16(b), the fifth-order rainbow formed by the CEC might separate from the second-order rainbow at certain incident angles, which makes it possible to observe the fifth-order rainbow. But for the circular cylinder, the fifthorder rainbow is not within the Alexander's dark band and is overwhelmed by the secondorder rainbow, whose intensity is much greater. This gives insight into the conditions under which one might be able to observe a natural fifth-order rainbow [START_REF] Edens | Photographic observation of a natural fifth-order rainbow[END_REF].

In Sec. 3.4, the scattering patterns of the CECs whose cross sections take the shapes of natural raindrops are calculated in the whole scattering region. The effects of shape deformation, refractive index and the direction of incident wave on the scattering characteristics, especially on the rainbows, are quantitatively analyzed. It has been found that

• The first-order rainbows formed by CECs of diverse deformations are distinctly different in scattering angle, which may explain why the formation of the twinned rainbows in nature requires two groups of raindrops of different shapes.

• The variation of atmospheric temperature (5 • C ∼ 30 • C) only has a slight effect on the rainbow position.

• The Alexander's dark band of the CEC whose cross section is deformed like a natural raindrop is wider than that of a circular cylinder.

• The elevation angle of the sun affects the emergence or submergence of the fifth-order rainbow.

Although these phenomena are found in the scattering patterns of deformed cylinders, they are capable of explaining certain natural phenomena produced by deformed raindrops. This is because the relation between a deformed CEC and a nonspherical raindrop in nature is much like the relation between a circular cylinder and a spherical particle.

Summary

This chapter reported the extension of VCRM allowing to account, in the high-frequency limit, for the direction, polarization, curvature of wavefront, amplitude, phase and scattered intensity of the light rays interacting with an infinite cylinder of arbitrary while smooth cross section. Based on the proposed method, a numerical study has been performed on the scattering patterns of composite elliptical cylinders (CECs), whose cross sections can take various shapes ranging from circular, elliptical to highly-deformed. The effects of shape deformation, refractive index and the direction of incident wave on the scattering patterns of the CECs whose cross sections approximate the shapes of natural raindrops were investigated and quantitatively analyzed, which provided insight into how these factors affect the appearance of a natural rainbow.

Being flexible and numerically efficient (a full scattering diagram is obtained in few seconds on a laptop computer (Intel i7-8550U @1.80GHz)), the proposed method is thought to have important applications in calculating and analyzing the scattering characteristics of light by cylindrical objects of cross sections ranging from simple to complex.

Chapter 4 Calculation method for the 3D light scattering by a large nonspherical particle

In the preceding chapter, the light scattering by an infinite cylinder of arbitrarily smooth cross section under normal incidence was solved in the framework of VCRM, where light rays stayed in a single plane and the scattering was just a two-dimensional (2D) problem. In this chapter, an algorithm based on VCRM is proposed to calculate the three-dimensional (3D) scattered intensity by a large nonspherical particle of any smooth surface.

Overview of the difficulties in solving 3D scattered intensity

Once the scattered rays are spread into the 3D space, several difficulties are encountered in the calculation for the scattered intensity.

1. The first difficulty is the tedious computation in ray tracing. The conventional GOA method which tries to compute the incident and refracted angles may already be frustrating here. But with VCRM, where the directions of light rays are described by their wave vectors, one can carry out the ray tracing in 3D space much more easily.

2. The second one is the calculation for the local curvatures of the concerned nonspherical particle at each interaction point; and the following calculation for the principal curvatures and principal directions of the wavefront after refraction or reflection.

3. The third difficulty is to deal with cross polarization effect, considering that the polarization state is usually varied at different points of interaction.

4. The fourth one is to compute the phase shifts when a light ray propagates in 3D space.

5. The last and the bottle-neck one is to account for the interference effect, since the light rays may be scattered into any direction and may be even folded in the 3D space.

These problems encountered in the 3D light scattering by large nonspherical particles will be addressed one by one in the following parts of this chapter.

Ray tracing in 3D scattering

In Sec. 3.2.1, the tracing for the light rays propagating within a single plane was described, where the plane of incidence Σ inc , consisting of the wave vectors and the normal vector, was consistent at all points of interaction. As soon as a 3D scattering problem is involved, the light rays are no longer in the same plane and the plane of incidence Σ inc might be altered for two successive points of interaction. We need therefore to determine the Σ inc firstly, so that the vectorial Snell's law could be applied. At one interaction, the wave vector of the incident ray is denoted as ⃗ k i and the unit normal vector of the dioptric surface is calculated as n (the formula for calculating the normal vector will be given in next section). The Σ inc at this point is the plane containing the ⃗ k i and n as shown in Fig. 4.1. The unit vector perpendicular to Σ inc is defined by

ê⊥ = ⃗ k i × n ⃗ k i × n (4-1)
Then, the unit vector τ tangent to the dioptric surface and located in Σ inc is determined by:

τ = n × ê⊥ (4-2)
The normal component and the tangent component of the incident wave vector, k n and k τ , are then obtained by the inner products ⃗ k i • n and ⃗ k i • τ, respectively. According to Eqs. (3-3) and (3-2), the normal component and the tangent component of the refracted or reflected wave vector can be acquired.

In numerical implementation, the expression of a wave vector from the local coordinate system (τ, n, ê⊥ ) to the Cartesian coordinate system can be transformed easily. For example, the refracted wave vector is given in ( x, ŷ, ẑ) by

⃗ k t =           k (t) x k (t) y k (t) z           =           k (t) τ τ x + k (t) n n x k (t) τ τ y + k (t) n n y k (t) τ τ z + k (t) n n z           (4-3)
where n x , n y and n z are the three components of the unit normal vector n along x, ŷ and ẑ, respectively; while τ x , τ y and τ z the three components of the unit tangent vector τ.

A light ray propagates straightforward in a homogeneous medium. The coordinates of two successive interaction points (x p , y p , z p ) and (x p+1 , y p+1 , z p+1 ) are related by

          x p+1 y p+1 z p+1           =           x p y p z p           + η •           k x k y k z           (4-4)
where (k x , k y , k z ) is the wave vector between the two interactions. η is called the distance factor (η > 0).

Suppose that the surface function of particle is f (x, y, z) = 0. The distance factor η between two successive points can be solved by equation

f (x p + ηk x , y p + ηk y , z p + ηk z ) -f (x p , y p , z p ) = 0 (4-5)
Thus, the search for the position of the next interaction is now boiled down to solving this one-variable equation. It should be noted that there may be several solutions for η, and the minimum positive one corresponds to the actual position of the next interaction. The reason for this has been discussed in Sec. 3.2.1. Then, the coordinates of the next interaction point are obtained by substituting the value of η into Eq. (4-4).

Through the procedure described above, the coordinates of each interaction point and the direction of each scattered ray can be calculated step by step. The following parts will address the amplitude and the phase. In the free propagation within a homogeneous medium, the curvature radii of a converging or diverging wavefront are changed gradually with the propagation distance. After a distance d, the curvature radii can be directly calculated according to

Curvature of wavefront

R ′ j = R j + d (4-7)
where j = 1 or 2. A positive curvature is attributed to a diverging wavefront, as illustrated in Fig. 3.4. (û 2 , ⃗ k), respectively. If F 1 F 2 , the wavefront is astigmatic and will have a saddle surface between F 1 and F 2 , where R ′ 1 < 0 but R ′ 2 > 0. After passing F 2 and F 1 , at the point O ′ for example, the local wavefront becomes diverging with R ′ 1 > 0 and R ′ 2 > 0. As to the two principal directions û1 and û2 , they have no change in a homogeneous medium.

2 û k k 1 û 2 û 1 û O 1 R 2 R 2 ' R ' O 2 F 1 ' R 1 F converging saddle diverging
Besides the gradual variation in the free propagation within a homogeneous medium, the wavefront is subjected to a sudden change upon being refracted or reflected. For example, when a planar wavefront is refracted or reflected by a curved interface, it may be converged or diverged dependent on the two principal curvatures of the particle surface and the relative refractive index. To calculate the variation of wavefront at a refraction or at a reflection, the shape information of the particle, more precisely, the principal curvatures and the principal directions of the particle's local surface at the refraction or reflection are required. The computation for the principal curvatures and the principal directions of a 3D surface involves the knowledge of Differential Geometry. This chapter presents the calculation method for 3D particles/objects of parameterized surface X = X(α, β). For those 3D particles/objects of implicit surface f (x, y, z) = 0, the calculation method will be presented in the following Chapter 6, where the iso-phase surface of elliptical Gaussian beam has an implicit function.

Let X = X(α, β) be the parametric expression of particle surface, where X is a vectorvalued function of the parameters α and β1 . More knowledge about parameterized surface can be found in Differential Geometry books and literatures, such as [START_REF] Pressley | Elementary Differential Geometry[END_REF][START_REF] Koenderink | Differential Geometry of Surfaces in Three-Dimensional Euclidean Space[END_REF][START_REF] Tapp | Differential Geometry of Curves and Surfaces[END_REF]. The first-order derivatives of X are calculated as X ′ α = ∂X/∂α and X ′ β = ∂X/∂β, while the secondorder derivatives are X ′′ αα = ∂ 2 X/∂α 2 , X ′′ αβ = ∂ 2 X/∂α∂β and X ′′ ββ = ∂ 2 X/∂β 2 . The X ′ α and X ′ β are two vectors tangent to the particle surface, indicating that the unit vector normal to the particle surface can be obtained by

n = X ′ α × X ′ β ∥X ′ α × X ′ β ∥ (4-8)
The coefficients of the first fundamental form [START_REF] Pressley | Elementary Differential Geometry[END_REF][START_REF] Koenderink | Differential Geometry of Surfaces in Three-Dimensional Euclidean Space[END_REF][START_REF] Tapp | Differential Geometry of Curves and Surfaces[END_REF] (E, F, G) are calculated by

             E = X ′ α • X ′ α F = X ′ α • X ′ β G = X ′ β • X ′ β (4-9)
while the coefficients of the second fundamental form (L, M, N) by

             L = X ′′ αα • n M = X ′′ αβ • n N = X ′′ ββ • n (4-10)
Then, the two principal curvatures and the corresponding principal directions of the particle surface can be found in the Weingarten equations in matrix form [START_REF] Gallier | Geometric Methods and Applications[END_REF]

J weingarten = 1 EG -F 2 [ MF -LG NF -MG LF -ME MF -NE ] (4-11)
The eigenvalues of J weingarten are -1/ρ 1 and -1/ρ 2 , where 1/ρ 1 and 1/ρ 2 are the two principal curvatures of the particle surface (ρ 1 and ρ 2 are the principal curvature radii). The eigenvectors of J weingarten , ⃗ v 1 and ⃗ v 2 , correspond to the principal directions of the particle surface. The curvature matrix of the particle surface is then obtained as

C = [ 1/ρ 1 0 0 1/ρ 2 ] (4-12)
A careful reader may have noticed that this calculation method for the principal curvatures and principal directions of particle surface requires the second-order derivatives of particle surface. It is because of this requirement that the current VCRM is restricted to particles of smooth surface.

It should be noted that the two principal directions ⃗ v 1 = [v (1) α , v (1) β ] and ⃗ v 2 = [v (2) α , v (2) β ] are two orthogonal vectors in the space expanded by X ′ α and X ′ β . They can be transformed into the Cartesian coordinate system ( x, ŷ, ẑ) through:

       ⃗ v 1 = v (1) α X ′ α + v (1) β X ′ β ⃗ v 2 = v (2) α X ′ α + v (2) β X ′ β (4-13)
Their normalized vectors are denoted as v1 and v2 , respectively.

From the calculation above, the curvature matrix C (whose diagonal elements are the two principal curvatures 1/ρ 1 and 1/ρ 2 ) and the two principal directions v1 and v2 at a local point of particle surface are obtained.

On the other hand, suppose that the two principal curvature radii of the incident wavefront are respectively R (i) 1 and R (i) 2 , and the corresponding principal directions are û(i) 1 and û(i) 2 . The curvature matrix of the incident wavefront is denoted as

Q i = [ 1/R (i) 1 0 0 1/R (i) 2 ] (4-14)
The projection matrix between (û (i) 1 , û(i) 2 ) and (v 1 , v2 ) is defined by [START_REF] Deschamps | Ray techniques in electromagnetics[END_REF]

Θ i = [ û(i) 1 • v1 û(i) 1 • v2 û(i) 2 • v1 û(i) 2 • v2 ] (4-15)
To describe the wavefront of the refracted ray, an orthogonal coordinate system (û (t) 1 , û(t) 2 ) in the tangent plane of the refracted wavefront is needed. For simplicity, û(t)

1 is simply set as the ê⊥ which has already been defined in Eq. (4-1), while û(t)

2 is defined by

û(t) 2 = ê⊥ × ⃗ k t ê⊥ × ⃗ k t (4-16)
The projection matrix between (û (t) 1 , û(t) 2 ) and (v 1 , v2 ) is

Θ t = [ û(t) 1 • v1 û(t) 1 • v2 û(t) 2 • v1 û(t) 2 • v2 ] (4-17)
Then, the curvature matrix of the refracted wavefront Q t is related to that of the incident wavefront Q i by [START_REF] Ren | Vectorial complex ray model and application to twodimensional scattering of plane wave by a spheroidal particle[END_REF][START_REF] Deschamps | Ray techniques in electromagnetics[END_REF][START_REF] Ren | Vectorial Complex Ray Model for light scattering of nonspherical particles[END_REF]:

k t Θ T t Q t Θ t = k i Θ T i Q i Θ i + ( k (t) n -k (i) n ) C (4-18)
where k i and k t are the wave numbers of incident and refracted waves, respectively. k (i) n and k (t) n are the normal components of the incident and refracted wave vectors, respectively. The letter T indicates the transposition of matrix.

The curvature matrix Q t contains the information about the shape of the wavefront after the refraction. The eigenvalues of Q t are the two principal curvatures 1/R (t) 1 and 1/R (t) 2 of the refracted wavefront, while the eigenvectors are the corresponding principal directions. It should be noted that the eigenvectors of Q t are two vectors in the local coordinate system (û (t) 1 , û(t)

2 ), and they can be transformed into the ( x, ŷ, ẑ) system in the same way as did in Eq. (4-13). While for the reflected wavefront, its principal curvatures and principal directions are obtained in the similar way.

After the refraction or reflection, the curvature radii of wavefront are again subjected to incremental variation as given in Eq. (4-7) when the light ray propagates in a homogenous medium.

Through the procedure given above, the principal curvatures and principal directions of the refracted and the reflected wavefronts at each interaction point with a particle of any smooth surface can be calculated step by step. At the interaction point W j ( j = 0, 1, 2, 3...), the two principal curvature radii of the refracted wavefront are denoted as R (t) 1, j and R (t) 2, j , while the principal curvature radii of the reflected wavefront as R (r) 1, j and R (r) 2, j . Different from the divergence factor of a cylindrical wavefront as deduced in Sec. 3.2.2, the calculation for the divergence factor of a general wavefront in 3D scattering involves two principal curvature radii. For a p-order scattered ray at a distance r from the exit point, the divergence factor which accounts for the overall variation in light intensity due to the divergence or convergence of wavefront is given by

D p,r =                          R (r) 1,0 R (r) 2,0 (R (r) 1,0 + r)(R (r) 2,0 + r) , p = 0 R (t) 1,0 R (t) 2,0 (R (t) 1,0 + s 0 )(R (t) 2,0 + s 0 )         p-1 ∏ j=1 R (r) 1, j R (r) 2, j (R (r) 1, j + s j )(R (r) 2, j + s j )         R (t) 1,p R (t) 2,p (R (t) 1,p + r)(R (t) 2,p + r) , p ≥ 1 (4-19)
where s j is the geometric length from W j to W j+1 ( j = 0, 1, ..., p -1), and r is the geometric length from the exit point W p to the observation point.

In the far-field case, i.e. r → ∞, the divergence factor for a p-order scattered ray is calculated as

D p =                    R (r) 1,0 R (r) 2,0 , p = 0 R (t) 1,0 R (t) 2,0 (R (t) 1,0 + s 0 )(R (t) 2,0 + s 0 )         p-1 ∏ j=1 R (r) 1, j R (r) 2, j (R (r) 1, j + s j )(R (r) 2, j + s j )         R (t) 1,p R (t) 2,p , p ≥ 1 (4-20)
The divergence factor D p calculates the variation in light intensity caused by the diver-gence or convergence of wavefront for a p-order scattered ray, thus the value √ D p is the corresponding variation in amplitude. In the following part, we will discuss another factor affecting the amplitude.

Cross polarization

In the special cases when the incident plane Σ inc at each interaction remains unchanged, the polarization state does not change if the electric vector is either perpendicular or parallel to Σ inc . The scattering by a sphere [1,[START_REF] Xu | Extension of geometrical-optics approximation to on-axis Gaussian beam scattering. I. By a spherical particle[END_REF], by an infinite cylinder under normal incidence [START_REF] Adler | Rainbow scattering by a cylinder with a nearly elliptical cross section[END_REF][START_REF] Lock | Amplification of high-order rainbows of a cylinder with an elliptical cross section[END_REF][START_REF] Jiang | Scattering from an elliptical cylinder by using the vectorial complex ray model[END_REF][START_REF] Jiang | Scattering of a Gaussian beam by an elliptical cylinder using the vectorial complex ray model[END_REF][START_REF] Marcuse | Light scattering from unclad fibers: ray theory[END_REF][START_REF] Świrniak | Approximate solution for optical measurements of the diameter and refractive index of a small and transparent fiber[END_REF], or by a spheroid with end-on incidence [START_REF] Xu | Extension of geometrical-optics approximation to on-axis Gaussian beam scattering. II. By a spheroidal particle with end-on incidence[END_REF] is among these special cases. These configurations where polarization state is preserved have one basic requirement: the incident plane Σ inc is consistent for different points of interaction so that the electric vectors are always perpendicular or parallel to the Σ inc .

However, in general cases, the incident planes at successive interactions might not coincide and consequently the cross polarization occurs [START_REF] Lock | Ray scattering by an arbitrarily oriented spheroid. II. Transmission and crosspolarization effects[END_REF][START_REF] Adler | High-order interior caustics produced in scattering of a diagonally incident plane wave by a circular cylinder[END_REF]. For an electric vector of vibration direction perpendicular (or parallel) to the Σ inc at one interaction, the vibration direction at the next interaction where Σ inc is different will be no longer perpendicular (or parallel) to the incident plane. In such cases, a linearly-polarized electric vector will change its direction of vibration after the refraction or reflection because the Fresnel coefficients for the perpendicular component and the parallel component are usually different. Furthermore, if total reflection occurs, after the reflection, the linear polarization becomes elliptical polarization since the phase difference between the perpendicular and the parallel components may reach π/6 when the relative refractive index is 1/1.3322.

This part covers the solution for the polarization state, namely the vibration direction and the amplitude of the electric vector, when a light ray is scattered by a 3D nonspherical particle.

Suppose the Jones vector of the incident electric vector at one interaction is given by

⃗ E i = [ E 1 e iϕ 1 E 2 e iϕ 2 ] (4-21)
where the complex amplitudes E 1 e iϕ 1 and E 2 e iϕ 2 are the components of the incident electric vector along ê1 and ê2 , respectively, with ê1 and ê2 being two orthogonal base vectors within the vibration plane of the incident elector vector. The phases ϕ 1 and ϕ 2 discussed in this section only refer to the phase shifts caused by reflection, since those phase shifts due to optical path and focal lines are independent of polarization.

To apply the Fresnel formulas, one needs the parallel and perpendicular components of ⃗ E i with respect to the incident plane Σ inc at this interaction, i.e. the component along the ê(i) ∥ which is parallel to the Σ inc , and the component along the ê⊥ which is perpendicular to Σ inc .

The direction ê⊥ has been defined in Eq. (4-1), while ê(i)

∥ can be determined by ∥ , ê⊥ ) within the vibration plane of incident electric vector.

ê(i) ∥ = ê⊥ × ⃗ k i ê⊥ × ⃗ k i
In the vibration plane of the incident electric vector as illustrated in Fig. 4.5, the coordinate transformation matrix from (ê 1 , ê2 ) to (ê (i) ∥ , ê⊥ ) is calculated by2

M 1 = [ (ê (i) ∥ • ê1 ) (ê (i) ∥ • ê2 ) (ê ⊥ • ê1 ) (ê ⊥ • ê2 ) ] (4-23)
Then, the Jones vector of the incident electric vector in the coordinate system (ê (i) ∥ , ê⊥ ) is obtained as

       E (i) ∥ e iϕ (i) ∥ E (i) ⊥ e iϕ (i) ⊥        = M 1 ⃗ E i (4-24)
One noteworthy fact about Eq. (4-24) is that the elements of Jones vector depend on the choice of the coordinate system. If these elements (complex amplitudes) are known in one coordinate system, they can be determined in another coordinate system by using coordinate transformation matrix [START_REF] Saleh | Fundamentals of Photonics[END_REF].

According to Fresnel formulas, the Jones vector of the electric vector after refraction can be obtained by

⃗ E t =        E (t) ∥ e iϕ (t) ∥ E (t) ⊥ e iϕ (t) ⊥        = [ t ∥ 0 0 t ⊥ ]        E (i) ∥ e iϕ (i) ∥ E (i) ⊥ e iϕ (i) ⊥        = J t M 1 ⃗ E i (4-25)
where J t is the Jones matrix for refraction, and its diagonal elements are the Fresnel refraction coefficients. On the other hand, the Jones vector of the electric vector after reflection is calculated by

⃗ E r =        E (r) ∥ e iϕ (r) ∥ E (r) ⊥ e iϕ (r) ⊥        = [ r ∥ 0 0 r ⊥ ]        E (i) ∥ e iϕ (i) ∥ E (i) ⊥ e iϕ (i) ⊥        = J r M 1 ⃗ E i (4-26)
where J r is the Jones matrix containing the Fresnel coefficients for reflection. The Fresnel reflection and refraction coefficients r ∥ , r ⊥ and t ∥ , t ⊥ are calculated directly by the normal components of the incident and refracted wave vectors as given in Eq. (3-9). ∥ , ê(r) ⊥ ), respectively, with ê(t) ⊥ = ê(r) ⊥ = ê⊥ , while ê(t) ∥ and ê(r) ∥ are defined in the same way as did in Eq. (4-22) just by replacing ⃗ k i with ⃗ k t and with ⃗ k r , respectively. Fig. 4.6 presents an illustration of these base vectors. The wave vectors ⃗ k i , ⃗ k r and ⃗ k t lie in the same plane, the incident plane Σ inc at this interaction. The base vector ê⊥ is perpendicular to the Σ inc , while ê(i)

∥ , ê(r) ∥ and ê(t) ∥ are within the Σ inc . As to the vibration directions of the electric vectors, for example, the refracted electric vector vibrates within the plane constituted by ê⊥ and ê(t)

∥ . The refracted or the reflected light ray may then act as the incident ray at the next interaction. If so, the procedure from Eq. (4-21) to Eq. (4-26) is to be repeated. Otherwise, it emerges out of the particle and contributes to the scattering field. For a scattered ray of wave vector ⃗ k s (s = r when the reflected ray is scattered out of the particle, while s = t if the refracted ray is scattered out), two orthogonal unit vectors in the vibration plane of ⃗ E s (electric vector of the scattered ray ⃗ k s ) are introduced as

êφ = ẑ × ⃗ k s ẑ × ⃗ k s and êθ = ⃗ k s × êφ ⃗ k s × êφ . (4-27)
For a fixed Cartesian coordinate system ( x, ŷ, ẑ), the two base vectors êθ and êφ are only related to the direction of the scattered ray ⃗ k s /∥ ⃗ k s ∥. It indicates the coordinate system (ê θ , êφ ) for describing the electric vector is consistent for those scattered rays emerging in the same direction.

In the vibration plane of ⃗ E s , the coordinate transformation matrix from the coordinate system (ê (s) ∥ , ê⊥ ) to (ê θ , êφ ) is calculated as

M 2 =       (ê θ • ê(s) ∥ ) (ê θ • ê⊥ ) (ê φ • ê(s) ∥ ) (ê φ • ê⊥ )       (4-28)
Then, the Jones vector of ⃗ E s expressed in (ê θ , êφ ) is obtained:

⃗ E s =       E θ e iϕ θ E φ e iϕ φ       = M 2 J s M 1 ⃗ E i (4-29)
The complex amplitudes E θ e iϕ θ and E φ e iϕ φ are the two orthogonal components of the electric vector of scattered ray along êθ and êφ , respectively. For any two scattered rays propagating in the same direction, their electric vectors are expressed in the same coordinate system.

Thus, the calculation method described above enables the straightforward addition of two scattered electric vectors for superposition. Besides, the calculation for the incident and refracted angles (as did in [START_REF] Adler | High-order interior caustics produced in scattering of a diagonally incident plane wave by a circular cylinder[END_REF]), which could be rather tedious especially in 3D scattering, is avoided. The method proposed here simplifies considerably the calculation for the polarization states of the light rays in 3D scattering when cross polarization is involved.

Phase shift in 3D scattering

In the preceding Sec. 3.2.3, the factors which affect the phase shift in the 2D scattering of plane wave by a cylinder of any smooth cross section has been discussed in detail. For the 3D scattering of plane wave by a nonspherical particle of any smooth surface, some of those factors need further discussions.

As to the phase shift due to optical path in the 3D scattering of plane wave by a particle of any shape, the formula proposed in Eq. to implement. It only needs the direction of the incident ray ki at the incident point W 0 , the direction of the scattered ray kp , and the coordinates at each interaction W 0 , W 1 , ...W p . More importantly, the proposed formula is applicable to a 2D or 3D particle of any shape.

Due to the presence of cross polarization, the electric vector of a 3D scattered ray is resolved into two components, E θ e iϕ θ and E φ e iϕ φ as given in Eq. . The arguments ϕ θ and ϕ φ record the phase shifts due to reflection for the two polarizations.

In the 3D scattering by a nonspherical particle, the wavefront of a scattered ray in the near field 3 is usually astigmatic or spherical, presenting as a curved surface in the 3D space.

Compared to a cylindrical wavefront, an astigmatic or spherical wavefront has two nonzero principal curvatures. The calculation for the two principal curvatures of a scattered wavefront can be found in Sec. 4.3. If any of the two curvatures is changed from negative to positive between two successive interaction points or between the exit point and infinity, it indicates that the wavefront has undergone a conversion from converging to diverging in the corresponding principal direction. We can then infer that a focal line has been encountered in this principal direction. By counting the total number N of the focal lines that have been encountered in both of the principal directions, the phase shift due to focal lines is then calculated by Nπ/2.

The phase difference between the incident rays is zero for plane wave incidence. The case with shaped beam incidence will be discussed in detail in the Chapter 6.

Interpolation of amplitude and phase in 3D scattering

When a light wave is scattered by a nonspherical particle, the scattered rays are usually irregularly distributed. For a concerned scattering direction, the calculation for the scattered intensity usually involves the interpolation of the scattered rays nearby. However, the interpolation in 3D scattering is not an easy task since:

1. Both amplitude and phase are bivariate functions.

2. The light rays from the same scattering order or from different orders may be folded one or more times, indicating that the amplitude and the phase may also be multi-valued functions. The superposition of light in the folding area must be taken into account.

Sadeghi et al. [6] used quadrilateral patches, each of which was composed of four rays, to evaluate the amplitude and phase shift due to the convergence/divergence of wavefront; then they implemented bilinear interpolation in each patch to calculate the amplitude and the phase in the query directions. Besides the complexity in implementation, their method was limited to the scattering region near the rainbows because they encountered obstacles to calculate the phase in other scattering directions. To the best knowledge of the author, there is still no solution yet applicable to the full scattering field in 3D space by a nonspherical particle of size much larger than the wavelength of light.

As discussed in the preceding sections of this chapter, the direction, amplitude and phase for each 3D scattered ray can be calculated in the framework of VCRM. In this section, an algorithm is proposed to address the bottle-neck problem for VCRM to calculate 3D scattered intensity, namely, the interpolation of the bivariate and multiple-valued scattered data for the final interference intensity in 3D space. For a scattered ray in the 3D space, its direction is described by the azimuth angle φ and the elevation angle ψ as shown in Fig. 4.8. φ is measured in the horizontal plane from that the grid points represent the query directions in the 3D space, while the sample points represent the calculated directions of the scattered rays.

t t M \ 1 2 3 1 2 3 4 5 B A … … … … t t t t t (1) (2)
… Firstly, we need to determine which grid points are enclosed by the sample points. As illustrated in Fig. 4.11, for two adjacent layers of incident rays, say sets A 0 and B 0 , we have two sets of scattered data, say A and B. They are irregularly distributed in terms of scattering direction (φ, ψ). The band between these two sets A and B is meshed incrementally with triangles as shown in Fig. 4.11 (3). In generating a new triangle, for the sake of simplicity and topological correctness, the candidate vertexes for the front edge A i B j (i, j = 1, 2, 3, ...) are always chosen in {A i+1 , B j+1 }, and the successful one satisfies the Delaunay criterion [START_REF] Lee | Two algorithms for constructing a Delaunay triangulation[END_REF][START_REF] Lee | Generalized Delaunay triangulation for planar graphs[END_REF][START_REF] Kuo | A delaunay-based region-growing approach to surface reconstruction from unorganized points[END_REF] as shown in Fig. 4.12. In the Delaunay criterion for triangulation, the circumcircle of the meshed triangle should not enclose other vertexes. By this way, the long and thin triangles, which may deteriorate the interpolation accuracy, can be avoided to the maximum extent possible.

As shown in Fig. 4.11(3), the meshed triangles form an area enclosed by the sample points. Those grid points which lie inside the meshed triangles will be disposed of for their amplitudes and phases.

If a non-linear interpolation method is to be applied, some Delaunay refinement algorithms [START_REF] Chew | Guaranteed-quality mesh generation for curved surfaces[END_REF][START_REF] Ruppert | A Delaunay refinement algorithm for quality 2-dimensional mesh generation[END_REF][START_REF] Shewchuk | Delaunay refinement algorithms for triangular mesh generation[END_REF] are needed to refine the meshed triangles so that each triangle is approximate to an equilateral triangle, which improves the accuracy of interpolation. The scattered rays from the same order p or from different orders may be folded in certain scattering directions. As a result, the triangles meshed with the representing sample points may intersect with each other. If two or more triangles intersect at one grid point (φ, ψ) as shown in Fig. 4.14, the interference of light must be taken into account. For one of the polarizations along êθ and êφ as defined in Eq. (4-27), assume the amplitude and the phase in the direction (φ, ψ) interpolated by the j-th ( j = 1, 2...) intersecting triangle are a j and ϕ j , respectively, then the sum of complex amplitudes accounts for the interference effect:

Ẽ(φ, ψ) = ∑ j=1 a j exp(iϕ j ) (4-34)
Considering the cross polarization effect, the complex amplitude Ẽ(φ, ψ) in Eq.(4-34) should be calculated separately for the two polarizations. The total electric field in the scattering direction (φ, ψ) can then be obtained by

⃗ E(φ, ψ) =            ∑ j=1 a (1) j exp(iϕ (1) j ) ∑ j=1 a (2) j exp(iϕ (2) j )            (4-35)
where the element ∑ j=1 a (1) j exp(iϕ (1) j ) is the summed complex amplitude along the polarization direction êθ , while ∑ j=1 a (2) j exp(iϕ (2) j ) the summed complex amplitude along the polarization direction êφ .

The final scattered intensity in the scattering direction (φ, ψ) is calculated by the inner product of ⃗ E(φ, ψ) with its complex conjugate:

I(φ, ψ) = ⃗ E(φ, ψ) • ⃗ E(φ, ψ) * (4-36)
The triangulation-based interpolation algorithm proposed above is summarized as follows:

1. The incident light wave is simulated by layers of light rays.

2. The calculation is initialized for the rays from the first layer to generate the first set of scattered rays. Each scattered ray is marked as a sample point in the (φ, ψ) plane. 5. The band between these two adjacent sets A and B is meshed incrementally with triangles.

For the sake of simplicity and topological correctness, the candidate vertexes for the front edge A i B j (i, j = 1, 2, 3, ...) are always chosen in {A i+1 , B j+1 }, and the successful one satisfies the Delaunay criterion. In this way, the long and thin triangles, which may deteriorate the interpolation accuracy, can be avoided.

6. Through iteration, the triangulation of the sample points from any two adjacent layers of incident rays is carried out.

7. In each of the meshed triangles, linear interpolation is carried out for the amplitudes and the phases at the grid points which are enclosed by the concerned triangle.

8. The intersecting triangles at one grid point account for the interference effect in the represented scattering direction.

Summary

This chapter reported the first realization of a calculation method in the framework of VCRM for the 3D scattered intensity by a large nonspherical particle of any smooth surface. The ray tracing, divergence factor, phase shifts due to focal lines and optical path, and cross polarization were addressed by an elegant way using vectorial rays and the wavefront curvature. A triangulation-based interpolation algorithm has been developed to break through the bottle-neck problem for VCRM to account for the interference of scattered rays in 3D space.

The following chapter will present an application of the proposed method to the 3D scattering of plane wave by a real liquid jet and the experimental examination.

Chapter 5 The 3D scattering of plane wave by a real liquid jet and experimental examination

In principle, the method proposed in the preceding chapter can be applied to the scattering in all directions by a large and smooth particle of any shape. A real liquid jet which has a complex geometry is taken as an example here. Having a non-zero curvature along the jet axis, the geometric profile of a real liquid jet is much more complex than that of an infinite cylinder. Besides, the radius of jet's cross section usually ranges from tens to hundreds of microns, which is much larger than the wavelength of the incident laser beam (size parameter ranging approximately from 500 to 5000). For these two reasons, the calculation for the light scattering by a real liquid jet poses challenge or even an impossible task for the existing analytical theories and numerical methods. In this chapter, the calculation for the 3D far-field scattered intensity of plane wave by a real liquid jet is achieved with the proposed calculation method, which is based on VCRM while allows to calculate the scattered intensity in 3D space.

Introduction to a real liquid jet

A real liquid jet usually presents a complex profile with varying radius [START_REF] Clanet | Transition from dripping to jetting[END_REF][START_REF] Amini | Effect of gravity on capillary instability of liquid jets[END_REF][START_REF] Ambravaneswaran | Dripping-jetting transitions in a dripping faucet[END_REF][START_REF] Utada | Dripping to jetting transitions in coflowing liquid streams[END_REF][START_REF] Eggers | Physics of liquid jets[END_REF][START_REF] Duan | Influence of spatial curvature of a liquid jet on the rainbow positions: Ray tracing and experimental study[END_REF]. Even for a jet without artificial disturbance, due to the effect of surface tension and gravitational force, its transversal radius shows a continuous contraction till the breakup.

The liquid jet used here is a capillary water jet ejected in the air from a long stainless-steel needle with a circular orifice (inside and outside diameters are respectively 0.95 and 1.32 mm). A pressure of 2.1 bar is imposed on the distilled water within the reservoir; besides, a flowmeter regulates and stabilizes the flow rate at 39 ml/min. At 19.0 • C, the refractive index of the distilled water is 1.3322 for the light with wavelength λ = 632.8 nm.

The profile of the water jet near the orifice is obtained through image edge detection and data fitting as shown in Fig. 5.1. The transversal radius r as function of the distance from the orifice h is given by r(h) = 0.4997exp(-0.0183h) + 0.1377exp(-1.7990h) (

The R-square (coefficient of determination) is 0.998, so the fitting is satisfactory. By attaching the origin of Cartesian coordinate system to the exit of jet (h = 0) and setting the z axis coincident with the jet axis (z = -h), the jet profile can be expressed in the Cartesian coordinate system as r(z) = 0.4997exp(0.0183z) + 0.1377exp(1.7990z) (

where z ranges from 0 to -5.5 mm. Then the surface function of the liquid jet is given by

f (x, y, z) = x 2 + y 2 -r(z) 2 = 0 (5-3)
According to Eq. (5-1), the stream-wise curvature of the jet surface along the jet axis can be calculated. Fig. 5.2 shows the variation of the curvature with the falling height h.

One can see that the stream-wise curvature decreases exponentially from 0.4 mm -1 near the exit to 10 -3 mm -1 when h = 3.5 mm, which means the curvature radius in the stream-wise direction is very small (2.5 mm) at the beginning but is very large (several meters) when the falling distance h is more than 3.5 mm.

Due to the presence of stream-wise curvature, the scattering field by a real liquid jet is different from that by an infinite cylinder. Firstly, the light rays scattered by a liquid jet are no longer within one transversal section, instead, they are spread into the 3D space, as illustrated in Fig. 5.3. Secondly, since the scattered light of different orders are elevated differently from the horizontal plane, the scattering field by a liquid jet shows an interference pattern against common sense as shown in Fig. 5.4. The scattered light of different orders are separated from each other, but interference still exists within the single-order scattered light of p ≥ 2.

In the following, the 3D scattering field of plane wave by the capillary water jet modeled in Eq. (5-1) will be simulated by the calculation method proposed in the preceding Chapter 4.

Simulation of the 3D scattering field

The calculation for the 3D scattered intensity by a real jet is far beyond the capabilities of most electromagnetic approaches because:

1. A real jet has a geometric profile being so complex that the separation of variable methods such as LMT and Debye cannot be applied.

2. The size parameter of the concerned jet is nearly 5000, which requires unimaginable computer memories and computation time for the numerical methods such as FDTD, T-matrix and DDA.

With the calculation method proposed in Chapter 4, which is based on VCRM while allows to calculate the scattered intensity in 3D space, the simulation for the 3D scattering field of plane wave by a real jet is successfully achieved here.

Before showing the intensity distribution of the scattering field, the result of ray tracing is presented firstly to offer readers an intuitive understanding of the ray paths.

Based on the ray tracing method discussed in Sec. 4.2, the tracing of the light rays scattered by the jet is shown in Fig. 5.5. The refractive index of jet is 1.3322; and the surface function has been given in Eqs.

(5-1) and (5-3). As a comparison, the tracing for the light rays scattered by an infinite circular cylinder is also presented. The infinite cylinder has the same refractive index as the jet, but its cross section has a constant radius of r 0 (the initial radius of the jet). The direction of the incident rays is perpendicular to the jet/cylinder axis. For clarity, only one layer of the incident rays and the corresponding scattered rays of p ≤ 2 are shown.

From Fig. 5.5, one can see that for the infinite cylinder, its scattered rays stay in the same plane as the incident rays; while for the real jet, due to the non-zero curvature along the jet axis, the emergent rays are scattered into 3D space. More specifically, because the normal vectors of jet surface point downwards, the externally reflected rays (p = 0) are tilted downwards, while the scattered rays of order p ≥ 1 are tilted upwards. Moreover, the higher the order is, the greater the elevation from the horizon will be (as illustrated in Fig. 5.3).

Through ray tracing, one can obtain the directions of the scattered rays. Furthermore, according to Secs. 4.3-4.5, the corresponding amplitudes and phases for these 3D scattered rays can be calculated. And, through the triangulation-based interpolation algorithm proposed in Sec. 4.6, one can acquire the complex amplitude and the final interference intensity in each scattering direction of the 3D space. One the other hand, it is found that as the incident beam moves downwards, i.e. when h is increased, the separation of scattered light along the ψ axis is lessen as shown in Fig. The simulation result by the calculation method based on VCRM for 3D scattering will be examined by an experiment in the following section.

Experimental setup

To validate the calculation method proposed in the preceding Chapter 4 and the simulation for the 3D scattered intensity by a real jet, an experiment is carried out. The setup is shown in Fig. 5.9, while the schematic diagram is given in Fig. 5.10. The experimental setup mainly contains three parts:

• Light source with beam shaping module;

• Generation and control of the liquid jet;

• Image acquisition part. in the far field), and is then recorded by the camera #1 (HAMAMATSU, C9100-022 ). The recorded scattered intensity near the first-order and second-order rainbows is shown in Fig.

Comparison of simulated result with that by experiment

The following part covers the comparison of the calculated 3D scattered intensity with the one measured by experiment.

Since the intensity of first-order rainbow is much greater than that of second-order rainbow, the scattering patterns near the two rainbows are measured separately with different exposure time. The comparison between the simulation and the experiment for the 3D scattered intensity near the first-order rainbow is given in Fig. 5.12(a), while the comparison near the second-order rainbow is given in Fig. 5.12(b). For clarity, near the first-order rainbow, the examination is made in three slices: ψ = 1.9 • , 2. Two notable features are observed: (1) the scattered light of different orders are naturally separated but the "low-frequency" Airy fringes still exist because of the interference of the scattered rays from the same order; (2) the scattered intensity varies also with the elevation angle ψ, which is attributed to the variation of the surface curvature along the jet axis. In this aspect, the jet acts much like a concave lens (see the jet profile shown in Fig.

5.1).

We can see that a good agreement is found not only in the direction of φ, where the simulated peaks and troughs (rainbow fringes) match the experimental results very well, but also in the direction of ψ, where the simulated intensity as function of the elevation angle ψ also fits well with that by experiment.

The main discrepancy is at the rainbow angles (caustics). It is caused by the intrinsic defect of ray model, which predicts an infinite intensity there [6,[START_REF] Hovenac | Assessing the contributions of surface waves and complex rays to far-field Mie scattering by use of the Debye series[END_REF][START_REF] Onofri | Experimental validation of the vectorial complex ray model on the inter-caustics scattering of oblate droplets[END_REF][START_REF] Laven | Supernumerary arcs of rainbows: Young's theory of interference[END_REF]. To improve this, the physical optics [START_REF] Sentis | Scattering of light by large bubbles: Coupling of geometrical and physical optics approximations[END_REF][START_REF] Nye | Optical caustics in the near field from liquid drops[END_REF][START_REF] Onofri | Physical-optics approximation of near-criticalangle scattering by spheroidal bubbles[END_REF] is to be resorted to. Handling with this flaw existing in ray optics needs much more effort and is not within the scope of the present thesis.

Summary

The calculation method proposed in Chapter 4, which is based on VCRM while allows to calculate the scattered intensity in 3D space, has been successfully applied to solving the 3D scattered intensity by a real liquid jet. Taking advantage of the ability of VCRM for interpreting the scattering mechanism, a systematic analysis has been made for the scattered light of different orders, in regard to their separation or interference in 3D space.

The computation was completed on a laptop (Intel i7-8550U @1.80GHz), and the calculation for obtaining the 3D scattering field of p ≤ 3 took about 20 minutes. Though the performance is not the focus of this thesis, the triangulation and interpolation processes might speed up in future by applying parallel computation, which I believe will reduce the time to minutes, or even seconds.

Experiment has also been carried out to verify the proposed method and to verify the simulated results. An incident wave of limited width 1 mm which has approximately uniform intensity and planar wavefront has been achieved in experiment, thus complying with the incident plane wave with sheet spot used in simulation. Good agreements have been found not only on the spatial separation of different scattering orders but also on the interference intensity and the angular extension ranges.

In the following chapter, the ray model of incident elliptical Gaussian beam will be established in the framework of VCRM, and the 3D scattered intensity of elliptical Gaussian beam by a real liquid jet will be solved and analyzed.

Chapter 6 The 3D scattering of elliptical Gaussian beam by a real liquid jet

The light beam emitted from a laser device generally has a Gaussian profile, thus making important the research on the scattering of Gaussian beam by particles. Many researchers have contributed to the development and application of the GLMT for the scattering of a Gaussian beam by particles. However, the GLMT is applicable only to those particles of simple forms such as spherical or spheroidal particles, or circular or elliptical cylinders.

In this Chapter, the calculation for the 3D scattered intensity of elliptical Gaussian beam by a large nonspherical particle of any smooth shape is achieved by VCRM.

Ray model of elliptical Gaussian beam

Compared to a plane wave, the light rays which characterize an elliptical Gaussian beam1 are different in direction, amplitude, phase and wave-front curvatures. In this section, the calculation for these properties will be addressed.

Consider an elliptical Gaussian beam with simple astigmatism [START_REF] Arnaud | Gaussian light beams with general astigmatism[END_REF][START_REF] Rohani | A fast Gaussian beam tracing method for reflection and refraction of general vectorial astigmatic Gaussian beams from general curved surfaces[END_REF][START_REF] Shishegar | Scattering of a Gaussian beam by a spherical dielectric object using the vectorial Gaussian beam tracing method[END_REF][START_REF] Rohani | A Fast Hybrid Method For Analysis and Design of Photonic Structures[END_REF][START_REF] Shabani | Vectorial Gussian beam expansion for high-frequency wave propagation[END_REF] which propagates along the ŵ axis of the coordinate system (û, v, ŵ) as shown in Fig. 6.1, the incident electric field can be expressed by

⃗ E G (u, v, w) = ⃗ E 0 E G exp(iϕ), (6-1)
where

⃗ E 0 = E 0u û + E 0v v.
The ratio E 0u : E 0v determines the vibration direction. For example, the case with E 0u = 0 indicates the wave is polarized along v. E G is the normalized amplitude, while ϕ is the phase. They are given respectively by [START_REF] Jiang | Theoretical study of light scattering by an elliptical cylinder[END_REF][START_REF] Arnaud | Gaussian light beams with general astigmatism[END_REF][START_REF] Rohani | A fast Gaussian beam tracing method for reflection and refraction of general vectorial astigmatic Gaussian beams from general curved surfaces[END_REF][START_REF] Shishegar | Scattering of a Gaussian beam by a spherical dielectric object using the vectorial Gaussian beam tracing method[END_REF][START_REF] Rohani | A Fast Hybrid Method For Analysis and Design of Photonic Structures[END_REF][START_REF] Shabani | Vectorial Gussian beam expansion for high-frequency wave propagation[END_REF]]

E G (u, v, w) = √ r 0u r u √ r 0v r v exp ( - u 2 r u 2 - v 2 r v 2 ) (6-2) and ϕ(u, v, w) = - 2π λ ( w + u 2 2R u + v 2 2R v ) + 1 2 tan -1 ( w -O u l u ) + 1 2 tan -1 ( w -O v l v ) (6-3)
where r 0u is the beam-waist radius of the elliptical Gaussian beam in the plane containing ŵ O v present as saddle surfaces. Only by passing both O u and O v , the wavefronts are diverging in both the ( ŵ, û) and the ( ŵ, v) planes.

In Fig. 6.2(b), the vectors normal to the wavefronts are calculated by Eq. (6-9), representing the propagation directions of the wave at different points. Different from a plane wave, the wave directions of a shaped beam cannot be characterized by a bundle of straight rays in a strict sense. Instead, they are calculated at the concerned points, for example at the incident points.

The expression of k in the Cartesian coordinates sytem ( x, ŷ, ẑ) can be obtained by

k = - ϕ ′ u û + ϕ ′ v v + ϕ ′ w ŵ ∥ϕ ′ u û + ϕ ′ v v + ϕ ′ w ŵ∥ (6-10)
To obtain the two principal curvatures as well as the corresponding principal directions of incident elliptical Gaussian beam at a concerned point (u, v, w), two irrelevant base vectors in the tangent plane of iso-phase surface are defined as follows:

α = v × k/∥v × k∥ (6-11) β = k × α (6-12)
Further, three scalars are defined as: -15) where the superscript T indicates the transposition. H(ϕ) is the Hessian matrix of the phase function defined by

ϕ αα = α • H(ϕ) • αT ; (6-13) ϕ αβ = α • H(ϕ) • βT ; (6-14) ϕ ββ = β • H(ϕ) • βT . ( 6
H(ϕ) =           ϕ ′′ uu ϕ ′′ uv ϕ ′′ uw ϕ ′′ vu ϕ ′′ vv ϕ ′′ vw ϕ ′′ wu ϕ ′′ wv ϕ ′′ ww           (6-16)
These second-order derivatives of the phase function with respect to u, v and w are calculated iso-phase surface.

Till now, the amplitude, phase, direction, principal curvatures and principal directions of each light ray which characterizes a local wavefront of an elliptical Gaussian beam can be obtained. These light rays illuminate the concerned particle; and then, they are subjected to continual reflections and refractions by the particle.

It should be noted that the coordinate system (û, v, ŵ) which describes the beam may be different from the coordinate system which describes the particle. Suppose the particle is described in the Cartesian coordinate system ( x, ŷ, ẑ) and centered at O : (0, 0, 0), while coordinate system (û, v, ŵ) which describes the Gaussian beam is centered at O

G : (x G , y G , z G )
as shown in Fig. 6.4. The position vector of a point Q : (x 0 , y 0 , z 0 ) in the system of ( x, ŷ, ẑ) is expressed as --→ OQ = x 0 x + y 0 ŷ + z 0 ẑ (6)(7)(8)(9)(10)(11)[START_REF] Li | Measurement of transient evaporation of an ethanol droplet stream with phase rainbow refractometry and high-speed microscopic shadowgraphy[END_REF][START_REF] Li | Simultaneous measurement of refractive index, diameter and colloid concentration of a droplet using rainbow refractometry[END_REF][START_REF] Langley | Critical-angle scattering of laser light from bubbles in water: measurements, models, and application to sizing of bubbles[END_REF][START_REF] Hansen | Mie scattering as a technique for the sizing of air bubbles[END_REF][START_REF] Arnott | Unfolded optical glory of spheroids: backscattering of laser light from freely rising spheroidal air bubbles in water[END_REF][START_REF] Zhang | Influence of bubbles on scattering of light in the ocean[END_REF][START_REF] Mougin | Path instability of a rising bubble[END_REF][START_REF] Lindken | A novel PIV technique for measurements in multiphase flows and its application to two-phase bubbly flows[END_REF][START_REF] Fujiwara | Bubble deformation and flow structure measured by double shadow images and PIV/LIF[END_REF][START_REF] Li | Influence of characteristics of micro-bubble clouds on backscatter lidar signal[END_REF][START_REF] Onofri | Near-criticalangle scattering for the characterization of clouds of bubbles: particular effects[END_REF][START_REF] Cubaud | Bubble dispenser in microfluidic devices[END_REF][START_REF] Yang | Gas content effect on bubble motion[END_REF] where x 0 , y 0 and z 0 are the projections of --→ OQ on x, ŷ and ẑ, respectively. Then, the vector

from O G to the point Q, ---→ O G Q, is calculated in the system of ( x, ŷ, ẑ) by --→ OQ - ---→ OO G . That is: ---→ O G Q = (x 0 -x G ) x + (y 0 -y G )ŷ + (z 0 -z G )ẑ (6-25)
Then, the components of ---→ O G Q along û, v and ŵ are obtained through:

             ---→ O G Q • û = (x 0 -x G ) x • û + (y 0 -y G )ŷ • û + (z 0 -z G )ẑ • û ---→ O G Q • v = (x 0 -x G ) x • v + (y 0 -y G )ŷ • v + (z 0 -z G )ẑ • v ---→ O G Q • ŵ = (x 0 -x G ) x • ŵ + (y 0 -y G )ŷ • ŵ + (z 0 -z G )ẑ • ŵ (6-26)
Finally, the point Q of coordinates (x 0 , y 0 , z 0 ) in the particle system ( x, ŷ, ẑ) is expressed in the beam system (û, v, ŵ) as:

          u 0 v 0 w 0           =           (û • x) (û • ŷ) (û • ẑ) (v • x) (v • ŷ) (v • ẑ) ( ŵ • x) ( ŵ • ŷ) ( ŵ • ẑ)                     x 0 -x G y 0 -y G z 0 -z G           (6-27)
On the other hand, suppose the direction of a light ray in the beam system is:

k = k u û + k v v + k w ŵ (6-28)
where k u , k v and k w are the projections of the normalized wave vector k on û, v and ŵ, respectively. In the particle system ( x, ŷ, ẑ), the unit vector which has the same direction as k is obtained by

             k • x = k u û • x + k v v • x + k w ŵ • x k • ŷ = k u û • ŷ + k v v • ŷ + k w ŵ • ŷ k • ẑ = k u û • ẑ + k v v • ẑ + k w ŵ • ẑ (6-29)
Then, the direction of the concerned ray in the particle system ( x, ŷ, ẑ) is obtained as

k = k x x + k y ŷ + k z ẑ (6-30)
where the k x , k y and k z are calculated according to

          k x k y k z           =           ( x • û) ( x • v) ( x • ŵ) (ŷ • û) (ŷ • v) (ŷ • ŵ) (ẑ • û) (ẑ • v) (ẑ • ŵ)                     k u k v k w           (6-31)

Phase shift due to optical path

For a particle illuminated by an elliptical or a circular Gaussian beam, the phase shifts of a scattered ray due to reflections and focal lines can be calculated in the same way as that for the scattering of a plane wave, which have been discussed in detail in Secs. 3.2.3 and 4.5, and will not be repeated here. However, it is worth noting that for a shaped beam, the incident light rays which characterize the beam are different in phase. According to Eq.

(6-3), the phase is dependent on the spatial position (u, v, w), indicating that the formula for the phase shift due to optical path proposed in Eq. (3-25) needs an amendment when a shaped beam is used as the light source.

The Gouy shift terms tan -1 ( w-O u l u ) and tan -1 ( w-O v l v ) in Eq. (6-3) can be interpreted as a phase difference between the elliptical Gaussian beam and a plane wave of equal wavelength and identical propagation direction [START_REF] Boyd | Intuitive explanation of the phase anomaly of focused light beams[END_REF]. Thus, the reference ray is chosen such that it has a planar wavefront propagating in the ŵ direction (axial direction of the Gaussian beam); after passing through the particle center, the reference ray emerges out in the same direction as the concerned scattered ray. On the particle surface, suppose an illuminated point W 0 : (x 0 , y 0 , z 0 ) as shown in Fig. 6.5. The coordinates of W 0 in the beam system, (u 0 , v 0 , w 0 ), can be obtained according to the Eq. . Then, the phase of the incident ray which characterizes the elliptical Gaussian beam at W 0 is calculated by substituting (u 0 , v 0 , w 0 ) into Eq. (6-3):

V V V V V V V V V V V V V V V V V V V V V
ϕ i = ϕ(u 0 , v 0 , w 0 ) (6-32)
The phase shift due to optical path for a p-order scattered ray at its emergent point is calculated according to the phase at the incident point ϕ i and the optical path it has experienced inside the particle, that is

ϕ p = ϕ i - 2π λ m p-1 ∑ j=0 s j (6-33)
For an externally reflected ray (p = 0) at its emergent point, ϕ p is simply equal to the ϕ i , since it has no extra optical path compared to the incident ray. s j is the geometric length from the interaction point W j to W j+1 . The negative sign indicates that the phase lags for a longer optical path. The emergent point of a p-order scattered ray is denoted as W p , while the emergent direction is denoted as kp .

The phase of the reference ray is calculated at N p , where N p is a point on the reference ray and satisfies ----→ W p N p ⊥ kp , as

ϕ re f = - 2π λ σ p (6-34)
where σ p = ---→ OW p • kp , as illustrated in Fig. 6.5.

Finally, the phase shift due to optical path of a p-order scattered ray relative to the reference ray is obtained as

ϕ p,OP = ϕ p -ϕ re f = ϕ(u 0 , v 0 , w 0 ) + 2π λ         σ p -m p-1 ∑ j=0 s j         (6-35)

Comparison with GLMT for a spherical particle

Before the calculation for particles of complex geometry, the proposed method is firstly examined by comparing with GLMT for the scattering by a spherical particle.

The radius of the spherical particle is set as 25 µm, while its refractive index relative one needs only to calculate the scattering of the light rays inside one scattering plane, for example, the scattering in the xy plane. The transverse electric (TE) polarization is under consideration, that is, the incident electric vector is linearly polarized along the z axis (vertical to the paper plane), with E 0z = 1 and E 0y = 0. In VCRM, the scattered rays of orders p ≤ 7 are taken into account. One can see that the agreement is rather satisfying. Since the effect of diffraction is not within the scope of the current thesis, so comes in the difference in the forward region (θ ≤ 2.5 • ). The size parameter of the spherical particle in Fig. 6.7 is only about 250; for larger particles, the agreement would be better.

Besides, it should be noted that for such a tightly focused beam, the characterization by rays may bring about inaccuracy, especially near the beam waist. For this reason, in the following calculations, the beam waist in the plane where the beam is tightly focused is set far away from the particle.

Although the VCRM provides only an approximation solution to the light scattering, it has distinct advantage over the GLMT in the aspect of nonspherical particles: the VCRM only requires the particle surface to be smooth. The following parts of this chapter deal with the 3D scattering of elliptical Gaussian beam by a real liquid jet of complex shape.

Scattering of an elliptical Gaussian beam by a real jet

Here we take the liquid jet modeled in Sec. 5.1 as the scatterer. As has already been discussed, the light scattering by a real liquid jet can hardly be solved by the analytical models or the existing numerical methods. In the following parts of this chapter, the scattering of elliptical Gaussian beam by a real liquid jet in 3D space will be simulated and analyzed based on VCRM.

The research on the scattering of elliptical Gaussian beam by a liquid jet is motivated primarily by the need for measuring the jet properties (size, refractive index or instability) by optical method. The beam-waist radius of the incident Gaussian beam in the plane containing the jet axis is preferred to be small, so that the scattering light from different parts of the jet can be distinguished. A typical configuration for such an elliptical light spot is presented in Fig. 6.8(a).

The jet axis coincides with the z axis of Cartesian coordinate system, while the beam axis is defined as the x axis (perpendicular to the paper plane in Fig. 6.8(a)). The incident beam is linearly polarized along the z axis (E 0z = 1 and E 0y = 0) 3 . The relative refractive index of the jet is 1.3322 for wavelength λ = 632.8 nm. The beam-waist radius along the jet axis, r 0z , is 0.4 mm, while the beam-waist radius along the ŷ axis, r 0y , is 1.6 mm. The divergence angle in the xz (vertical) plane, γ z = λ/πr 0z , is 0.029 • , while the divergence angle in the xy (horizontal) plane, γ y = λ/πr 0y , is 0.007 • . The beam centers in the horizontal and vertical planes, O y and O z are set at x = 0 (on the jet axis) and x = -30 cm, respectively.

The distance of the center of light spot from the exit, denoted as h, is 3 mm as shown in Fig. 6.8(a). At the incident points, the directions, phases and amplitudes of wave, and the principal curvatures and principal directions of wavefront are determined according to Sec. 6.1. Then, based on the calculation method proposed in Chapter 4, the direction, phase, amplitude, polarization and divergence factor for each scattered ray and, the final interference intensity of all scattered rays in the 3D space can be obtained. Since the radial range of a Gaussian beam has no boundary in theory, we need to truncate the radial range and just retain the main part of energy in calculation. Here, the radial range of the incident Gaussian beam is truncated at 1.5 times of the beam-waist radius,

where the light intensity falls to 1.11% of its axial value. This truncation sets a limit to the observation range (the scattered light has a limited extent along the ψ axis), but has little effect on the precision of results since only the trivial part of energy is cut off.

As shown in Fig. 6.8(b), the elevation angle ψ of the p = 0 scattered light ranges approximately from -1.2 • to -0.9 • . Besides, because of the Gaussian profile of incident beam, the energy of the p = 0 scattered light is mainly distributed at ψ = -1.0 • .

Considering that the divergence angle of the incident beam is only 0.029 • in xz plane, the elevation angle of scattered light is mainly related to the surface curvature of jet along the z axis. In this aspect, the scattering of an elliptical Gaussian beam with very small divergence angle is analogous to the scattering of a wave of planar wavefront. However, there is notable difference here because an elliptical Gaussian beam, even with a near-zero divergence angle, has non-uniform intensity. The Gaussian profile of intensity leads to the Gaussian-like distribution of the rainbow fringes along the axis of ψ, as shown in Fig. 6.8(b).

Please notice the difference from the rainbow fringes for a wave of planar wavefront and uniform intensity given in Fig. 5.8.

As the incident Gaussian beam illuminates the jet at a different position, the corresponding scattering field shall be altered because the stream-wise curvature of jet surface is changed. Fig. 6.9 shows the scattering field when the incident position moves down to h = 3.5 mm. It is found that as the distance of beam center from the exit of orifice is increased, the deviation of scattered light from the horizontal plane (ψ = 0 • ) becomes smaller.

For example, the p = 3 scattered light for h = 3.0 mm (shown in Fig. 6.8(b)) lies in the region where ψ ranges approximately from 2.6 • to 3.6 • , while the p = 3 scattered light for h = 3.5 mm (shown in Fig. 6.9) is less elevated with ψ ranging approximately from 2.5 • to 2.9 • . This is attributed to the decrease of the stream-wise curvature of jet surface. Besides, as the expansion along the ψ axis becomes shrunken, the intensity of the scattered light is greater. For example, even without amplification, the second-order rainbows in Fig. 6.9 can be observed clearly.

Furthermore, the 3D scattered intensity of a circular Gaussian beam by the jet is calcu-beam of small divergence. In this section, the influence of divergence angle on the scattering field is investigated.

Firstly, consider an elliptical Gaussian beam which propagates along the x axis has a considerable divergence angle in the plane containing the jet axis (xz plane). The beamwaist radius in the xz plane, r 0z , is set as small as 8 µm (divergence angle γ z = 1.44 • ).

The beam profile in this plane is shown in Fig. 6.11(a). To focus on the influence of the divergence angle in the vertical xz plane, the divergence angle in the horizontal xy plane, γ y , is fixed as 0.007 • , the same as that in Figs. rays which characterize the elliptical Gaussian beam are of considerable divergence in the xz plane, that is, the z-components of wave vectors become significant. Consequently, the elevation angles of scattered rays, which are directly related to the z-components of the scattered wave vectors, become much different from those in the scattering of an elliptical Gaussian beam with small divergence angle. When compared to the results in Fig. 6.8(b) and Fig. 6.9 where γ z = 0.029 • , the scattering field for γ z = 1.44 • as shown in Fig. 6.12 has a notable feature: the scattered light is expanded along the ψ axis and the scattered light of orders p = 2 and p = 3 overlap with the light of p = 0. The comparison of Fig. 6.12 with Figs. 6.8(b) and 6.9 indicates that the superposition or separation of the scattered light of different orders can be regulated by changing the divergence angle of the incident beam. To further prove this, Fig. 6.13 presents the calculated scattering fields for two different divergence angles. In Fig. 6.13(a), the beam-waist radius in xz plane r 0z is 6 µm and the corresponding divergence angle γ z is about 1.92 • , while in Fig. 6.13(b) the beam-waist radius r 0z is 16 µm and γ z is decreased to 0.72 • . The other parameters are the same as those in Fig. 6.12. angle γ z is lessened, the scattered light of different orders are separated in the way similar to the case when a plane wave illuminates the jet, where the stream-wise curvature of jet surface becomes the dominant factor for the separation.

From the comparison in

Spatial characteristics of the 3D scattering field

It is well known that the angular distribution of the scattered intensity near the firstorder rainbow of a spherical droplet can be decomposed into two components: the "lowfrequency" component (referred to as the Airy structures) and the "high-frequency" component (Ripple structures), which are attributed to the interference of the p = 2 scattered rays themselves and the interference of the p = 0 rays with the p = 2 rays, respectively. This decomposition can be achieved by using a femtosecond laser [START_REF] Laven | Time domain analysis of scattering by a water droplet[END_REF] where the scattered light of different orders are separated in the time domain. With the calculation method for 3D scattering (Chapter 4), the calculation for the scattered light by a large nonspherical particle can be achieved in the 3D space, which enables us to separate the scattered light of different orders in the space domain.

As discussed in the preceding Secs. 6.4 and 6.5, by adjusting the divergence angle of the incident elliptical Gaussian beam, the scattered light of different orders from a real liquid jet may be separated from or superpose with each other in the 3D space. Here, the scattering field given in Fig. 6.13(a) is taken for an example. The zoomed view of the scattered intensity near the first-order rainbow is shown in Fig. 6.14. For the convenience of the readers, the parameters for the elliptical Gaussian beam are tabulated in Table 6.6.

Table 6.6 Parameters of the elliptical Gaussian beam in Figs. 6.13(a) and 6.14. In the region where ψ > 2.0 • , the rainbow fringes result from the interference within the p = 2 scattered rays. Thus, only the Airy structures are formed in this region. The intensity along the sampled line (a) where ψ = 2.2 • is shown in Fig. 6.16(a). One can see that the first-order rainbow contains only the "low-frequency" Airy structures for ψ = 2.2 • .

Beam axis

As to the overlapping region, the rainbow fringes are formed by the superposition between the p = 2 and the p = 0 scattered rays. The intensity along the sampled line (b)

where ψ = 1.0 • is shown in Fig. 6.16(b). The notable difference from the interference fringes in Fig. 6.16(a) is the appendage of the "high-frequency" Ripple structures on the "low-frequency" Airy structures. But, the Ripple structures are still not distinct because at this elevation angle, the intensity of the p = 0 scattered rays is still much lower compared to the p = 2 rays as given in Fig. 6.15. As we go downward further till ψ = 0 • , the intensity of the p = 2 rays is decreased, while the intensity of the p = 0 rays is increased. Consequently, their interference effect becomes more distinct as shown in Fig. 6.16(c). Furthermore, according to Fig. 6.15, the intensity of the p = 2 rays is roughly equal to that of the p = 0 rays when ψ = -1.0 • . Along the sampled line (d) where ψ = -1.0 • , the rainbow fringes are presented in Fig. 6.16(d). One can see that the interference effect between p = 2 and p = 0 rays is the most distinct for ψ = -1.0 • .

It concludes from the Figs. 6.14 and 6.16 that the separation and the superposition of the scattered light of different orders can be coexisting in the scattering of elliptical Gaussian beam by a real liquid jet. This analysis in the space domain helps to have a better understanding of the 3D scattering of shaped beam by a large nonspherical particle.

Summary

In this chapter, the ray model for incident elliptical Gaussian beam has been proposed in the framework of VCRM. Not only the propagation direction, amplitude and phase, but also the principal curvatures and principal directions at each point of the incident wavefront have been successfully characterized by light rays. It is now possible to calculate the 3D scattered intensity of laser beam by a large particle of any smooth surface based on VCRM.

After the examination by comparing with GLMT for the scattering by a spherical particle of different refractive indices, the proposed method was then applied to solving and analyzing the 3D far-field scattered intensity of elliptical Gaussian beam by a real liquid jet of complex shape. It has been found that by changing the divergence angle of the incident elliptical Gaussian beam, the scattered light by a real jet for different orders p can be separated from or interfere with each other, resulting in different scattering patterns. A quantitative analysis has also been made for the spatial characteristics of the scattering field when a tightly focused beam illuminates the jet.

It is wished to promote the understanding of the 3D scattering of laser beam by a large nonspherical particle and to provide methodology for solving the 3D scattered intensity of other shaped beams by large nonspherical particles under the framework of VCRM.

Chapter 7 Conclusions and perspectives

This thesis is devoted to the development of vectorial complex ray model (VCRM) for solving the 2D and 3D scattered intensity of plane wave or shaped beam by large particles of any smooth surface. This chapter draws the conclusions of the current thesis and gives perspectives in future studies.

Conclusions

For a nonspherical particle of size much larger than the wavelength (size parameter ranging from several hundreds to thousands or even larger), the calculation for the scattered intensity of light is a formidable or even impossible task using the exact numerical methods such as DDA, FDTD and T-matrix. One feasible way is to approximate the incident light wave as a bundle of rays and then, to study the interaction of rays with the particle. In a GOA method, the wave-front curvature of light wave is not included, making it difficult to calculate the shift in amplitude and the shift in phase due to the divergence and convergence of wavefront. In the VCRM, developed recently in the Coria institute & University of Rouen Normandy, the concept of wave-front curvature is introduced to the ray model of light to describe the wave-front divergence or convergence. Moreover, the ray directions and the Fresnel coefficients are calculated by the wave vectors and their components, which considerably simplifies the calculation when a nonspherical particle is involved.

For the flexibility and efficiency merits of VCRM, this thesis aimed to extend the VCR-M to solving the following issues and has achieved substantial progress:

1. The 2D light scattering by an infinite cylinder of arbitrarily smooth cross section. In the high-frequency limit, this thesis has provided a solution, based on VCRM, to the whole scattering field by an infinite cylinder of any cross section. The proposed method is flexible and, in principle, only requires the cross section of cylinder to be smooth. It permits now to calculate and analyze the scattering characteristics of light by an infinite cylinder with cross section ranging simple to complex.

Based on the proposed method, a numerical study has been performed on the light scattering by composite elliptical cylinders (CECs), whose cross sections approximate the shapes of natural raindrops. The effects of shape deformation, refractive index and the direction of incident wave on the scattering patterns of CECs have been investigated and quantitatively analyzed, which provided insight into how these factors affect the appearance of a natural rainbow.

2. The 3D light scattering by a large particle of any smooth surface. This thesis reported the first realization of a calculation method in the framework of VCRM for the 3D scattered intensity by a large nonspherical particle. The problems such as the ray tracing, the calculation of the divergence factor, the phase shifts due to focal lines and optical path, and the cross polarization effects encountered in calculating the 3D scattering field have been addressed. Moreover, a triangulation-based interpolation algorithm has been proposed

to resolve the bottle-neck problem for VCRM to address the interference of light in 3D

space. The proposed calculation method is directly applicable to calculating the scattered intensity in any direction of 3D space by a large particle of any smooth shape.

3. The simulation of the 3D scattering field by a real liquid jet and experimental examination. In this thesis, the calculation for the 3D far-field scattered intensity of plane wave by a real liquid jet has been achieved with the proposed calculation method, which is based on VCRM while allows to calculate the scattered intensity in 3D space. It has been found

that due to the stream-wise curvature of jet surface, the scattered light of different orders are naturally separated in the 3D space. Consequently, the interference field shows a very interesting pattern that has not been found in the scattering by an infinite cylinder or a

sphere.

An experiment has also been carried out to verify the calculation method proposed for 3D scattering and to verify the calculation result. By using a 4 f optical system and a spatial filter, an incident wave of limited width 1 mm and of approximately uniform intensity and planar wavefront has been generated, which complied with the incident wave used in simulation. Good agreements between experiment and simulation have been found not only on the spatial separation of different scattering orders but also on the scattered intensity and the angular extension ranges.

4. The 3D scattering field of elliptical Gaussian beam by a real liquid jet. A ray description method for incident elliptical Gaussian beams has been proposed in the framework of VCRM, thus making it possible to calculate the 3D scattered intensity of laser beam by a large particle of any smooth surface. The propagation direction, amplitude, phase, principal curvatures and principal directions at each point of the incident elliptical Gaussian beam are characterized with a light ray. Then, all the light rays that represent the incident elliptical Gaussian beam illuminate the particle and are scattered to all possible directions of 3D space. After the validation by comparing with GLMT for the light scattering by spherical particles, the 3D far-field scattered intensity of elliptical Gaussian beam by a real jet which has a complex shape has been simulated and analyzed.

The influence of the beam divergence angle on the spatial distribution of the scattering field near the first-and second-order rainbows has been investigated. It has been found

that by changing the divergence angle of the incident elliptical Gaussian beam, the scattered light by a real jet for different orders p can be separated from or superpose with each other, resulting in different interference patterns. A quantitative analysis has also been made for the spatial characteristics of the scattering field when a tightly focused beam illuminates the jet.

Perspectives

One major innovation of the current thesis is the first realization of an algorithm in the framework of the VCRM to calculate the 3D scattered intensity, in the high-frequency limit, by a smooth particle of any shape. Although the proposed method was devoted to the 3D scattered intensity by a real liquid jet in this thesis, the calculation for the 3D scattered intensity of plane wave or elliptical Gaussian beam by other nonspherical particles, natural raindrops/droplets for example (see Fig. 7.1), should not be too difficult. Another practical application is the calculation for the temporal-spatial distribution of the light rays transmitted from sea into air or from air into sea. In the presence of wind, the sea surface shows a wavy profile of ups and downs as shown in Fig. 7.2. The JONSWAP spectrum [START_REF] Yang | Refraction correction of airborne lidar bathymetry based on sea surface profile and ray tracing[END_REF] is used here. The wind region is set as 20 km. For wind speed being 4. More work is needed in future to obtain the polarization and the 3D intensity distribution of the coherent/incoherent light transmitted from sea surface. The author believes these issues could be well solved with the calculation method proposed in this thesis.

The calculation method proposed in this thesis for the 2D and 3D scattered intensity by large nonspherical particles may serve as the foundations for addressing the inverse problems, namely the measurements for the size, geometry and refractive index of relevant particles.

Besides the perspective applications, the fundamental researches on the physical problems of the ray model of light (including the GOA and VCRM) are also of significant importance. One challenging problem in the ray model of light is to calculate the light intensity near caustics. To achieve this, the combination of geometrical optics with physical optics is one interesting and significant work in the future. Another problem to be solved is the inclusion of surface waves. The author believes that the breakthrough of these two bottle-neck

Fig 2 .

 2 Fig 2.1 Schematic diagram of an infinite circular cylinder illuminated by a plane

1

 1 Fig 2.1 Schematic diagram of an infinite circular cylinder illuminated by a plane wave. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Fig 2 . 2 Fig 2 . 3 Fig 2 . 4

 222324 Fig 2.2 Ray path in the transversal section of circular cylinder. . . . . . . . . . . . Fig 2.3 Schematic diagram of the divergence of a ray pencil interacting with a circular cylinder. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fig 2.4 Schematic diagram of the virtual rays (dotted) serving as the references in calculating the optical paths of scattered rays. . . . . . . . . . . . . . . . .

Fig 2 . 5

 25 Fig 2.5 Two types of focal lines defined by van de Hulst [1] in the scattering by a sphere. For a circular cylinder concerned in this section, only the focal lines of type A exist. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Fig 2 . 6

 26 Fig 2.6 Comparison of the scattered intensity calculated by GOA with that by LMT for different size parameters. m = 1.3322. . . . . . . . . . . . . . . . . . .

Fig 2 . 7

 27 Fig 2.7 The total and the single-order scattered intensity by a circular cylinder for two polarizations. For clarity, only the scattered light of p = 0, 1, 2 and 3 are presented. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Fig 2 . 8 Fig 2 . 9

 2829 Fig 2.8 First-order rainbow calculated by GOA for an infinite circular cylinder. The calculation parameters are the same with those in Fig. 2.7(a). . . . . . . . . Fig 2.9 Comparison of the scattered intensity calculated by GOA with that by LMT for three spherical particles of different radii. m = 1.3322. . . . . . . . . .

Fig 2 .

 2 Fig 2.10 Surface wave contributions for two spherical particles of different radii. . .

Fig 2 . 11

 211 Fig 2.11 Zoomed view of Fig. 2.9(c) in linear scale for the first-order rainbow calculated by GOA and LMT. . . . . . . . . . . . . . . . . . . . . . . . . . . .

Fig 3 . 1

 31 Fig 3.1 An infinite cylinder of cross section f (x, y) = 0 under normal incidence.The cylinder axis (z axis) is perpendicular to the paper plane. . . . . . . . .

Fig 3 . 2

 32 Fig 3.2 Local coordinate system (n, τ) to describe the wave vectors at one interaction.

Fig 3 . 3 Fig 3 . 4 Fig 3 . 5

 333435 Fig 3.3 Notations for the order, the exit point and the direction of a scattered ray. . Fig 3.4 Schematic diagram of the incremental variation of wavefront in homogeneous medium. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fig 3.5 Variation of the wavefront of a ray pencil refracted by a circular cylinder. .

Fig 3 . 6 Fig 3 . 7 Fig 3 . 8 Fig 3 . 9 Fig 3 . 10

 36373839310 Fig 3.6 Schematic diagram in counting the optical path of a p = 1 ray relative to the reference ray. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fig 3.7 The focal line for a converging wave of cylindrical wavefront. . . . . . . . Fig 3.8 Geometry of the composite elliptical cylinder (CEC). . . . . . . . . . . . . Fig 3.9 Comparison of VCRM with LMT, DSE and GOA for the scattered intensity of plane wave by a circular cylinder. The curves for DSE, LMT and GOA have been offset respectively by 10 -4 , 10 -2 and 10 2 for clarity. . . . . . . . Fig 3.10 Comparison of the calculated scattering diagram of an elliptical cylinder (b 1 = b 2 = b = 50 µm, a = 1.2b) with that of a long ellipsoid (a = 60 µm, b = 50 µm, c = 5000 µm). For clarity, diffraction is not included. . . . . . Fig 3.11 The variation of the axis ratios ξ i (b i /a) with the equivalent radius r 0 . . . . Fig 3.12 The scattered intensity of plane wave by five CECs of diverse deformations.

Fig 3 .

 3 Fig 3.13 The scattering diagrams of a CEC as the temperature T goes up from 5 • C to 30 • C. The cross section of the CEC takes the shape of a raindrop of equivalent radius 0.75 mm. For clarity, the curves for T = 10 • C, 15 • C,...and 30 • C have been shifted by 10 2 , 10 4 ,...and 10 10 , respectively. . . . . . . . .

Fig 3 . 14

 314 Fig 3.14 Detail for the variation of the first-order rainbow as the temperature goes up.

Fig 3 . 15 Fig 3 . 16

 315316 Fig 3.15 Variation of the scattering diagram of a CEC as the incident angle θ 0 varies from -20 • to 20 • in steps of 10 • . The cross section of the CEC is the same as the one used in Fig. 3.13. m = 1.333. . . . . . . . . . . . . . . . . . . . Fig 3.16 The scattering angles of the first-order, the second-order and the fifth-order rainbows for different incident angles. (a) a circular cylinder of radius 0.75 mm; (b) a deformed CEC of equivalent radius 0.75 mm (the one used in Fig. 3.15). m = 1.333. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Fig 4 . 1 Fig 4 . 2

 4142 Fig 4.1 Local coordinate system (τ, n, ê⊥ ) for describing the incident wave vector ⃗ k i , the reflected ⃗ k r and the refracted ⃗ k t . . . . . . . . . . . . . . . . . . . . . Fig 4.2 Schematic diagram of a local wavefront. . . . . . . . . . . . . . . . . . . .

Fig 4 . 3

 43 Fig 4.3 Evolution of a local wavefront within a homogeneous medium. . . . . . . .

Fig 4 . 4

 44 Fig 4.4 Schematic diagram of the sudden change of wavefront caused by refraction.

Fig 4 . 5 Fig 4 . 6

 4546 Fig 4.5 Illustration of the two coordinate systems (ê 1 , ê2 ) and (ê (i) ∥ , ê⊥ ) within the vibration plane of incident electric vector. . . . . . . . . . . . . . . . . . . Fig 4.6 Definition of the base vectors (ê (i) ∥ , ê⊥ ), (ê (r) ∥ , ê⊥ ) and (ê (t) ∥ , ê⊥ ) for describing the parallel and the perpendicular components of the incident, reflected and refracted electric vectors. . . . . . . . . . . . . . . . . . . . . . . . . . . .

Fig 4 . 7

 47 Fig 4.7 Schematic diagram in calculating the optical path of a scattered ray relative to the reference ray. For clarity, an example for a scattered ray of order p = 1 is shown. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Fig 4 . 8

 48 Fig 4.8 Definition of the direction (φ, ψ) of a scattered ray in the 3D space. . . . . .

Fig 4 . 9 Fig 4 . 10

 49410 Fig 4.9 Schematic diagram of the 3D scattering space discretized uniformly with grid points. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fig 4.10 Schematic diagram of two groups of scattered rays (marked as the sample points) which are folded with each other. . . . . . . . . . . . . . . . . . .

Fig 4 . 11

 411 Fig 4.11 Procedure of the triangulation. (1) In each iteration, only two adjacent layers of incident rays are considered; (2) The scattered rays are irregularly distributed in terms of the scattering direction (φ, ψ); (3) The triangles meshed in the triangulation of the sample points. . . . . . . . . . . . . . . . . . . .

Fig 4 .

 4 Fig 4.12 The Delaunay criterion used in triangulation. (a) A i+1 is not a successful candidate since the circumcircle of triangle △A i B j A i+1 contains other sample points; (b) B j+1 with A i B j forms successfully a triangle. . . . . . . . . . . .

Fig 4 . 13

 413 Fig 4.13 Inclusion of a grid point inside a meshed triangle. The amplitude and the phase at the enclosed grid point, referred to as the query point (φ q , ψ q ), are to be interpolated. In the right illustrates this query scattering direction and the three sampled rays around it. . . . . . . . . . . . . . . . . . . . . . . .

Fig 4 . 14 Fig 5 . 1

 41451 Fig 4.14 Schematic diagram of two intersecting triangles. A grid point inside the intersecting area is an interference point, where the interference of light must be taken into account. . . . . . . . . . . . . . . . . . . . . . . . . . .

Fig 5 . 2

 52 Fig 5.2 The stream-wise curvature of jet surface as function of the falling height h.

Fig 5 . 3

 53 Fig 5.3 Schematic diagram of the ray path for a light ray interacting with a liquid jet.

Fig 5 . 4

 54 Fig 5.4 Experimental observation of the scattering field near the first-, the second-, the fifth-and the sixth-order rainbows by a capillary water jet. . . . . . . .

Fig 5 . 5

 55 Fig 5.5 Tracing of the scattered rays: (a) by an infinite circular cylinder; (b) by a real jet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Fig 5 . 6

 56 Fig 5.6 The configuration and the corresponding 3D scattered intensity by the liquid jet. (a) the beam profile incident on the jet; (b) the scattering field near the first-and second-order rainbows. For clarity, the intensity of p = 0 and p = 3, I 0 and I 3 , have been amplified by a factor of 10. . . . . . . . . . . .

Fig 5 . 7 86 Fig 5 . 8 87 Fig 5 . 9 88 Fig 5 . 10 89 Fig 5 . 11 91 Fig 6 . 1 94 Fig 6 . 2 96 Fig 6 . 3 99 Fig 6 . 4 99 Fig 6 . 5 101 Fig 6 . 6 103 Fig 6 . 7 104 XIFig 6 . 8 Fig 6 . 9 Fig 6 . 10 Fig 6 .Fig 6 .Fig 6 . 14 Fig 6 . 16

 57865887598851089511916194629663996499651016610367104686961066614616 Fig 5.7 The scattered intensity near the first-and the second-order rainbows (φ = [125.0 • , 140.5 • ]) for an infinite cylinder and a real jet. . . . . . . . . . . . 86

Fig 7 . 2 AFig 7 . 3

 7273 Fig 7.2 A wavy sea surface when the wind speed is 4.4 m/s. . . . . . . . . . . . .Fig 7.3 Preliminary result for the temporal-spatial distribution of the light rays transmitted from the sea at the height of a satellite (485 km). The colormap illustrates the received photon numbers (relative values). The wind speeds for the left column and the right column are 4.4 m/s and 12.3 m/s, respectively.

γ

  : divergence angle of incident beam ε : variation of amplitude due to reflection and refraction coefficients ζ : curvature of wavefront η : distance factor θ : scattering angle in 2D scattering Θ : projection matrix λ : wavelength of light ρ : curvature radius of particle surface Σ inc : incident plane τ : unit vector tangent to the particle surface ϕ : phase φ : azimuth angle in 3D scattering ψ : elevation angle in 3D scattering Subscripts or superscripts X : ⊥ or ∥ respectively for perpendicular or parallel polarization ⊥ : the component of electric field perpendicular to the incident plane ∥ : the component of electric field parallel to the incident plane p : for p-order scattered ray (i) : parameters for incident ray (r) : parameters for reflected ray (t) : parameters for refracted ray n J r , J t : Jones matrices respectively for reflection and refraction k : wave number ⃗ k : wave vector k n : normal component of wave vector k τ : tangent component of wave vector m : relative refractive index of particle M t : refractive index of the refracted medium relative to the incident medium M : coordinate transformation matrix n : unit vector normal to the particle surface O u , O v : positions of the beam waists respectively in ( ŵ, û) and ( ŵ, v) planes p : the order of scattered ray Q : curvature matrix of wavefront r : geometric length from the exit point to the observation point r 0u , r 0v : beam-waist radii respectively in ( ŵ, û) and ( ŵ, v) planes r 0y , r 0z : beam-waist radii respectively in ( x, ŷ) and ( x, ẑ) planes r ⊥ , r ∥ : Fresnel reflection coefficients

  extended the GLMT to the scattering of Gassian beam by infinite circular cylinders, while Gouesbet and Méès contributed to the extension of GLMT to the scattering by infinite elliptical cylinders [72, 73]. Han and Wu et al. devoted to the scattering of shaped beams by spheroidal particles [74-76]. Wang and Han et al. applied the GLMT to the scattering of shaped beam by a sphere with an eccentrically located spherical inclusion [77, 78]. Later, Wang et al. proposed an algorithm in the framework of GLMT to address the internal field distribution of a

  Ren et al. applied the VCRM to the light scattering in the transversal plane of an ellipsoidal particle in 2012 [139]. Jiang et al. applied it to the scattering of plane wave and Gaussian beam by elliptical cylinders in 2012 and 2013 [140, 141]. Onofri et al. retrieved the evolution of the principal curvature radii and the refractive index of an oblate droplet with a minimization method that involves VCRM predictions and experimental light scattering patterns in 2015 [138]. Sun et al. improved the three-dimensional (3D) ray tracing by considering the wave-front distortion and phase shift in the scattering by a spheroid in 2016 [142].

Fig 2 . 1 Fig 2 . 2

 2122 Fig 2.1 Schematic diagram of an infinite circular cylinder illuminated by a plane wave.

Fig 2 . 3

 23 Fig 2.3 Schematic diagram of the divergence of a ray pencil interacting with a circular cylinder.

4

 4 Fig 2.4 Schematic diagram of the virtual rays (dotted) serving as the references in calculating the optical paths of scattered rays.

)Fig 2 . 5

 25 Fig 2.5 Two types of focal lines defined by van de Hulst[1] in the scattering by a sphere. For a circular cylinder concerned in this section, only the focal lines of type A exist.

Fig. 2 .

 2 7 shows an example for decomposing the total interference field into several single-order scattered light.

Fig 2 . 6 Fig 2 . 7

 2627 Fig 2.6 Comparison of the scattered intensity calculated by GOA with that by LMT for different size parameters. m = 1.3322.

Fig 2 . 8

 28 Fig 2.8 First-order rainbow calculated by GOA for an infinite circular cylinder. The calculation parameters are the same with those in Fig. 2.7(a).

Fig. 2 .Fig 2 . 9

 229 Fig. 2.9(a).But for a spherical particle of a larger radius (a = 300 µm) as presented in Fig.2.10(b), the surface waves of orders p = 1, 2,... and 6 have much lower intensity in most of the affected regions. On the other hand, according to the amplitude formula given in Eq. (2-30),

Fig 2 .

 2 Fig 2.11 Zoomed view of Fig. 2.9(c) in linear scale for the first-order rainbow calculated by GOA and LMT.

Fig 3 . 1

 31 Fig 3.1 An infinite cylinder of cross section f (x, y) = 0 under normal incidence. The cylinder axis (z axis) is perpendicular to the paper plane.

2

 2 Fig 3.2 Local coordinate system (n, τ) to describe the wave vectors at one interaction.

Fig. 3 .Fig 3 . 4 Fig 3 . 5

 33435 Fig 3.4 Schematic diagram of the incremental variation of wavefront in homogeneous medium.

Fig 3 . 6

 36 Fig 3.6 Schematic diagram in counting the optical path of a p = 1 ray relative to the reference ray.

Fig 3 . 7

 37 Fig 3.7 The focal line for a converging wave of cylindrical wavefront.

Fig 3 . 8

 38 Fig 3.8 Geometry of the composite elliptical cylinder (CEC).

( b 2 k 2 x + a 2 k 2 y) ( 3 - 33 )

 22333 ). The two successive interactions are on the same (upper or lower) part. (x p , y p ) and (x p+1 , y p+1 ) satisfy the same equation, Eq. (3-30) or Eq. (3-31), and η is solved as η = -2 ( b 2 k x x p + a 2 k y y p ) / where b = b 1 if they are on the upper part, otherwise b = b 2 .

Fig 3 . 9 Fig 3 . 10

 39310 Fig 3.9 Comparison of VCRM with LMT, DSE and GOA for the scattered intensity of plane wave by a circular cylinder. The curves for DSE, LMT and GOA have been offset respectively by 10 -4 , 10 -2 and 10 2 for clarity.

Fig 3 . 11

 311 Fig 3.11 The variation of the axis ratios ξ i (b i /a) with the equivalent radius r 0 .

Fig 3 . 12

 312 Fig 3.12 The scattered intensity of plane wave by five CECs of diverse deformations.

Fig 3 . 13

 313 Fig 3.13 The scattering diagrams of a CEC as the temperature T goes up from 5 • C to 30 • C. The cross section of the CEC takes the shape of a raindrop of equivalent radius 0.75 mm. For clarity, the curves for T = 10 • C, 15 • C,...and 30 • C have been shifted by 10 2 , 10 4 ,...and 10 10 , respectively.

Fig 3 . 14

 314 Fig 3.14 Detail for the variation of the first-order rainbow as the temperature goes up.

1 bFig 3 . 15

 1315 Fig 3.15 Variation of the scattering diagram of a CEC as the incident angle θ 0 varies from -20 • to 20 • in steps of 10 • . The cross section of the CEC is the same as the one used in Fig. 3.13. m = 1.333.

Fig. 3 .

 3 Fig.3.[START_REF] Arnott | Unfolded optical glory of spheroids: backscattering of laser light from freely rising spheroidal air bubbles in water[END_REF] shows the detail for the scattering angles of the first-order, the second-order and the fifth-order rainbows. The incident angle θ 0 varies from -20 • to 20 • in steps of 2 • . For comparison, the rainbow angles of a circular cylinder are also presented. The circular cylinder has a radius of 0.75 mm (the equivalent radius of the CEC), while the other parameters are same with the CEC.

A le x a n d e r 's d a r k b a nFig 3 . 16

 316 Fig 3.16 The scattering angles of the first-order, the second-order and the fifth-order rainbows for different incident angles. (a) a circular cylinder of radius 0.75 mm; (b) a deformed CEC of equivalent radius 0.75 mm (the one used in Fig. 3.15). m = 1.333.

Fig 4 . 1

 41 Fig 4.1 Local coordinate system (τ, n, ê⊥ ) for describing the incident wave vector ⃗ k i , the reflected ⃗ k r and the refracted ⃗ k t .

  To characterize the divergence or convergence of wavefront, the local curvature of wavefront is integrated as a new property of a light ray in VCRM. In a general case, a local wavefront can be any shape, dependent on the values of its two principal curvatures ζ 1 and ζ 2 . A planar wavefront (ζ 1 = 0 and ζ 2 = 0), a cylindrical wavefront (ζ 1 = 0 or ζ 2 = 0) as discussed in Chapter 3 and a spherical wavefront (ζ 1 = ζ 2 0) are among the special cases.

Fig. 4 .Fig 4 . 2

 442 Fig. 4.2 illustrates a local wavefront whose principal curvature radii are R 1 and R 2 (reciprocals of the principal curvatures, R 1 = 1/ζ 1 and R 2 = 1/ζ 2 ). The two tangent directions û1 and û2 , along which the wavefront has the minimum and the maximum curvatures (principal curvatures ζ 1 and ζ 2 ), are called the principal directions of wavefront. The direction of the representing light ray is denoted as ⃗ k (wave vector). The two principal directions û1 and

Fig 4 . 3

 43 Fig 4.3 Evolution of a local wavefront within a homogeneous medium.

Fig. 4 .

 4 Fig. 4.3 shows an example, in the framework of VCRM, for a light ray of converging wavefront propagating in a homogeneous medium. At O, the ray has a converging wavefront

Fig. 4 .

 4 4 illustrates the sudden change of wavefront at a refraction by the interface of particle.

Incident

Fig 4 . 4

 44 Fig 4.4 Schematic diagram of the sudden change of wavefront caused by refraction.

Fig 4 . 5

 45 Fig 4.5 Illustration of the two coordinate systems (ê 1 , ê2 ) and (ê (i)∥ , ê⊥ ) within the vibration plane of incident electric vector.

Fig 4 . 6

 46 Fig 4.6 Definition of the base vectors (ê (i)∥ , ê⊥ ), (ê (r) ∥ , ê⊥ ) and (ê (t) ∥ , ê⊥ ) for describing the parallel and the perpendicular components of the incident, reflected and refracted electric vectors.

   holds too. As shown in Fig.4.7, the reference ray arrives at an artificially defined point O and then exits in the direction of the concerned scattered ray. For simplicity, the position of O is better set at the center of particle.The method proposed in Eq. for calculating the phase due to optical path is very easy

Fig 4 . 7

 47 Fig 4.7 Schematic diagram in calculating the optical path of a scattered ray relative to the reference ray. For clarity, an example for a scattered ray of order p = 1 is shown.

Fig 4 . 8

 48 Fig 4.8 Definition of the direction (φ, ψ) of a scattered ray in the 3D space.

Fig 4 . 11

 411 Fig 4.11 Procedure of the triangulation. (1) In each iteration, only two adjacent layers of incident rays are considered; (2) The scattered rays are irregularly distributed in terms of the scattering direction (φ, ψ); (3) The triangles meshed in the triangulation of the sample points.

Fig 4 .

 4 Fig 4.12 The Delaunay criterion used in triangulation. (a) A i+1 is not a successful candidate since the circumcircle of triangle △A i B j A i+1 contains other sample points; (b) B j+1 with A i B j forms successfully a triangle.

Fig 4 . 14

 414 Fig 4.14 Schematic diagram of two intersecting triangles. A grid point inside the intersecting area is an interference point, where the interference of light must be taken into account.

3 . 4 .

 34 In each iteration, a new set of rays is added from the next layer to obtain a new set of scattered rays. For two adjacent layers of incident rays, say sets A 0 and B 0 , we have two sets of scattered rays, say A and B. They are irregularly distributed in terms of scattering direction (φ, ψ) as shown in Fig. 4.11(2).

Fig 5 . 1

 51 Fig 5.1 The image and the extracted profile of the water jet near the orifice.

Fig 5 . 2 Fig 5 . 3

 5253 Fig 5.2 The stream-wise curvature of jet surface as function of the falling height h.

  Azimuth angle [deg.]

Fig 5 . 4

 54 Fig 5.4 Experimental observation of the scattering field near the first-, the second-, the fifth-and the sixth-order rainbows by a capillary water jet.

Fig 5 . 5

 55 Fig 5.5 Tracing of the scattered rays: (a) by an infinite circular cylinder; (b) by a real jet.

Fig. 5 .

 5 Fig.5.6(a) presents the configuration in calculation. A monochrome plane wave of wavelength λ = 632.8 nm is used. The direction of the incident rays is perpendicular to the jet axis, and it is set as the x axis of Cartesian coordinate system, while the incident electric vector is linearly polarized along the jet axis (z axis). The incident light has a planar wavefront of uniform intensity, i.e. a plane wave (see later for the experimental realization of this beam). The incident light spot has a thickness of 1.0 mm along the z axis, and its center on the jet is set as h = 2.4 mm. The relative refractive index of the jet m is set as 1.3322.

Fig. 5 .Fig 5 . 6

 556 Fig. 5.6(b) shows the simulated result of the 3D scattered intensity by the jet. φ and ψ are the azimuth angle and the elevation angle, respectively. For clarity, only the scattered light of p = 0, 2 and 3 near the first-and second-order rainbows (φ = [125 • , 140.5 • ]) are shown. The first-order rainbow arises from the interference of the p = 2 scattered light itself, while the second-order rainbow from the interference of the p = 3 scattered light. Due to the cross-polarization effect, the scattered light generally has a polarization state different from the incident light. And the scattered intensity are calculated according to Eq. (4-36). The hot map shows clearly the variation of intensity with the azimuth angle φ and the elevation angle ψ. Considering that the intensity of the p = 0 and p = 3 light are much weaker than that of the p = 2 light in this area, their intensity I 0 and I 3 are multiplied by a factor of 10 in calculating the hot map. We note firstly that the scattered light of different orders by the real jet are deviated out of the horizontal plane (ψ = 0 • ). And, the elevation angles for the scattered light of different orders are different. These p = 2 and p = 3 scattered light are located at positive elevation angles (ψ p=2,3 > 0 and ψ p=3 > ψ p=2 ), while p = 0 scattered light lies in the region below the horizontal plane where ψ p=0 < 0. As a result, the scattered rays of different orders do not interfere with each other. But, the scattered rays of same order p might superpose with each other and form the interference fringes as shown in Fig. 5.6(b).

Fig 5 . 7

 57 Fig 5.7 The scattered intensity near the first-and the second-order rainbows (φ = [125.0 • , 140.5 • ]) for an infinite cylinder and a real jet.

Fig 5 . 8

 58 Fig 5.8 The 3D scattered intensity by the water jet when h = 3.0 mm. The other parameters are the same with those in Fig. 5.7(b).

Fig 5 . 9

 59 Fig 5.9 Experimental setup for measuring the scattered intensity by a liquid jet.

Fig 5 . 10

 510 Fig 5.10 Schematic diagram of the experimental setup.

Fig 5 . 11

 511 Fig 5.11 Experimental observation of the scattered intensity of plane wave by the water jet near the first-order and second-order rainbows. For clarity, the intensity of the p = 0 and p = 3, I 0 and I 3 , have been amplified by a factor of 10. h = 2.4 mm.

3 •

 3 and 2.7 • . In a slice of fixed elevation angle ψ, the scattered intensity varies only with the azimuth angle φ. For the second-order rainbow, the examination is made in the slices of ψ = 3.0 • and 4.0 • . Besides, the p = 0 scattered light is recorded with the same exposure time as the second-order rainbow. As shown in Fig.5.12b, the examination of the p = 0 scattered intensity is made in two slices of two fixed azimuth angles (φ = 125.7 • and 128.7 • ), since the p = 0 intensity in this area varies mainly with the variable ψ.

  Fig 5.12 Experimental examination of the simulation results for the 3D scattered intensity near the firstorder and second-order rainbows of a real jet. The hot maps present the relative intensity in the images captured by experiment.

Fig 6 . 2

 62 Fig 6.2 Variation of the intensity and the iso-phase surface of an elliptical Gaussian beam: (a) the intensity profiles at different axial positions; (b) three iso-phase surfaces sampled respectively at w = -125 µm, w = -25 µm and w = 75 µm. The vectors normal to the iso-phase surfaces depict the propagation directions.

Fig 6 . 3

 63 Fig 6.3 The calculated principal directions for the saddle wavefront in Fig. 6.2(b).

Fig 6 . 4

 64 Fig 6.4 Two coordinate systems established respectively for the particle and the beam.

Fig. 6 .

 6 5 illustrates an example for calculating the optical path of a p = 1 ray relative to the reference ray. p p p p p p p p p p p p p p

WFig 6 . 5

 65 Fig 6.5 Schematic diagram for calculating the optical path of a scattered ray relative to the reference ray (dotted blue line) in the 3D scattering of a shaped beam.

Fig 6 . 6

 66 Fig. 6.7 shows the comparison results for the scattered intensity by the dielectric spherical particle of different refractive indices. Due to symmetry, the scattering diagrams in any two scattering planes (the planes containing the beam axis ŵ = x) are the same. Thus,

1 Fig 6 . 7

 167 Fig 6.7 Comparison of VCRM with GLMT for the scattered intensity of circular Gaussian beam by a spherical particle of different refractive indices m. For clarity, the intensity by VCRM has been shifted by a factor of 10.

Fig. 6 .Fig 6 . 8

 668 Fig.6.8(b) shows the simulated result for the scattered intensity of the elliptical Gaus-

6 .8 and 6 . 9 .

 669 The beam-waist center in the xz plane, O z , is at x = -20 mm, while the beam-waist center in the xy plane, O y , is set at x = 0 (coincides with the jet axis). The intensity profiles of this elliptical Gaussian beam at different transversal sections are presented in Fig. 6.11(b). They are sampled respectively at x = -20 mm (beam-waist center in xz plane), x = -10 mm and x = 0 (the position of the jet).

Table 6 . 5 Fig 6 . 11

 65611 Fig 6.11 An elliptical Gaussian beam with a considerable divergence angle in the vertical xz plane: (a) the beam diameter in the vertical xz plane; (b) the intensity profiles in the transversal yz planes sampled at three positions.

Fig 6 .

 6 Fig 6.12 Single-order scattered intensity of the tightly focused elliptical Gaussian beam modeled in Fig. 6.11 by the liquid jet. For clarity, I 0 and I 3 have been amplified by a factor of 10.

Fig. 6 .Fig 6 .

 66 Fig 6.13 Simulated 3D scattered intensity for two elliptical Gaussian beams of different divergence angles in xz plane.

From Fig. 6 .

 6 Polarization r 0z /mm r 0y /mm O z /mm O y /mm λ/nm 14, one can see that the elevation angle ψ of the p = 0 scattered light ranges from -3.8 • to 2.0 • , while the elevation angle of the p = 2 scattered light ranges from about -1.5 • to 4.4 • . Along the sampled line (A) where φ = 139.0 • , the intensity distribution of the p = 0 scattered light is compared with that of the p = 2 light in Fig. 6.15. In the overlapping region (-1.5 • ≤ ψ ≤ 2.0 • ), the scattered light of p = 2 and p = 0 superpose and interfere with each other. Consequently, this first-order rainbow shows different interference patterns in the two regions: ψ > 2.0 • and ψ ≤ 2.0 • .

Fig 6 . 14

 614 Fig 6.14 Zoomed view of the interference fringes near the first rainbow.

Fig 6 .

 6 Fig 6.15 Single-order intensity of the p = 0 and p = 2 scattered light along the sampled line (A).

Fig 7 . 1

 71 Fig 7.1 The shapes of natural raindrops with increasing radii (from [6]).

4 m

 4 /s (class-3 wind), the power spectral density is mainly concentrated upon the angular frequencies ranging from 1.41 to 9.39 rad/s. The corresponding wavelength of the wavy sea surface varies approximately from 0.70 m to 30.97 m. Since the roughness of sea surface is generally much larger than the wavelength of the commonly used 532 nm green light, it is reasonable to describe the incident laser beam by a bundle of light rays.

Fig 7 . 2 A

 72 Fig 7.2 A wavy sea surface when the wind speed is 4.4 m/s.

Fig. 7 .

 7 Fig. 7.3 shows the preliminary results for the temporal-spatial distribution of the light rays transmitted from the sea, observed at the height of a satellite (the turbulence of air is not included). The spot of the incident beam on the sea surface is of 2 × 2 m 2 . The refractive index of the sea water relative to the air is 1.333 and the wavelength of light in air is 532 nm.For comparison, the right column of Fig.7.3 presents the results when wind speed reaches 12.3 m/s (class-6 wind). One can see that the wavy sea surface has a significant influence on the spatial distribution of the transmitted light spot; and the spatial distribution varies with the time t. As the wind speed is increased, the transmitted light spot has a larger deviation from the center of receiving region, as the pentagram shows. It indicates that the acceptable probability of optical signal drops down when wind speed is increased.

Fig 7 . 3

 73 Fig 7.3 Preliminary result for the temporal-spatial distribution of the light rays transmitted from the sea at the height of a satellite (485 km). The colormap illustrates the received photon numbers (relative values). The wind speeds for the left column and the right column are 4.4 m/s and 12.3 m/s, respectively.

Table 3 .

 3 1 Geometric parameters of the CEC at different equivalent radius r 0 (unit: mm) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 3 .

 3 2 Refractive indices of pure water for λ = 0.6328µm at different temperatures T ( • C) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 6 .

 6 1 Parameters of the elliptical Gaussian beam in Fig. 6.2. . . . . . . . . . . .

Table 6 .

 6 2 Parameters of the circular Gaussian beam in Fig. 6.7. . . . . . . . . . . . .

Table 6 .

 6 3 Parameters of the elliptical Gaussian beam in Fig. 6.8. . . . . . . . . . . .

Table 6 .

 6 4 Parameters of the circular Gaussian beam in Fig. 6.10. . . . . . . . . . . .

	Table

Table 3 .

 3 1 Geometric parameters of the CEC at different equivalent radius r 0 (unit: mm)

	r 0 0.50 0.75 1.00	1.25 1.50
	a 0.503 0.760 1.026 1.301 1.585
	b 1 0.499 0.748 0.995 1.238 1.477
	b 2 0.490 0.711 0.905 1.070 1.209

Table 3 .

 3 2 Refractive indices of pure water for λ = 0.6328µm at different temperatures T ( • C)

	T 5	10	15	20	25	30
	m 1.3330 1.3328 1.3325 1.3321 1.3316 1.3311

Table 6 .

 6 2 Parameters of the circular Gaussian beam in Fig. 6.7.

	x	ẑ	10	10	0	0	632.8

Beam axis Polarization r 0z /µm r 0y /µm O z /µm O y /µm λ/nm

Table 6 .

 6 3 Parameters of the elliptical Gaussian beam in Fig. 6.8. Beam axis Polarization r 0z /mm r 0y /mm O z /mm O y /mm λ/nm

	x	ẑ	0.4	1.6	-300	0	632.8

Size parameter depicts the size of particle compared to the wavelength of light. It is calculated as

2πa/λ with a being the equivalent radius of particle and λ being the wavelength of light.

By smooth surface, it means that the particle surface is differentiable to the second order or more.

For example, the surface function (x/a)

+ (y/b) 2 + (z/c) 2 = 1 for an ellipsoid can be expressed by two parameters α and β as X(α, β) = (a cos α cos β, b cos α sin β, c sin α), with -π/2 ≤ α ≤ π/2 and 0 ≤ β < 2π.

In the thesis, a bold variable indicates a 2 × 2 matrix, while a variable with arrowhead means a 2D or

3D vector.

A scattered wave in the far field approximates a spherical wave.

https://www.thorlabs.us/thorproduct.cfm?partnumber=ID20

http://meyerinst.com/digital-cameras/hamamatsu/C9100-02.pdf

A circular Gaussian beam is the special case of elliptical Gaussian beam.

The ( ŵ, û) plane is expanded by base vectors ŵ and û, while the ( ŵ, v) plane is expanded by ŵ and v.

Analogous to the TE incidence in the scattering by a cylinder, here the polarization of incident beam is set along the jet axis (z axis).

I want to acknowledge the financial support from the China Scholarship Council (CSC) for my study in France as a joint PhD. student.

the x axis, while ψ is the elevation angle relative to the horizontal plane. They are calculated according to

where k x , k y and k z are x-, y-and z-component of the wave vector ⃗ k s , respectively. The directions of scattered rays are described in the (φ, ψ) coordinate system as sample points, which are always distributed irregularly (as illustrated in Fig. 4.10). The light rays from the same scattering order or from different orders may be folded with each other. In the folded area, the superposition of light occurs. The calculation for the scattered intensity at each grid point requires the interpolation of the sample points around it. It is worth noting

For a meshed triangle, its three vertexes are denoted as v j = (φ j , ψ j ) with j = 1, 2, 3.

If a grid point is inside the triangle, the coordinates of the grid point v q = (φ q , ψ q ) need to satisfy:

To accelerate the procedure of finding the grid points that are inside a concerned triangle, only the grid points which lie within the rectangular area bounded by the three vertexes are to be searched.

For a grid point v q = (φ q , ψ q ) inside the concerned triangle, its amplitude and phase are to be interpolated according to the amplitudes and phases of the three sample points (the vertexes of triangle). As have been discussed in the preceding parts of this chapter, the amplitude and phase for each 3D scattered ray can be determined in VCRM. Here we note the amplitudes (or the phases) at the sample points v j as f j ( j = 1, 2, 3). Then, the amplitude (or the phase) f q at the grid point v q = (φ q , ψ q ) can be acquired through linear interpolation:

where λ 1 and λ 2 are determined by solving

The physical meaning of (4-32) and (4-33) is that the amplitude and the phase in a query scattering direction (φ q , ψ q ) can be obtained through the interpolation of the surrounding three sampled rays. Fig. 4.13 presents a schematic diagram for this.

Grid points

Query point Fig 4 .13 Inclusion of a grid point inside a meshed triangle. The amplitude and the phase at the enclosed grid point, referred to as the query point (φ q , ψ q ), are to be interpolated. In the right illustrates this query scattering direction and the three sampled rays around it.

and û, while r 0v is the beam-waist radius in the plane containing ŵ and v. The r j ( j = u, v) is the characteristic radius of beam at which the field amplitude falls to 1/e of the axial value, and it is given by

The coordinate system (û, v, ŵ) used in describing an elliptical Gaussian beam of beam-waist radii r 0u and r 0v . For a circular Gaussian beam, r 0u = r 0v and

The origin of beam system (û, v, ŵ) is denoted as O G : (0, 0, 0). The positions of the beam waists in ( ŵ, û) and ( ŵ, v) planes 2 are (0, 0, O u ) and (0, 0, O v ), respectively. The Rayleigh distance l j indicates the axial position where the intensity is one half of the peak intensity at the beam waist. It is related to the beam-waist radius r 0 j and the wavelength λ by

The curvature radius R j of the wavefront at (0, 0, w) is calculated as

The gradient of the phase function is defined as

where ϕ ′ u , ϕ ′ v and ϕ ′ w are the first-order derivatives of the phase function with respect to u, v and w, respectively. They are calculated as:

The dependence of ϕ ′ u , ϕ ′ v and ϕ ′ w on the values of u, v and w indicates that the gradient (ϕ ′ u , ϕ ′ v , ϕ ′ w ) at different points of an iso-phase surface (a wavefront) is varied. At a point of coordinates (u, v, w), the unit vector which defines the direction of the local wave vector is

k is the normalized wave vector normal to the local wavefront and characterizes the propagation direction of the wave at (u, v, w). The negative sign indicates that the wave direction is opposite to the gradient of phase function, considering that the phase lags along the propagation direction (the phase ϕ is decreased with w as given in Eq. (6-3)). After the passage of O u , the beam becomes diverging in the ( ŵ, û) plane. Since the beam still converges in the orthogonal ( ŵ, v) plane, the shapes of the wavefronts between O u and as

Then the Gauss curvature ζ G of the local wavefront at (u, v, w) is obtained by [START_REF] Albin | Computational assessment of curvatures and principal directions of implicit surfaces from 3D scalar data[END_REF] ζ

while the mean curvature ζ M by

where

, the two principal curvatures ζ 1 and ζ 2 can then be determined by

The sign function χ equals 1 or -1, the criterion for choosing which can be found in [START_REF] Albin | Computational assessment of curvatures and principal directions of implicit surfaces from 3D scalar data[END_REF].

The corresponding principal directions in Cartesian coordinate system are calculated by From Fig. 6.10(b), one can see that for this circular Gaussian beam, the rainbows have much lower intensity than those in Fig. 6.8(b). This is because the first-and secondorder rainbows are formed by the incident light adjacent to the edge of particle [START_REF] Laven | Supernumerary arcs of rainbows: Young's theory of interference[END_REF][START_REF] Han | Surface waves/geometrical rays interferences: numerical and experimental behaviour at rainbow angles[END_REF][START_REF] Adam | The mathematical physics of rainbows and glories[END_REF], where the incident circular Gaussian beam with a beam-waist radius comparable to the particle radius has only 1/e 2 of its central intensity.

Effect of divergence angle on the scattering field

The preceding section covers the results for the 3D scattered intensity of an elliptical Gaussian beam whose divergence angle is less than 0.03 • in both the vertical and the horizontal planes. For a beam of a larger divergence angle, the propagation directions of the incident rays will be distinctly different, and the scattering field might present some interesting phenomena which do not exist in the scattering of a plane wave or an elliptical Gaussian problems, which exist in all ray models of light, will make the VCRM a powerful tool in solving the light scattering by particles/targets ABSTRACT In the framework of vectorial complex ray model (VCRM), this thesis aims to solve the three-dimensional (3D) scattered intensity of plane wave or shaped beam by a large particle of any smooth surface. The main work and achievements are summarized as follows:

As the first step, the calculation method based on VCRM for the 2D scattered intensity of plane wave by a cylinder of any smooth cross section is proposed. And the proposed method is applied to solving the scattered intensity of plane wave by a composite elliptical cylinder (CEC), whose cross section can take various shapes ranging from circular, elliptical to highly-deformed. The effects of shape deformation, refractive index and incident direction on the scattering fields, especially on the rainbows, are quantitatively analyzed.

Based on VCRM, the ray tracing, the phase shifts due to focal lines and optical path, the divergence and convergence of wavefront, and the cross polarization in 3D scattering are addressed. An interpolation algorithm based on triangulation has been developed which permits to take into account the interference of 3D scattered rays, thus breaking through the bottle-neck problem for VCRM in the extension to 3D scattering.

The proposed method, which is based on VCRM while allows to calculate 3D scattering field, is applied to simulating the 3D scattered intensity of plane wave by a real liquid jet. Furthermore, taking advantage of the ability of VCRM for interpreting the scattering mechanism, a systematic analysis is made for the scattered light of different orders, in regard to their separation or interference in 3D space. An experiment is carried out to verify the proposed method for 3D scattering and to examine the simulated results.

In the framework of VCRM, a ray description method for incident elliptical Gaussian beam is proposed, thus providing one feasible way to calculate the 3D scattered intensity of elliptical or circular Gaussian beam by a large particle of any smooth surface. The calculation for the 3D far-field scattered intensity of elliptical Gaussian beam by a real liquid jet is successfully achieved. The scattering fields near the first-and second-order rainbows for incident beams of different divergence angles are investigated in 3D space. These results as well as the proposed method open a promising way to characterize finely the structure of a real liquid jet and particles of other complex surfaces.

Keywords: light scattering, vectorial complex ray model, three-dimensional scattering, large nonspherical particles, geometrical optics approximation, liquid jets