# Extension of Vectorial Complex Ray Model and its application to the metrology of nonspherical particles

Zelong Ma

► To cite this version:

Zelong Ma. Extension of Vectorial Complex Ray Model and its application to the metrology of nonspherical particles. Classical Physics [physics.class-ph]. Université de Rouen Normandie; CORIA, 2018. English. NNT: 2018NORMR053. tel-03835054

HAL Id: tel-03835054

https://hal.science/tel-03835054

Submitted on 31 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# THESIS

For the degree of Doctor of Philosophy

Speciality : Energy

Prepared at the University of Rouen

## Title of the thesis

Extension of Vectorial Complex Ray Model and its application to the metrology of non-spherical particles

Presented and defended by
Zelong MA

Thesis defended publicly on January 31, 2018
before the jury composed of

<table><tr><td>Fabrice ONOFRI</td><td>Director of research at CNRS, IUSTI, Marseille, France</td><td>Reviewer</td></tr><tr><td>Bernard POULIGNY</td><td>Director of research at CNRS, CRPP, Bordeaux, France</td><td>Reviewer</td></tr><tr><td>Denis LEBRUN</td><td>Professor at Rouen University, CORIA-UMR 6614, Rouen, France</td><td>Examiner</td></tr><tr><td>Loïc MÉÈS</td><td>Research fellow at CNRS, LMFA, Ecole centrale de Lyon, France</td><td>Examiner</td></tr><tr><td>Claude ROZÉ</td><td>Professor at Rouen University, CORIA-UMR 6614, Rouen, France</td><td>Co-supervisor</td></tr><tr><td>Kuan Fang REN</td><td>Professor at Rouen University, CORIA-UMR 6614, Rouen, France</td><td>Supervisor</td></tr></table>

Thesis directed by Kuan Fang REN and Claude ROZÉ, laboratory CORIA-UMR 6614
<figure><img src="image_2.png" /></figure>
# THESE

Pour obtenir le diplôme de doctorat

Spécialité : Energétique

Préparée au sein de l'université de Rouen

## Titre de la thèse

Extension du modèle de Tracé de Rayons Vectoriels Complexes et application à la caractérisation d'une particule non-sphérique

Présentée et soutenue par
Zelong MA

Thèse soutenue publiquement le 31 Janvier, 2018
devant le jury composé de

<table><tr><td>Fabrice ONOFRI</td><td>Directeur de recherche au CNRS, IUSTI, Marseille, France</td><td>Rapporteur</td></tr><tr><td>Bernard POULIGNY</td><td>Directeur de recherche au CNRS, CRPP, Bordeaux, France</td><td>Rapporteur</td></tr><tr><td>Denis LEBRUN</td><td>Professeur à l'Université de Rouen, CORIA-UMR 6614, Rouen, France</td><td>Examinateur</td></tr><tr><td>Loïc MÉÈS</td><td>Chargé de recherche au CNRS, LMFA, Ecole centrale de Lyon, France</td><td>Examinateur</td></tr><tr><td>Claude ROZÉ</td><td>Professeur à l'Université de Rouen, CORIA-UMR 6614, Rouen, France</td><td>Codirecteur de thèse</td></tr><tr><td>Kuan Fang REN</td><td>Professeur à l'Université de Rouen, CORIA-UMR 6614, Rouen, France</td><td>Directeur de thèse</td></tr></table>

Thèse dirigée par Kuan Fang REN et Claude ROZÉ, laboratoire CORIA-UMR 6614
<figure><img src="image_3.png" /></figure>
*My goal is simple. It is a complete understanding of the universe, why it is as it is and why it exists at all.*

Stephen Hawking
<figure><img src="image_4.png" /></figure>
Acknowledgements

I would like to pay special thankfulness, warmth and appreciation to the persons be-
low who made my research successful and assisted me at every point to cherish my goal:

Firstly, I would like to express my deep and sincere gratitude to my supervisor Professor Kuan Fang REN for giving me the opportunity to do research and providing invaluable guidance throughout of my research and writing of this thesis. He has taught me the methodology to carry out the research and to present the research works as clearly as possible. And his dynamism, vision, sincerity and motivation have deeply inspired me. Ever since, Professor REN has supported me not only by providing a research assistantship, but also academically and emotionally thorough the rough road to finish my thesis. It was a great privilege and honor to work and study under his guidance. I am extremely grateful for what he has offered me and could not have imagined having a better supervisor and mentor for my Ph.D study.

Besides, I would like to thank my co-supervisor Professor Claude ROZÉ, for his en-
couragement and supports. His creative idea and excellent insight into the physical
phenomenon bring me a open mind to do research work. The confidence and freedom
he gave me has boosted all my interests in light scattering studies.

I sincerely wish to express my special gratitude to Dr. Said Idlahcen, Engineer of Research in CORIA. He has also given me much instructions and great help in the experiment preparation and realizations. Without his precious support it would not be possible to conduct this research.

Thanks are also due to the members of the jury, Dr. Fabrice ONOFRI, Dr. Bernard POULIGNY, Professor Denis LEBRUN and Dr. Loïc MÉÈS for many useful suggestions regarding the text. I would express my respect for them for their high scientific competence, their carefulness in reading the manuscript, and their very precious comments. The manuscript has been revised according to their remarks.

Last but not the least, I would like to thank my family: my parents and my fiancée for supporting me spiritually throughout writing this thesis and my life in general.
<figure><img src="image_5.png" /></figure>
<table><tbody><tr><td>3.2.4</td><td>Numerical results and comparisons</td><td>53</td></tr><tr><td>4</td><td>Rainbow of a spheroidal particle</td><td>57</td></tr><tr><td>4.1</td><td>Revisit of Airy theory</td><td>57</td></tr><tr><td>4.1.1</td><td>Comparison of Airy theory with Debye theory</td><td>57</td></tr><tr><td>4.1.2</td><td>Comparison between Airy theory and VCRM</td><td>58</td></tr><tr><td>4.1.3</td><td>Discrepancy on the particle sizing</td><td>60</td></tr><tr><td>4.2</td><td>Rainbows of spheroidal particles</td><td>66</td></tr><tr><td>4.2.1</td><td>Amplitude of scattered field near rainbow angle</td><td>68</td></tr><tr><td>4.2.2</td><td>Peak intensity of supernumerary bows</td><td>74</td></tr><tr><td>4.2.3</td><td>Incident angles of rays arriving at bow peaks</td><td>76</td></tr><tr><td>4.2.4</td><td>Rainbow intensity ratio of a prolate spheroid</td><td>77</td></tr><tr><td>4.2.5</td><td>Hyperbolic umbilic foci of an oblate particle</td><td>83</td></tr><tr><td>4.3</td><td>Conclusion</td><td>85</td></tr><tr><td>5</td><td>Experimental study of scattering by a pendant drop</td><td>87</td></tr><tr><td>5.1</td><td>Experimental setup</td><td>87</td></tr><tr><td>5.1.1</td><td>Consideration of the measurement system</td><td>87</td></tr><tr><td>5.1.2</td><td>Composition of the measurement system</td><td>88</td></tr><tr><td>5.1.3</td><td>Acquisition and image processing software</td><td>90</td></tr><tr><td>5.2</td><td>Experimental observations</td><td>91</td></tr><tr><td>5.3</td><td>Image processing</td><td>92</td></tr><tr><td>5.4</td><td>Results and conclusions</td><td>96</td></tr><tr><td>6</td><td>Conclusions and Perspectives</td><td>99</td></tr><tr><td>6.1</td><td>Conclusions</td><td>99</td></tr><tr><td>6.2</td><td>Perspectives</td><td>101</td></tr><tr><td colspan="2">Appendices</td><td>103</td></tr><tr><td colspan="2">A Interpolation of scattered field</td><td>105</td></tr><tr><td colspan="2">B Macros of ImageJ for image processing</td><td>107</td></tr><tr><td>B.1</td><td>Macro for Batch background noise reduction</td><td>107</td></tr><tr><td>B.2</td><td>Macro for the treatment of images of drops</td><td>109</td></tr><tr><td>B.3</td><td>Macro for scattering pattern</td><td>112</td></tr><tr><td colspan="2">Bibliography</td><td>115</td></tr><tr><td colspan="2">Résumé</td><td>125</td></tr><tr><td colspan="2">Abstract</td><td>126</td></tr></tbody></table>
# List of Abbreviations and Symbols

## Abbreviations

<table><tr><td>GO</td><td>Geometrical Optics</td></tr><tr><td>HU</td><td>hyperbolic umbilic</td></tr><tr><td>LMT</td><td>Lorenz-Mie Theory</td></tr><tr><td>VCRM</td><td>Vectorial Complex Ray model</td></tr></table>

## Roman Symbols

<table><tr><td>a</td><td>radius of a spherical particle<br>semi-axis of a spheroidal particle, i.e. equatorial radius<br>semi-axis in x direction of an ellipsoidal particle</td></tr><tr><td>A</td><td>amplitude of an emergent ray at a distance r from the particle center</td></tr><tr><td>b</td><td>semi-axis in y direction of an ellipsoidal particle</td></tr><tr><td>c</td><td>semi-axis in z direction of an ellipsoidal particle</td></tr><tr><td>C</td><td>curvature matrix of a dioptric surface</td></tr><tr><td>d</td><td>distance between two successive interaction points within the particle</td></tr><tr><td>D</td><td>divergence factor of an emergent ray in VCRM</td></tr><tr><td>D<sub>s</sub></td><td>divergence factor of an emergent ray scattered by a sphere</td></tr><tr><td>ϑ</td><td>distance free divergence factor</td></tr><tr><td>E<sub>i</sub></td><td>incident electric field</td></tr><tr><td>E<sub>1</sub></td><td>electric field inside the sphere</td></tr><tr><td>E<sub>s</sub></td><td>scattered electric field</td></tr><tr><td>F</td><td>dimensionless function (intensity) of scattering light</td></tr></table>
$H_i$ incident magnetic field

$H_1$ magnetic field inside the sphere

$H_s$ scattered magnetic field

$i$ distance free intensity of scattering light

$I$ intensity of scattered light at a distance $r$ from the particle center

$I_0$ intensity of incident light

$k$ wave number

$\mathbf{k}$ wave vector of rays in the surrounding media

$\mathbf{k}'$ wave vector of rays in the particle

$m$ refractive index of particle

$\mathbf{n}$ normal vector of the dioptric surface at incident point

$p$ order of ray

$Q$ wave front curvature matrix of incident ray

$Q'$ wave front curvature matrix of reflected or refracted ray

$r$ Fresnel reflection coefficient

$R$ wave front curvature radius before interaction

$R'$ wave front curvature radius after interaction

$R''$ wave front curvature radius of reflected ray

$S$ complex amplitude of an emergent ray at the distance $r$ from the particle center

$x$ particle size parameter

## Greek Symbols

$\alpha$ incident angle of a ray relative to the normal of a dioptric surface

$\beta$ refracted angle of a ray relative to the normal of a dioptric surface

$\epsilon$ relative amplitude of ray calculated by Fresnel coefficients

$\gamma$ slope of the ray

$\tau$ unit vector tangent to the interaction surface and the incident plane

$\lambda$ wavelength of ray

$\rho$ principal curvature radius of ellipsoid
$\theta$ scatting angle

$\theta_p$ total deviation angle from the original direction

$\theta_{Airy,K}$ angular position of Airy maxima

$\theta_{rg}$ rainbow angle in the geometrical optics

$\Phi$ phase of ray

$\phi$ azimuth angle

$\kappa$ aspect ratio of spheroid

$\chi$ relative intensity ratio

Subscripts or superscripts

$X$ state of polarization, 1 for perpendicular polarization and 2 for parallel polarization

1 used in the curvature radius in the symmetric plane

2 used in the curvature radius in the direction perpendicular to the symmetric plane
List of Figures

<table>
  <tbody>
    <tr>
      <td>1.1</td>
      <td>Rainbow in nature (photographed by the author in Rouen's city center).</td>
      <td>17</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>Scattering geometry, spherical coordinate system.</td>
      <td>24</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Debye series model of scattering by a homogeneous sphere.</td>
      <td>27</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>Light paths of ray interaction with a sphere.</td>
      <td>29</td>
    </tr>
    <tr>
      <td>2.4</td>
      <td>The geometric rainbow angle for first, second and third order rainbows as function of refractive index.</td>
      <td>34</td>
    </tr>
    <tr>
      <td>2.5</td>
      <td>Coordinate system for Airy's theory of the rainbow.</td>
      <td>35</td>
    </tr>
    <tr>
      <td>2.6</td>
      <td>Intensity distribution calculated by Airy theory for a sphere of refractive index m = 1.333 and radius a = 50 µm</td>
      <td>37</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Schema of an ellipsoid illuminated by a plane wave</td>
      <td>45</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Amplitude of scattered light near the second rainbow angle for a plane wave λ = 0.6328 µm scattered by a sphere of refractive index m = 1.333 and radius a = 50 µm</td>
      <td>52</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Scattering diagram predicted by VCRM and LMT for a sphere of radius a = 50 µm illuminated by a plane wave of wavelength λ = 0.6328 µm. Refractive index of particle m = 1.333. The results of VCRM are shifted by 10<sup>2</sup> for clarity.</td>
      <td>53</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Same parameter setting as in Fig. 3.3(a) except particle radius a = 30 µm.</td>
      <td>54</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Same parameter setting as in Fig. 3.3(a) except particle radius a = 10 µm.</td>
      <td>54</td>
    </tr>
  </tbody>
</table>
3.6 Scattering diagrams computed by VCRM for a prolate spheroid ($b = 65$ µm and aspect ratio $κ = 1.3$) with that of a sphere (radius $a = 50$ µm). The plane wave of wavelength $λ = 0.6328$ µm propagates along $z$ axis polarized in $y$ axis. The observe plane is fixed in $xz$ plane. The results of VCRM are shifted by $10^2$ for clarity. 55

3.7 Scattering diagrams computed by VCRM for an oblate spheroid ($b = 45$ µm and aspect ratio $κ = 0.9$) with that of a sphere (radius $a = 50$ µm). The plane wave of wavelength $λ = 0.6328$ µm propagates along $z$ axis polarized in $y$ axis. The observe plane is fixed in $xz$ plane. The results of VCRM are shifted by $10^2$ for clarity. 56

4.1 The comparison of main peak and supernumerary bows calculated by Debye series and Airy theory for a sphere of radius $a = 50$ µm and refractive index $m = 1.333$. 58

4.2 First and second order rainbows predicted by Debye theory, VCRM and Airy theory for a sphere whose radius $a = 50$ µm and refractive index $m = 1.333$. 59

4.3 Same parameter setting as in Fig. 4.2 except particle radius $a = 100$ µm. 60

4.4 Same parameter setting as in Fig. 4.2 except particle radius $a = 500$ µm. 60

4.5 Same parameter setting as in Fig. 4.2 except refractive index $m = 1.2$. 61

4.6 Same parameter setting as in Fig. 4.2 except refractive index $m = 1.5$. 61

4.7 Angular positions of main peak and subsidiary supernumerary bows predicted by Airy theory, Debye and VCRM as function of particle radius. The rainbow angle of GO is also added for comparison. 64

4.8 Angular differences of main peak and supernumerary bows between VCRM and Debye, and those between Airy and Debye, as function of particle radius. 64

4.9 Angular differences of main peak and supernumerary bows between VCRM and Debye, and those between Airy and Debye, as function of refractive index for a droplet whose radius $a = 50$ µm. 65

4.10 Size error introduced by Airy theory and VCRM respect to Debye series, as function of particle radius. 67

4.11 Scattering angle of rays $p = 2, 3, 4, 5$ as function of incident angle, and angular positions of main supernumerary peaks for a spheroid of refractive index $m = 1.33$) and symmetric semi-axis $a = 50$ µm. 77
<table>
  <tbody>
    <tr>
      <td>4.12</td>
      <td>Scattering angle of the first rainbow (p = 2) as function of incident angle, as well as its geometric optics rainbow angle and the angular positions of the main peak for spheroid/sphere of m = 1.33, a = 20, 50, 500 and 1500 µm. The incident wavelength is 0.6328 µm.</td>
      <td>78</td>
    </tr>
    <tr>
      <td>4.13</td>
      <td>Same as Fig. 4.12 but for the second rainbow (p = 3).</td>
      <td>78</td>
    </tr>
    <tr>
      <td>4.14</td>
      <td>Same as Fig. 4.12 but for the third rainbow (p = 4).</td>
      <td>79</td>
    </tr>
    <tr>
      <td>4.15</td>
      <td>Same as Fig. 4.12 but for the second rainbow (p = 5).</td>
      <td>79</td>
    </tr>
    <tr>
      <td>4.16</td>
      <td>Comparison of scattering diagrams for a sphere and a spheroid illuminated by a plane wave of λ = 0.6328 µm in perpendicular polarization.</td>
      <td>81</td>
    </tr>
    <tr>
      <td>4.17</td>
      <td>First and second order rainbow of an spheroidal water droplet with a = 50 µm and b = 50, 55, 60, 65 µm (i.e. κ=1.0, 1.1, 1.2, and 1.3) illuminated by a plane wave of λ = 0.6328 µm in perpendicular polarization.</td>
      <td>83</td>
    </tr>
    <tr>
      <td>4.18</td>
      <td>Comparison of VCRM and experimental normalized equatorial scattering diagrams of a Di-Ethyl-Hexyl-Sebacat droplet [1]</td>
      <td>84</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Experimental setup</td>
      <td>88</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Schema of experimental set-up</td>
      <td>89</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Interface of the image acquisition and camera control software for our linear CCD</td>
      <td>91</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Images of pendant drops and their scattering pattern near rainbow angle. The region between two yellow line is the equatorial region selected. The yellow and red curve along the rainbow pattern are respectively the outline extracted and curve fitted for the main peak of rainbow.</td>
      <td>93</td>
    </tr>
    <tr>
      <td>5.5</td>
      <td>image of scattering diagram, and its intensity distribution in equatorial plane got from experiment</td>
      <td>95</td>
    </tr>
    <tr>
      <td>5.6</td>
      <td>image of scattering diagram, and its intensity distribution in a selected window (wide = 40 pixels) near equatorial plane</td>
      <td>95</td>
    </tr>
    <tr>
      <td>5.7</td>
      <td>Intensity ratio of first to second order rainbow as function of aspect ratio.</td>
      <td>97</td>
    </tr>
    <tr>
      <td>B.1</td>
      <td>Scattering pattern near the rainbow angle before and after treatment.</td>
      <td>108</td>
    </tr>
    <tr>
      <td>B.2</td>
      <td>Scattering diagram in the equatorial plane before and after treatment.</td>
      <td>108</td>
    </tr>
    <tr>
      <td>B.3</td>
      <td>The image of drop before and after making binary</td>
      <td>109</td>
    </tr>
  </tbody>
</table>
# Chapter 1

# Introduction

## 1.1 Optical metrology

In the environment control, the biochemistry, the fluid mechanics and the combustion fields, we need to measure the characteristics of particles, like their size, shape, velocity, temperature etc. Among different kinds of measurement techniques, optical metrology is largely employed because of being accurate, reliable and non-intrusive.

A large effort has been devoted to develop powerful optical techniques for non-intrusive measurements of particle dynamics and characteristics since the invention of laser. According to the measurement principle and quantities to be measured, optical particle characterization techniques can be classified as direct imaging, intensity or intensity ratio, interferometry, time shift, pulse delay and inelastic scattering [2]. Among them, the intensity-based light-scattering methods are most attractive owing to their advantages of being simple in realization and rapid in the interpretation of the results. A rather large number of measurement techniques fall under this general classification, like Laser Doppler Velocimetry (LDV), Phase Doppler Anemometry (PDA), rainbow refractometry (RRF), diffraction and holography.

Laser Doppler Velocimetry and Phase Doppler Anemometry are two classical techniques in fluid mechanics. The former uses frequency of signals acquired from small inhomogeneities in the flow, usually tracer particles, to measure the flow velocity [3]. The latter is developed as an extension of LDA for measuring size and velocity of a moving particle simultaneously [4]. Its principle of velocity measurement is the same as that used in LDA. For particle size measurement, the scattered light is collected by two or three detectors. The diameter of the particle is deduced from the phase difference between the signals. A large efforts were made to their extensions and configurations like extended PDA [5], dual mode PDA [6], dual burst PDA [7]. However the particles have to be spherical, otherwise the diameter and consequently the mass flux cannot be determined correctly.
Rainbow refractometry is a versatile intensity-based light-scattering diagnostic technique for determining the refractive index and sizes for both, single liquid droplets [8] and in entire sprays [9, 10]. The principle of RRF is that the angular position of a rainbow is very sensible to the refractive index and its supernumerary structure depends mainly on the size of the particle. By analyzing the rainbow signal, the droplet size and refractive index, can be simultaneously extracted. However in the implementation of classical rainbow refractometry, two conditions must be satisfied: (i) the droplets must be spherical, (ii) the size of the droplet must be far greater than the wavelength of incident light [11].

Similarly, for the other optical techniques, such as diffraction and holography, the shape of the particles is usually considered to be spherical [12]. Therefore the characterization of non-spherical particles is still a veritable challenge in the development of optical metrology in very large scope of applications. It is mainly due to a lack of theoretical model to describe the interaction of light with large non-spherical particles.

## 1.2 Classical theories and models for light scattering

All the optical techniques listed above are based on the achievement of fundamental researches on light scattering theory. Various theories and models have been developed to describe the interaction between light and particles. They can be generally categorized into three categories: rigorous theories, numerical methods and approximate models. In practice, we can choose an adequate method among them according to the properties of the waves and the scatters.

### 1.2.1 Rigorous theories

Light is a high frequency electromagnetic wave. Theoretically, all kinds of phenomena concerning the interaction of light and electromagnetic waves can be solved by the famous Maxwell's equations together with the associated boundary conditions. During the past decades, many researchers have devoted themselves to the development of theories by solving the Maxwell's equations with a variable separation method. And the spherical particles, because of their simplifies, are widely used for the exploration of light scattering. Among these theories, the Lorenz-Mie Theory (LMT) developed by Gustav Mie [13] and Ludvig Lorenz [14] is one of the most classical theories.

Despite Lorenz-Mie Theory was originally developed to describe the interaction between a homogeneous isotropic sphere and an electromagnetic plane wave, to meet the requirements of more practical situations, it has been generalized along two lines in the past half century: (i) particle shape and (ii) beam shape. Most particles existing
in nature or produced in industrial processes are non-spherical, therefore LMT is extended to the scattering of a plane wave by other simple shape particles, such as spheroid [15], cylinder [16–18], elliptical cylinder [19, 20], multi-layered circular cylinders [21, 22], etc. On the other hand, with the emergence of laser and new light sources, the focused laser beam is widely used in many applications such as optical manipulation, optical trapping and optical metrology. The Lorenz-Mie theory was then generalized to the application of focused laser beams. One of the most famous extensions is the Generalized Lorenz-Mie Theory (GLMT) developed by Gouesbet et al [23]. Since its foundation in 1980s, more and more researches devoted themselves to its studies. Nowadays, GLMT has been extended from a single spherical particle to multi-layered spheres [24], spheroids [25, 26], infinity cylinders [27], infinity elliptical cylinders [28], pulsed laser scattering [29] and aggregates of spheres [30].

Lorenz-Mie theory and its extensions can be used to numerically compute the scattered intensity. However none of them can shed light on the physical mechanisms that cause various prominent features appearing in the intensity. In contrast, Debye series expansion (DSE) allows a detail view of the scattering by expanding scattering coefficients in Lorenz-Mie Theory in Debye series to study the contribution of different orders of rays in rigorous regime. Until now, remarkable work has been done for cylinders [31], spheres [32], coated [33] and multilayered spheres [34], and multilayered cylinders [35]. In recent years the Debye series was developed for light scattering by spheroidal particles [36], and then generalized to non-spherical particles with complex shape using the extended boundary condition method [37].

Generally, the advantage of the Lorenz-Mie theory and its extensions is that they provide a rigorous solution to the scattering problem. As a result, they are widely referred to validate the other methods [38]. However due to the difficulties in evaluating the special functions involved in such theories, except for the spherical particle and the infinite circular cylinder at normal or oblique incidence [35], the scatterer can hardly exceed a few tens of wavelengths.

### 1.2.2 Numerical methods

The scattering problem of arbitrary shaped particle can also be dealt with by direct numerical solution of Maxwell's equations. These methods are termed as numerical methods.

Among the numerical methods, Waterman's T-matrix method [39, 40], also referred to as either the extended boundary condition method or the null field method, is an accurate and powerful tool capable of yielding a highly accurate numerical solution for the scattering of light by non-spherical particles [18, 41–43]. In addition, the elements of the T-matrix are independent of the incident and scattered fields, depending only on the scatterer and its orientation with respect to the reference frame. Therefore once the T-matrix is calculated, it can be used directly in light scatter-
ing computations for any incidence direction [44]. T-matrix method is one of the
most widely used approaches to obtain highly accurate numerical optical properties
of morphologically complex particles with moderate aspect ratios and size parameters
ranging from zero to 200 [45].

The other numerical methods can be categorized into two groups by the form of
Maxwells equations they solve: differential equation methods and integral equation
techniques. The first group uses a differential form of the governing equations and
needs the discretization in the entire domain. Two of the most famous methods in
this group are the finite element method (FEM) [46] and the finite difference time
domain method (FDTD) [47]. But both FDTD and FEM are very costly in term of
computer resources and the size parameter of the scatter is also limited.

The integral equation techniques reformulate the original boundary value problems
by applying the integral equations over the boundary interfaces or over the entire
volume. And the integral equations are converted to mathematical equations with the
help of techniques like method of moments (MoM) [48] or its modifications like the
discrete dipole approximation (DDA) [49]. MoM is widely used for scattering by large
homogeneous particles because of its accuracy and efficiency. Various acceleration
methods are developed for MoM. Among these acceleration methods, the multilevel
fast multipole algorithm (MLFMA) [50, 51] can solve challenging problems and reduce
the time and the memory requirement of MoM. In recent years, MLFMA have already
applied for computing radiation pressure force and surface stress of large non-spherical
particles whose size parameter larger than 600 [52, 53].

However, even with the development of computer technology in recent year, the
CPU time and memory resources of numerical methods remains to be very consuming
and is its main constraint on development.

1.2.3 Approximate models

As noted above, even if the rigorous theories can predict the interaction between waves
and particles, due to the difficulties in evaluating the special functions involved, except
for sphere, and infinite circular cylinder, the scatterer can hardly exceed a few tens
of wavelengths. For the numerical methods, the size parameter of the particle is also
limited and they are very time consuming and require a lot of CPU memory. On the
other hand, approximate models, even if not fully rigorous, are very flexible and can
predict the scattering of light by complex shaped particles.

For a particle whose size is much smaller than the wavelength of light, its scattering
problem can be solved by the Rayleigh approximation [54] assuming that the incident
field inside and near the particle behaves almost as an electrostatic field and the
internal field is homogeneous. For large and complex shaped particles, high frequency
approximate methods, such as the geometrical optics (GO), the physical optics, the
geometrical theory of diffraction and the physical theory of diffractions, can be applied. Geometrical optics approximation, despite some inaccuracies, is most widely applied. Its main advantage is that it can be applied to any shape. Besides its powerful computational capabilities for large particles of complex shape, GO also allows to isolate the contributions of different orders of rays and the effects of interferences, which helps us to understand the mechanism of scattering.

Studies on geometric optics have been of interest to many researchers in the past several decades. Davis [55] derived the geometric divergence factor, as well as intensity of scattered light by an air bubble in water. Van de Hulst [56] gave a detailed discussion on the GO approximation of plane wave scattering by spherical particles in his classical monograph. The intensities and phases of emergent rays are studied. Since then many researchers have devoted themselves to the improvement of geometrical optics approximation. By introducing the interference between the diffraction, the reflected and refracted rays, Glantschnig and Chen [57] further calculated the scattered intensity of water droplets in the forward angular from 0° to 60°. Ungut et al [58] gave a detailed comparison between the GO and the Mie theory in the forward scattering angles (0–20°) for transparent spherical particles in a range of diameter from 1 µm to 100 µm. Hovenac [59] used GO to predict far field light scattering by particles which are symmetric about the optical axis. The errors in the backscatter direction and the limitation of GO approximation were also analyzed. Xu et al [60] applied GO to light scattering by a coated particle and compared their results to those obtained with the Aden-Kerker theory [61], which shows that the GO is valid for the coated particles in the near forward directions. Li et al. [62] gave the GO method for the computation of light scattering intensity within a forward angular range (0 – 60°) for spheres of gradient-index. Yu and Shen [63] presented the GO method for acceleration of the computation of the scattering intensity in all directions by lager air bubbles in water. Grynko and Shkuratov [64] studied the scattering characteristics of semi-transparent particles faceted with various shapes by means of GO method. Yang et al [65, 66] combined geometrical optics with the electromagnetic wave method.

However GO is rarely applied to light scattering of a three dimensional (3D) irregular shaped particle because of the difficulties in the calculation of local divergence/convergence and the phase shift due to focal lines. It is easy to calculate the phase due to focal lines for spherical particle or circular cylinder, in which, for a given incident ray, the deviation angles of all orders can be obtained analytically and the focal lines or focal points can also be counted accordingly. In contrast, for a particle of complex form, it's impossible to obtain an analytical expression of phase due to the focal line.

In this context, a new model of light-particle interaction, called Vectorial Complex Ray Model (VCRM) has been recently developed in the laboratory. Theoretical and numerical comparisons with other methods show that VCRM is very suitable for the description of the interaction of an arbitrary wave with an object of smooth surface and complex shape [1, 67].
1.3 Rainbow and caustics of non-spherical drops

1.3.1 Rainbow of spherical particles and its applications

Our human beings have been curious about nature since time immemorial. Rainbow
phenomena, one of the most beautiful optical phenomena in nature, has long been a
source of curiosity for scientists who have tried to explain it theoretically.

We can observe rainbows in nature when we face a rain shower with the sun shining
at our backs. Fig. 1.1 is a photo of rainbow in nature. Marvelously concentric coloured
bands can be distinguished on a part or parts of a circle arranged from the inner to
the outer border as violet, indigo, blue, green, yellow, orange and red. Sometimes,
higher in the sky, a secondary rainbow is seen which is much weaker in intensity
and has its colours reversed. The first one originates from rays suffering from one
internal reflection, while the second one results from rays that have experienced two
internal reflections. With each reflection some light is lost, that is why the secondary
rainbow is fainter than the primary one. The rainbow observed in the nature is
composed of different colors, because the sun light has a continuous spectrum of
wavelengths. The dependence of the refractive index on the wavelength of the sunlight
makes the different colors at different directions, so forming the multicolored arches.
The black region between the two orders is so called Alexander's dark band, after the
Greek philosopher Alexander of Aphrodisias, who first described it in about A.D. 200.
Theoretically an infinite number of rainbows exist. These higher-order rainbows are
never seen in nature because they are weaker than the background sky brightness.

Aristotle [68] is the first who attempt to rationally explain the appearance of the
rainbow. He explained correctly the circular shape of the bow and perceived that
rainbow is the light scattered into the eyes. In 1266 the angle formed by the rainbow
rays and the incident sunlight was first measured by Roger Bacon [69]. These angles
are named as the rainbow angle in the successive researches. In 1304, Theodoric [70]
showed that each drop is individually capable of producing a rainbow. By using a
spherical cup full of water, the mechanism, as well as the angular positions of first
and second orders of rainbow are given. Theodoric's findings were ignored for three
centuries, until they were rediscovered by Descartes [71] in the seventeenth century.
Descartes gave a very satisfactory mathematical explanation of the rainbows with his
newly established geometrical optics. With the help of the geometrical optics model,
we can explain the rainbow of different orders and approximatively determine the
rainbow angle.

When a single particle is illuminated by a monochromatic light, we can observe a
set of fringes, called supernumerary bows, near the geometrical optics rainbow angle
[72]. The optical effect underlying the supernumerary arcs was firstly explained in 1803
by Thomas Young [73]. He explained the supernumerary arcs based on a wave theory
of light and showed that light is capable of interference. In 1835, Potter calculated
<figure><img src="image_6.png" /><figcaption>Figure 1.1: Rainbow in nature (photographed by the author in Rouen’s city center).</figcaption></figure>

the intensity of scattered light near the rainbow angle by introducing the calculation of caustics. Based on the model of Potter and the wave propagation principle of Hugyens-Fresnel, in 1838 Airy [74] predicted quantitatively the intensity distribution taking into account the size and the refractive index of the particle. Compared with the results predicted by Descartes and Young, the angular position of the maximum intensity of rainbow, an infinite intensity at rainbow angle and complete darkness in Alexander’s dark band are revised in Airy’s prediction.

The difference between the geometrical rainbow angle and the first Airy bow, and those between two adjacent Airy bows are calculated by Wang et al [75]. They showed that the position of Airy bows depends on the particle size and the refractive index. So in principle we can measure the refractive index and the size of the particle from the angular positions of supernumerary bows. Based on this principle, a new light-scattering technique is investigated by Roth et al [76] for determining the refractive indices of transparent, homogeneous, spherical droplets. The angles of the first intensity peak and the maxima of supernumerary bows, are measured by a linear CCD, and then used to determine the refractive index and size of the particle. This measurement technique is our aforementioned rainbow refractometry (RRF).

In the implementation of classical rainbow refractometry, the following issues should be taken into account: the size of the droplet, the temperature/refractive
index gradients inside the droplet [77], the high-frequency oscillations superimposed
on the supernumerary bows, and the sphericity of droplets. In classical rainbow refrac-
tometry, the size of the droplet should be far greater than the wavelength of incident
rays. The existence of temperature or composition gradients inside the droplets, which
induce a modification of the relationship between rainbow angular location and re-
fractive index value. The high frequency oscillations found on the supernumerary
bows, called ripple structures, are the signature of interferences between the light
directly reflected at the surface of the particle and the light internally reflected. Rip-
ple structures can increase the uncertainty of the position determination especially
when the particle size is not much greater than the wavelength. The sphericity of
the droplets may modify the angular positions of the rainbows. Van Beek and Rieth-
muller [11, 78] indicated that a shape different from the perfect sphere will introduce
the modification of the relationship between refractive index and rainbow position in
their measurement of free-falling droplets.

1.3.2 Rainbow and caustics of levitated drops

Most of natural and artificial particles have non-spherical shapes or lack a spherically
symmetric internal structure. The scattering properties of non-spherical particles can
differ dramatically from those of spheres [79]. For example, many optical phenomena
observed for non-spherical particles, like halos, arcs, pillars and zenith-enhanced lidar
backscatter observed for ice crystals, are not remarkable for a spherically symmetric
particles. Therefore, the study of light scattering by non-spherical particles is essen-
tial for developing optical diagnostic techniques to characterize deformed droplets in
terms of refractive index, size, shape, and orientation, as well as for evaluating errors
arising from the spherical shape assumption by many scattering-based measurement
techniques.

Oblate droplets, first approximation of particles deformed, exist in many industrial
processes, e.g. fuel injection for combustion, spray cooling, and spray coating. In
the experiment of light scattering by levitated drops [1], the levitated drop can be
consider as a oblate spheroid. When a plane wave is incident on the drop in the side-
on orientation, the scattered intensity observed on a distant screen exhibits a number
of prominent and interesting optical caustics. For a range of shapes, the caustic
is made of two parts, an Airy caustic (which corresponds to that of the ordinary
rainbow) and a cusp caustic spread out transverse to the generally radial direction
of the outward propagating wavefront. These caustics are full of information and
experimental studies of light scattering by a levitated drop involve relatively simpler,
cheaper instrumentation, making them very effective to study the light scattering of
oblate particles and validate theoretical model.

Caustics from a spheroidal droplet in the vicinity of the primary rainbow region
were firstly referred by Marston et al [80, 81]. They found that the modulation of the
light scattered in the rainbow region reveals small changes in the drops shape, in their
early measurements of the response of acoustically levitated liquid drops to modulated
ultrasonic radiation pressure. To distinguish such a rainbow from those arising from a
spherical droplet, the pattern from spheroidal droplets is called a generalized rainbow
pattern in the successive researches. In 1984, Marston et Trinh [82] showed how
an acoustic levitated spheroidal water drop of millimeter illuminated by a parallel
beam of light forms a caustic in the far field. By changing the ultrasonic amplitude,
they observed the generalized rainbow patterns of different axis-ratio. They found a
horizontal V-shaped focal section of a caustic is formed for a particular aspect of the
drop. This caustic is characterized by the hyperbolic umbilic (HU) catastrophe – a
generalization of the primary rainbow formed by a spherical drop. Simultaneously,
Nye [83] computed the drop’s aspect ratio for a hyperbolic umbilic focal section. In
1985, the angular position of the cusp is calculated as function of the axis ratio of
an oblate spheroid [84]. In the calculation, Marston assumes that, the diameter of
the drop illuminated should be much greater than the wavelength of incident light
and a hyperbolic umbilic caustic is limited to illumination perpendicular to the short
axis of the spheroid. Besides, the transverse cusp diffraction catastrophes [85, 86] and
E6 diffraction catastrophe [87], from scattering of levitated drops by monochromatic
and collimated white-light illumination, have also been observed and discussed. The
opening rate of the transverse cusp diffraction catastrophe in light scattered by oblate
spheroidal drops is calculated [88].

In 1992, Nye [89] studied the essential features of the far-field caustics including HU caustic, lip events [85] and E6 catastrophe according to geometrical optics and the principles of catastrophe optics. In 2000, Lock et al [90] studied the formation of the rainbow caustic, transverse cusp and HU caustics of oblate spheroids in the numerical solution of the exact electromagnetic. In 2013, the optical caustics structures, including rainbow and HU fringes in the vicinity of primary rainbow of light scattering from oblate levitated water droplets, were simulated by Yu et al [91] with their new developed vector ray tracing (VRT) model based on geometric optics. The location of cusp caustics, as well as the changes of the optical caustic structures in response to shape deformation of oblate droplets were investigated. Recently, Onofri and Ren [1] compared the scattering patterns of large oblate droplets with the predictions of Vectorial Complex Ray Model. Through the comparisons, the authors showed that the VCRM predictions, with its high computational efficiency, fit the experimental scattering patterns in the inter-caustic region very well, therefore it is suitable for the characterization of large non-spherical particles.

The rainbow and caustics of levitated drops are widely used to study light scat-
tering by non-spherical particles. However it has also drawbacks. Firstly, for large
aspect ratios (e.g. the horizontal radius / the vertical radius > 1.23), the droplet
vibrates significantly in the levitator and the generalized rainbow pattern becomes
highly instable and blurry which makes it difficult to distinguish the rainbow fringe.
For this reason, the aspect radios of the levitated droplet are limited [91]. Besides
the heating effect by the acoustic trap is also introduced and the refractive index of droplet varies during the course of the experiment. Most of all, the method of acoustic levitation can only generate oblate droplets, while the scattering properties of oblate droplets is quite different from that of prolate droplets. Therefore the study of light scattering by prolate spheroids are needed.

### 1.3.3 Rainbow pattern and caustics of pendant drops

Pendant drop is another typical non-spherical particle. On the equatorial plane it can be considered as a prolate spheroid, complement to acoustic levitated oblate spheroid. Besides, it is a natural non-spherical drop, making it easy to obtain and control. The pendant drop obtained in experiment is also very stable. Most of all, pendant drop is much rich in the scattering pattern, like twisted rainbow structure and forward rainbows and Airy like bows due to reflexion and refraction light (i.e. zero order and primary order of rainbow). Therefore its study can reveal the influences of particle non-sphericity on the rainbow patterns. Its phenomena will be explored in this thesis using VCRM.

Lower and higher-order rainbows formed by pendant drops have been experimentally studied by a number of researchers. Sassen [92] has firstly observed rainbows up to the 6th order formed by a pendant drop. From the experiment he found that the near-spherical particles can generate also primary, secondary, and associated supernumerary bows, but unexpectedly strong higher order bows as well. The scattered light is increased in certain scattering angle intervals, especially for the higher-order rainbows. Sassen concluded that it is due to the decreased curvature in the central elongated regions, thus providing relatively more surface area available within the drop for generating rainbow rays. He also suggested that large raindrops having reduced curvature at their equators would also show an increased Airy brightness. Nevertheless, the study of Langley and Marston [93] showed that vertically focused rays outside the equatorial plane may contribute significantly to rainbow brightness for smaller drops that are oblate but have their greatest curvatures at the equator. Most of all, Sassen found that, compared with the primary order rainbow, the relative strength of the secondary order rainbow increases with increasing drop size, often resulting in that the secondary order rainbow has an equal or greater strength. He supposed that it is because of the generally reduced near-backscattering efficiency of the pendant drops.

After that, Walker [72] observed rainbows up to the 17th order using a pendant drop. Lock and Woodruff [94] observed an enhancement of the glare spot, virtual images of scattered light undergone internal reflection(s) inside and emerge from the droplet, associated with the 11th-order rainbow formed by a pendant water drop. Chan and Lee [95] observed as high as 32 order rainbows formed by a pendant drop whose radius in equatorial is 2.6 mm. Ng et his coworkers [96] observed higher-order
rainbows, up to the 200th order, formed by a pendant water drop. The experimental
set up is similar to that of Chan and Lee. Deliberate misalignment of 0.1 mm from the
equatorial plane of the laser beam was made to observer the very high order rainbow.
Then Ng et al [97] studied the interference of the 11th- and higher-order rainbows
formed by a pendant water drop. Time variations of the glare-point intensity, due to
vibration of the droplet, were compared with the numerical calculations in reference
[95].

However, to our best knowledge, except the prediction of the rainbow angle for a
slightly deformed sphere by van Beek [98] and Möbius [99], there was no reliable expla-
nation for the phenomenon observed in the experiment of Sassen [92], nor analytical
expression for the intensity ratio of first to second order rainbow.

1.4 Structure of the thesis

The body of the thesis is organized as follows. In Chapter 2 classical models for light scattering by a sphere are presented. Among them, Lorenz-Mie theory (LMT) and its series version - Debye theory are rigorous solutions of the Maxwell equations and will serve us as a reference in this thesis. Geometric optics approximation is an approximate method and the base of our VCRM. Its introduction provides fundamental concepts and formula in the following chapters. Besides, the Airy theory is also introduced. We'd like to reexamine it later in the framework of VCRM to show that the Airy theory can be further improved. Its idea serves us to deal with the rainbow structure of non-spherical particles in the framework of VCRM. In chapter 3, the Vectorial Complex Ray Model (VCRM) is introduced and applied to the scattering in a symmetric plane of an ellipsoid. The scattering diagrams of spherical particles by VCRM are compared with those of LMT. The effect of particle size on the precision of VCRM is examined. In addition, the scattering diagrams of spheroidal particles are presented to show the influence of particle non-sphericity on the rainbow pattern. In chapter 4, we apply VCRM first to revisit the Airy theory and then to investigate the properties of the rainbows formed by a prolate and an oblate spheroidal particles. The positions and the intensity of different order of rainbow of spheroidal particles illuminated by a plane wave in the symmetrical plane are calculated analytically in the framework of VCRM. The peak due to the focal effect of an oblate spheroidal particle is evaluated analytically by using the divergence factor of VCRM. A rainbow measurement system for light scattering by pendant drops is set up in the laboratory. The chapter 5 is devoted to the description of the system and its exploitation. The experimental investigation of light scattering by a pendant drop will be done in terms of VCRM. The influence of particle ellipticity on intensity distribution of scattered light is studied. The last Chapter is dedicated to the conclusions and perspectives of the work.
# Chapter 2

## Classical models for light scattering
by a sphere

In this chapter, two rigorous theories and two approximate methods for light scattering by a sphere are recalled. Lorenz-Mie theory is the rigorous solution of Maxwell equations which can be used to describe the plane wave scattering by an optically homogeneous spherical particle. It is widely used as a reference to validate the other models. Debye series expansion is developed for scattering of light by a homogeneous particle which could interpret the scattered intensity in terms of various physical processes. However, these rigorous theories cannot be applied to large non-spherical particles. Geometrical optics approximation is a very simple and intuitive method. It is flexible to deal with the scattering of the plane wave by a spherical particle in the neighborhood of the rainbow angles and the scattering of non-spherical particles. Airy theory is an approximate method, widely used in the optical particle characterization.

### 2.1 Lorenz-Mie theory and Debye series expansion

The principle of Lorenz-Mie theory is that

1. Express the electric and magnetic field of the incident plane wave, the scattered wave and those in the particle in a vector spherical harmonic expansion, which satisfies Maxwell equations.

2. Apply the boundary conditions at the sphere surface to obtain the expansion coefficients.

3. Calculate the physical quantities, such as the extinction and scattering efficiencies, and scattering diagram, etc.

Consider a homogeneous, isotropic sphere of radius $a$ and refractive index $m_1$ embedded in a dielectric medium of refractive index $m_2$. The sphere is located at the
<figure><img src="image_7.png" /><figcaption>Figure 2.1: Scattering geometry, spherical coordinate system.</figcaption></figure>

center *O* of a Cartesian coordinate system as shown in Fig. 2.1 and illuminated by a plane wave polarized in *x* direction and propagating in *z* direction. The relative refractive index of the sphere to the medium *m* = *m*₁/*m*₂. The scattered light is observed at point *P* in the far field with distance *r*. The incident electric field **E**ᵢ and magnetic field **H**ᵢ can be expanded in an infinite series of vector spherical harmonics [100]

$$ \mathbf{E}_i = \sum_{n=1}^{\infty} E_n(\mathbf{M}_{o1n}^{(1)} - i\mathbf{N}_{e1n}^{(1)}) \qquad (2.1) $$

$$ \mathbf{H}_i = \frac{-k}{\omega\mu} \sum_{n=1}^{\infty} E_n(\mathbf{M}_{e1n}^{(1)} - i\mathbf{N}_{o1n}^{(1)}) \qquad (2.2) $$

The electromagnetic field inside the sphere **E**₁ and **H**₁, and the scattered field **E**s and **H**s can also be expanded in the same way:

$$ \mathbf{E}_1 = \sum_{n=1}^{\infty} E_n(c_n \mathbf{M}_{o1n}^{(1)} - i d_n \mathbf{N}_{e1n}^{(1)}) \qquad (2.3) $$

$$ \mathbf{H}_1 = \frac{-k_1}{\omega\mu_1} \sum_{n=1}^{\infty} E_n(d_n \mathbf{M}_{e1n}^{(1)} - i c_n \mathbf{N}_{o1n}^{(1)}) \qquad (2.4) $$

$$ \mathbf{E}_s = \sum_{n=1}^{\infty} E_n(i a_n \mathbf{N}_{e1n}^{(3)} - b_n \mathbf{M}_{o1n}^{(3)}) \qquad (2.5) $$

$$ \mathbf{H}_s = \frac{k}{\omega\mu} \sum_{n=1}^{\infty} E_n(i b_n \mathbf{N}_{o1n}^{(3)} + a_n \mathbf{M}_{e1n}^{(3)}) \qquad (2.6) $$
where the expression coefficient is given by

$$E_n = i^n E_0 \frac{2n + 1}{n(n + 1)} \qquad (2.7)$$

$M_{oln}$, $M_{eln}$, $N_{oln}$ and $N_{eln}$ are vector spherical harmonics. The superscripts in the functions of $M$ and $N$ denote the kind of spherical Bessel function. $a_n$ and $b_n$ in Eqs. (2.5) and (2.6) are Mie coefficients to calculate the amplitudes of scattered field, and $c_n$ and $d_n$ in Eqs. (2.3) and (2.4) are those for the internal field. These coefficients are determined by the boundary conditions of the fields at the sphere surface

$$(\mathbf{E}_i + \mathbf{E}_s - \mathbf{E}_1) \times \hat{\mathbf{e}}_r = 0 \qquad (2.8)$$

$$(\mathbf{H}_i + \mathbf{H}_s - \mathbf{H}_1) \times \hat{\mathbf{e}}_r = 0 \qquad (2.9)$$

$\hat{\mathbf{e}}_r$ is the unit vector in the radial direction. By solving the equation system generated by Eqs. (2.8) and (2.9), we can obtain the expressions of the scattering coefficients $a_n$ and $b_n$ [58]:

$$a_n = \frac{\psi_n(x)\psi_n'(\mu) - m\psi_n(\mu)\psi_n'(x)}{\xi_n(x)\psi_n'(\mu) - m\psi_n(\mu)\xi_n'(x)} \qquad (2.10)$$

$$b_n = \frac{m\psi_n(x)\psi_n'(\mu) - \psi_n(\mu)\psi_n'(x)}{m\xi_n(x)\psi_n'(\mu) - \psi_n(\mu)\xi_n'(x)} \qquad (2.11)$$

where $x = 2\pi a/\lambda$ is the size parameter, and $\lambda$ is the wavelength of plane wave. $\mu = mx$ is the second size parameter. $\psi_n(z)$ and $\xi_n(z)$ are Ricatti-Bessel functions deformed by

$$\psi_n(z) = \sqrt{\frac{\pi z}{2}} J_{n+\frac{1}{2}}(z) \qquad (2.12)$$

$$\xi_n(z) = (-1)^n \sqrt{\frac{\pi z}{2}} J_{-n-\frac{1}{2}}(z) \qquad (2.13)$$

where $J_{n+\frac{1}{2}}(z)$ is the classical half-integer Bessel function. The prime ' indicates the derivative of the function relative to its augment.

The scattered wave in far field is spherical and its electric fields for the two polarized directions are respectively [101]:

$$E_{s,1} = \frac{i \exp(ikr)}{kr} S_1(\theta) E_{i,1} \qquad (2.14)$$

$$E_{s,2} = \frac{i \exp(ikr)}{kr} S_2(\theta) E_{i,2} \qquad (2.15)$$

$k = 2\pi/\lambda$ is the wave number. $E_{s,1}$ is the scattered far-field component in the scattering plane (OZP), and $E_{s,2}$ is the component perpendicular to this plane. $E_{i,1}$ and $E_{i,2}$ are incident electric field components of two polarized directions and we have

$$E_{i,1} = E_i \cos \varphi \qquad (2.16)$$
$$E_{i,2} = E_i \sin \varphi \qquad (2.17)$$

$S_1$ and $S_2$ are amplitude functions in the perpendicular and parallel polarizations and can be calculated by

$$S_1(\theta) = \sum_{n=1}^{\infty} \frac{2n+1}{n(n+1)} [a_n \pi_n(\cos \theta) + b_n \tau_n(\cos \theta)] \qquad (2.18)$$

$$S_2(\theta) = \sum_{n=1}^{\infty} \frac{2n+1}{n(n+1)} [a_n \tau_n(\cos \theta) + b_n \pi_n(\cos \theta)] \qquad (2.19)$$

$\pi_n$ and $\tau_n$ are Legendre functions defined by

$$\pi_n(\cos \theta) = \frac{P_n^1(\cos \theta)}{\sin \theta} \qquad (2.20)$$

$$\tau_n(\cos \theta) = \frac{dP_n^1(\cos \theta)}{d\theta} \qquad (2.21)$$

where $P_n^1(\cos \theta)$ are n-th order Legendre functions of first kind.

The corresponding far-field scattering intensities in the two polarized directions are [75]

$$I_1(\theta) = \frac{i_1}{k^2 r^2} \cos^2 \varphi I_0 \qquad (2.22)$$

$$I_2(\theta) = \frac{i_2}{k^2 r^2} \sin^2 \varphi I_0 \qquad (2.23)$$

where

$$i_1 = |S_1(\theta)|^2 \qquad (2.24)$$

$$i_2 = |S_2(\theta)|^2 \qquad (2.25)$$

$I_0$ is the incident intensity. In practice $i_1$ and $i_2$ are distance free intensities of scattered light and its plot is called the scattering diagram or scattering pattern.

Lorenz-Mie theory provides rigorous solutions to the Maxwell equations and are widely referred to validate the other methods. However, it cannot reveal the contribution of different scattering orders to the overall scattering pattern. To meet this demand, the Debye series is developed to study the contribution of different order of wave in rigorous regime. It rewrites each term of Mie series as another infinite series and clarifies the physical origins of many effects that occur in electromagnetic scattering [102, 103]. The classical Mie coefficients $a_n$ (TM amplitude) and $b_n$ (TE amplitude) are expanded in the Debye series to represent the separated contributions of the scattered waves:

$$\left. \begin{array}{l} a_n \\ b_n \end{array} \right\} = \frac{1}{2} \left[ 1 - R_n^{212} - \frac{T_n^{21} T_n^{12}}{1 - R_n^{121}} \right] \\ = \frac{1}{2} \left[ 1 - R_n^{212} - \sum_{p=1}^{\infty} T_n^{21} (R_n^{121})^{p-1} T_n^{12} \right] \qquad (2.26)$$
<figure><img src="image_8.png" /><figcaption>Figure 2.2: Debye series model of scattering by a homogeneous sphere.</figcaption></figure>

where as denoted in Fig. 2.2, $p$ is the order of emergent waves. $T_n^{21}$ and $T_n^{12}$ are partial-wave transmission coefficients, and $R_n^{21}$ and $R_n^{12}$ are partial-wave reflection coefficients. The superscripts 1 and 2 represent regions inside and outside the sphere respectively. These partial-wave transmission and reflection coefficients are defined by

$$T_n^{21} = -m \frac{2i}{D_n} \qquad (2.27)$$

$$R_n^{21} = \frac{\alpha \xi_n^{(2)'}(m_2ka) \xi_n^{(2)}(m_1ka) - \beta \xi_n^{(2)}(m_2ka) \xi_n^{(2)'}(m_1ka)}{D_n} \qquad (2.28)$$

$$T_n^{12} = -\frac{2i}{D_n} \qquad (2.29)$$

$$R_n^{12} = \frac{\alpha \xi_n^{(1)'}(m_2ka) \xi_n^{(1)}(m_1ka) - \beta \xi_n^{(1)}(m_2ka) \xi_n^{(1)'}(m_1ka)}{D_n} \qquad (2.30)$$

where

$$D_n = -\alpha \xi_n^{(1)'}(m_2ka) \xi_n^{(2)}(m_1ka) + \beta \xi_n^{(1)}(m_2ka) \xi_n^{(2)'}(m_1ka) \qquad (2.31)$$

$$\alpha = \begin{cases} 1 & \text{TE wave} \\ m & \text{TM wave} \end{cases} \qquad (2.32)$$

$$\beta = \begin{cases} m & \text{TE wave} \\ 1 & \text{TM wave} \end{cases} \qquad (2.33)$$

The superscripts (1) and (2) in the functions $\xi_n^{(1)}$ and $\xi_n^{(2)}$ denotes the kind of Ricatti-Hankel function.

Each term in the right-hand side of Eq. (2.26) has physical interpretation. The first term 1/2 corresponds to the contribution of the diffraction. The second term
$\frac{1}{2}R_n^{212}$ corresponds to the specular reflected wave. And the third term is the sum of refracted wave undergoing $p-1$ internal reflections. Therefore we have

$$ \left. \begin{array}{l} a_n(p=0) \\ b_n(p=0) \end{array} \right\} = \frac{1}{2} (1 - R_n^{212}) \qquad (2.34) $$

for the contribution from diffraction and the specular reflected wave, and

$$ \left. \begin{array}{l} a_n(p) \\ b_n(p) \end{array} \right\} = -\frac{1}{2} T_n^{21} (R_n^{121})^{p-1} T_n^{12} \qquad (2.35) $$

for the emergent wave of order $p$. When all the Debye series terms are added together, the results of Lorenz-Mie theory are recovered. This Debye series is integrated in ABSphere developed by K. F. Ren [104]. In the following of the thesis, all calculations concerning Debye series will be realized with ABSphere.

## 2.2 Geometrical optics

All rigorous solutions for scattering problems of a particle (spheroid, elliptical cylinder,...) are limited to the objects of simple shape that can be described in a mathematical coordinate system of the same geometry. But even in these simple cases, due to the difficulties in evaluating the special functions involved, the calculable size can hardly exceed a few tens of wavelengths except for a sphere or an infinite circular cylinder [67]. For particles whose dimension is much larger than the wavelength, geometrical optics (GO) is a good choice. It permits to deal with the scattering of particles with irregular shapes. We will recall firstly the basic concepts of geometrical optics. Then as a typical example, we will explain how to deal with the scattering problem with GO, which serves us as a basic knowledge for the further development of VCRM.

In GO, all waves are described by bundles of rays. When a ray enters the medium 1 from the medium 2 and arrives on the interface of the two media at an incident angle $\alpha$. It is partly reflected and refracted. The relationship between the incident angle $\alpha$ and refractive angle $\beta$ is given by Snell's law:

$$ \sin \alpha = m \sin \beta \qquad (2.36) $$

$m = m_1/m_2$ is the relative index of the medium 1 to the medium 2.

The direction of the reflected ray is determined by the law of reflection:

$$ \alpha = \alpha' \qquad (2.37) $$

$\alpha'$ is the angle of reflection. Both the reflected and refracted rays are in the same plane as the incident ray.
<figure><img src="image_9.png" /><figcaption>Figure 2.3: Light paths of ray interaction with a sphere.</figcaption></figure>

The ratio of the amplitude of the reflected wave to the incident wave is given by Fresnel formulas [56]:

$$r_1 = \frac{\cos \alpha - m \cos \beta}{\cos \alpha + m \cos \beta} \qquad (2.38)$$

$$r_2 = \frac{m \cos \alpha - \cos \beta}{m \cos \alpha + \cos \beta} \qquad (2.39)$$

where the subscript 1 and 2 of Fresnel coefficients stand for perpendicular and parallel polarization respectively.

When $\alpha + \beta = 90^\circ$, the Fresnel coefficient for parallel polarization $r_2$ is equal to zero, which means that the reflected wave at this angle is perpendicularly polarized. The corresponding incident angle $\alpha_B$ is called Brewster's angle and given by

$$\alpha_B = \tan^{-1} m \qquad (2.40)$$

It is $36.87^\circ$ for the relative refractive index $m = 0.75$ and $53.12^\circ$ for $m = 1.333$.

When the incident light propagates from an optically denser medium into the less dense medium, i.e. the relative index $m < 1$, and $\sin \alpha = m$, from Eq. (2.36) we know that $\sin \beta = 1$, i.e. $\beta = 90^\circ$. The ray emerges in a direction tangent to the boundary and when the incident angle is larger than this angle, the energy is totally reflected. The corresponding incident angle is called critical angle $\alpha_c$ and can be calculated by

$$\alpha_c = \sin^{-1} m \qquad (2.41)$$

In the following, we will apply the fundamental laws of GO (Snell's law and Fresnel formulas) to study the plane wave scattering by a homogeneous sphere. Fig. 2.3 presents light scattering by a homogeneous sphere of refractive index $m$ in the geometrical optics sense. A ray impinging on the sphere is partly reflected and partly
refracted. For each incident ray, there is an infinite series of emergent rays. And each emergent ray is characterized by the incident angle $\alpha$ and an integer $p$, called the order of ray. $p = 0$ stands for the ray reflected at the out surface of the sphere and $p = 1$ for the first order refracted ray. When $p \ge 2$, $p$ designates the order of rays undergoes $p - 1$ internal reflections.

### 2.2.1 Direction of emergent rays

The total deviation angle $\theta'$ from the original direction is given by Descartes [105]:

$$ \theta' = 2\tau - 2p\tau' \qquad (2.42) $$

or

$$ \theta' = 2p\beta - 2\alpha + (1-p)\pi \qquad (2.43) $$

where $\tau = \pi/2 - \alpha$, $\tau' = \pi/2 - \beta$. Because of the symmetry of the problem, in practice, we talk often about the scattering angle $\theta$ between 0 and $\pi$, such that:

$$ \theta = 2\pi l + q\theta' \qquad (2.44) $$

where $l$ is an integer and $q = +1$ or $-1$. They are chosen to make sure that $\theta$ is reported to the interval $[0, \pi]$.

### 2.2.2 Amplitude of scattered ray

The amplitude of scattered ray is affected by two factors in term of GO. One is the reflection and refraction on the particle surface and the other is the convergence factor of the wave on the curved surface of the particle. For a homogeneous sphere, both of them can be calculated analytically.

The amplitude ratio $\varepsilon_X$ of an emergent wave to the incident wave is calculated by the Fresnel formulas. For the external reflected ray ($p = 0$), this ratio is directly the Fresnel coefficient of reflection. For the internal reflection, the same fractions are reflected except the sign of $r_X$ is reversed. And the refracted part of the intensity is $1 - r_X^2$. Therefore the total intensity of a given incident pencil, is divided into the parts $r_X^2$ for the ray $p = 0$, $(1 - r_X^2)^2$ for the ray $p = 1$, $r_X^2(1 - r_X^2)^2$ for the ray $p = 2$, etc. All rays of order $p \ge 2$ undergo two transmissions and $p - 1$ internal reflections, so its amplitude ratio $\varepsilon_X$ can be calculated by.

$$ \varepsilon_X = \begin{cases} r_X & p = 0 \\ (1 - r_X^2)(-r_X)^{p-1} & p \ge 1 \end{cases} \qquad (2.45) $$

The subscript $X = 1$ or $2$ stands for the perpendicular or parallel polarization.
Consider a pencil of rays of the plane wave arriving at the sphere with an incident angle $\alpha$. The element of the illuminating surface is $dS_i = a^2 \cos \tau d\tau d\phi$. $\phi$ is the azimuth angle of incident ray. Let $I_0$ denote its intensity, then the incident energy flux is:

$$I_0 a^2 \cos \tau \sin \tau d\tau d\phi \qquad (2.46)$$

The emergent rays spread into a solid angle $d\Omega = \sin \theta d\theta d\phi$ which corresponds to a surface element $dS = r^2 d\Omega = r^2 \sin \theta d\theta d\phi$ in a distance $r$ from the sphere. According to the energy balance, the scattering intensity of order $p$ in GO is given by [56]

$$I_X(p, \alpha) = \frac{\varepsilon_X^2 I_0 a^2 \cos \tau \sin \tau d\tau d\phi}{r^2 \sin \theta d\theta d\phi} = \frac{a^2}{r^2} I_0 \varepsilon_X^2 \mathcal{D}_p \qquad (2.47)$$

$\mathcal{D}_p$ is usually called divergence factor, a parameter denoting the influence of the shape of a particle on the angular intensity distribution. $\mathcal{D}_p$ is given by

$$\mathcal{D}_p = \frac{\sin \alpha \cos \alpha}{\sin \theta |d\theta'/d\alpha|} \qquad (2.48)$$

According to Snell law $\sin \alpha = m \sin \beta$, we have

$$\frac{d\beta}{d\alpha} = \frac{\cos \alpha}{m \cos \beta} \qquad (2.49)$$

and the derivative of Eq. (2.43) yields

$$\frac{d\theta'}{d\alpha} = 2 \left( p \frac{d\beta}{d\alpha} - 1 \right) = 2 \left( p \frac{\cos \alpha}{m \cos \beta} - 1 \right) \qquad (2.50)$$

Therefore $\mathcal{D}_p$ can be calculated by:

$$\mathcal{D}_p = \frac{\sin \alpha \cos \alpha}{2 \sin \theta \left| p \frac{\cos \alpha}{m \cos \beta} - 1 \right|} = \frac{\sin \alpha \cos \alpha}{2 \sin \theta} \left| \frac{\sqrt{m^2 - \sin^2 \alpha}}{p \cos \alpha - \sqrt{m^2 - \sin^2 \alpha}} \right| \qquad (2.51)$$

For convenience, the amplitude and the phase of the scattered light can be expressed by means of the “amplitude functions” $S_1(\theta)$ and $S_2(\theta)$, which are firstly introduced in section 2.1 and used throughout this thesis. Therefore we can define the corresponding amplitude function of scattering light by:

$$S_{j,p}(\theta) = x \varepsilon_X \sqrt{\mathcal{D}_p} \exp(i\phi_p) \qquad (2.52)$$

The size parameter $x = 2\pi a/\lambda$. $\phi_p$ is the phase of emergent ray and will be introduced in next subsection.

The dimensionless intensity $i_1$ and $i_2$ of emergent ray of order $p$ in two polarizations are then:

$$i_{1,p} = |S_{1,p}(\theta)|^2 = x^2 \varepsilon_1^2 \mathcal{D}_p \qquad (2.53)$$

$$i_{2,p} = |S_{1,p}(\theta)|^2 = x^2 \varepsilon_2^2 \mathcal{D}_p \qquad (2.54)$$
2.2.3 Phase of emergent rays

The phase of emergent ray $\phi_p$ plays an important rule in the interference of different modes and the rays of the same mode undertaken different paths. According to van de Hulst [56], it is composed of three contributions:

1. the phase due to the reflection

2. the phase due to the optical path $\phi_{op}$

3. the phase due to the focal line $\phi_{fl}$

The reflection may change the sign of the amplitude of light and introduces a phase shift $\pi$. This possible change of the sign is already included in the definition of amplitude ratio $\varepsilon_X$ calculated by Fresnel formula. When total reflection occurs or for an absorbing medium, the Fresnel coefficients are complex and the phase shift can be calculated from the complex Fresnel formula.

The optical path is the product of the geometric length of the path of a ray, and the index of refraction of the medium through which it propagates. In the calculation of the phase due to the optical path, the emergent ray is compared to a reference ray, which arrives at the center of the particle in the same direction as the incident ray and goes out in the direction of emergent ray as if no particle exists. Referring to Fig. 2.3, we can get the difference in optical path of an emergent ray of order $p$ to the reference ray $\delta = 2a(\cos \alpha - pm \cos \beta)$. Therefore the phase due to the optical path $\phi_{op}$ can be calculated by

$$ \phi_{op} = k\delta = 2x(\cos \alpha - pm \cos \beta) \qquad (2.55) $$

According to van de Hulst, at the passage of any focal line, the phase advances by $\pi/2$. Thus to calculate the phase due to the focal line $\phi_{fl}$, the number of the focal lines encountered along the entire path must be counted. van de Hulst classified the focal lines in two types: type a for the crossing of two adjacent rays in the scattering plane and type b for the crossing of a ray with the central optical axis. It is found that the rays pass $p - (1 - s)/2$ focal lines of type a and $-2l + (1 - q)/2$ focal lines of type b, where s is equal to +1 or -1 depending to the sign of $d\theta'/d\tau$. l and q are defined in Eq. (2.44). The phase due to the focal line $\phi_{fl}$ is then given by:

$$ \phi_{fl} = \frac{\pi}{2} \left( p - 2l + \frac{s - q}{2} \right) \qquad (2.56) $$

The total phase of a ray is then

$$ \phi_p = \phi_{op} + \phi_{fl} \qquad (2.57) $$
### 2.2.4 Total complex amplitude

The total amplitude of scattered light at a given angle $\theta$ is calculated by the summation of the complex amplitude of all order rays as well as the diffraction:

$$S_j = S_d(\theta) + \sum_{p=0}^{\infty} S_{j,p}(\theta) \qquad (2.58)$$

$S_d(\theta)$ is the amplitude of diffraction which can be calculated by using the Fraunhofer diffraction of the projection circular disk of the object on a plane perpendicularly to the incident wave.

## 2.3 Rainbows and Airy theory

In the framework of geometrical optics, a rainbow corresponds to an extreme of deviation for rays experiencing internal reflections. This means the divergence factor $D_p$ defined in Eq. (2.51) trends to infinity, i.e.

$$\frac{d\theta'}{d\alpha} = 2 \left( p \frac{\cos \alpha}{m \cos \beta} - 1 \right) = 0 \qquad (2.59)$$

Combining with Snell's law, we can obtain incident angle $\alpha_{rg}$ of rainbow angle by

$$\sin \alpha_{rg} = \sqrt{\frac{p^2 - m^2}{p^2 - 1}} \qquad (2.60)$$

and the refraction angle $\beta_{rg}$ by

$$\sin \beta_{rg} = \sqrt{\frac{p^2 - m^2}{m^2(p^2 - 1)}} \qquad (2.61)$$

Their complement can be then calculated from:

$$\tau = \tan^{-1} \sqrt{\frac{m^2 - 1}{p^2 - m^2}} \qquad (2.62)$$

$$\tau' = \tan^{-1} \sqrt{\frac{p^2(m^2 - 1)}{p^2 - m^2}} \qquad (2.63)$$

By substituting Eqs. (2.62) and (2.63) into Eqs. (2.42) and (2.44), we can get the rainbow angle $\theta_{rg}$ by [75]

$$\theta_{rg}(p, m) = 2\pi l + q \left[ 2 \arctan \left( \sqrt{\frac{m^2 - 1}{p^2 - m^2}} \right) - 2p \arctan \left( \sqrt{\frac{p^2(m^2 - 1)}{p^2 - m^2}} \right) \right] \qquad (2.64)$$
<figure><img src="image_10.png" /><figcaption>Figure 2.4: The geometric rainbow angle for first, second and third order rainbows as function of refractive index.</figcaption></figure>

Then for the first order rainbow ($p = 2$), we have

$$ \theta_{rg} = 6 \tan^{-1} \sqrt{\frac{m^2 - 1}{4 - m^2}} \qquad (2.65) $$

and for the second order rainbow ($p = 3$)

$$ \theta_{rg} = 2\pi - 16 \tan^{-1} \sqrt{\frac{m^2 - 1}{9 - m^2}} \qquad (2.66) $$

For a sphere of relative refractive index $m = 1.333$, the primary and the secondary rainbow angles are 137.922° and 129.109° respectively.

From the above equations, we can find that rainbow angles in GO are function of refractive index and independent of particle size. Fig. (2.4) depicts the geometric rainbow angles of the first, second, and third order rainbows ($p = 2, 3$ and 4 respectively) as function of the refractive index $m$. From the figure we can find that when the refractive index $m$ is rather small, the angular position of the first order rainbow is smaller than that of the second order rainbow. With increase of the refractive index, the geometric rainbow angle of the first order rainbow increases while that of the second order rainbow decreases. When refractive index $m = 1.312$, they have the same value.

Geometrical optics (GO) permits to explain clearly the formation of rainbows of different orders and give a clear physical picture of mechanism of light scattering. However GO cannot predict the intensity at the rainbow angles. Furthermore, the position of the maximum intensity in a rainbow is offset slightly from the geometrical rainbow angle $\theta_{rg}$. Most of all it cannot explain the existence of series of maxima and minima associated with rainbows (known as supernumerary structure). To overcome these problems, Airy introduced the wave nature of light. The idea of Airy in dealing with this problem can be explained by Figure 2.5. A homogeneous sphere is illuminated by a plane wave. Five rays around the first rainbow angles are shown here and
<figure><img src="image_11.png" /><figcaption>Figure 2.5: Coordinate system for Airy’s theory of the rainbow.</figcaption></figure>

the path of rays is calculated by software VCRMell2D [106] developed in laboratory. The middle ray (in red) is the one with the minimum deflection and leaves the sphere at the geometric rainbow angle. All the other incident rays in its neighborhood (number 1, 2 in blue and 4, 5 in green) emerge from the drop in the same side of the extreme ray. Because these rays from two side of the extreme ray emergent in the same side and experience different optical path. They will interfere and the interference may be constructive or destructive. This leads to the supernumerary structure.

In order to understand what’s going on in simple terms and calculate the intensity of scattered light, Airy used a virtual cubic wavefront to replace the true wavefront, and assumes that the outgoing wave comes from this virtual plane perpendicular to the rainbow ray. The total field in a given direction is the integration of the field on this plane. To obtain the phase function, Airy used the Taylor series expansion of scattering angle in the neighborhood of the rainbow angle. He designed a Cartesian coordinate system *Ouv*. The *u*-axis is along the extreme rainbow ray and the *v*-axis is perpendicular to the rainbow ray and passing through the center of the sphere. Suppose that $\theta_{rg}$ is the rainbow angle in GO and $\alpha_{rg}$ is its incident angle which can be calculated from Eq. (2.60). $\tau_{rg} = \pi/2 - \alpha_{rg}$ is the complement of the incident angle. When the incident angle $\alpha = \alpha_{rg}$, we have $\theta = \theta_{rg}$. Because a rainbow corresponds to an extreme of the deviation angle, we have:

$$ \frac{d\theta}{d\tau} = 0 \qquad (2.67) $$

Then by further differentiation we have

$$ \frac{d^2\theta}{d\tau^2} = \frac{3}{2 \tan \tau_{rg}} \qquad (2.68) $$
Therefore the scattering angle in the neighborhood of the rainbow angle can be given in the sense of Taylor expansion by

$$ \theta - \theta_{rg} = \frac{3}{4 \tan \tau_{rg}} (\tau - \tau_{rg})^2 \qquad (2.69) $$

Besides, the slope of the wavefront in the neighborhood of $\theta_{rg}$ is

$$ \frac{du}{dv} = \theta - \theta_{rg} \qquad (2.70) $$

and

$$ v = a \cos \tau - a \cos \tau_{rg} \approx -a(\tau - \tau_{rg}) \sin \tau_{rg} \qquad (2.71) $$

Finally we have

$$ \frac{du}{dv} = \frac{3v^2}{4a^2 \sin^2 \tau_{rg} \tan \tau_{rg}} \qquad (2.72) $$

The equation of the cubic wavefront can be got by integration:

$$ u = \frac{hv^3}{3a^2} \qquad (2.73) $$

where h is a constant

$$ h = \frac{(p^2 - 1)^2 (p^2 - m^2)^{1/2}}{p^2 (m^2 - 1)^{3/2}} \qquad (2.74) $$

Airy further assumes that the amplitude along the wavefront is constant. The far field amplitude of the wave in the direction $\theta$ is the integral of the contributions at all the points on the cubic wavefront. He extend the integration on $v$ from $-\infty$ to $+\infty$. So the amplitude of the scattered wave is proportional to

$$ \int_{-\infty}^{\infty} e^{-ikv(\theta-\theta_{rg})+ikhv^3/3a^2} dv \qquad (2.75) $$

By combining the integrands for positive and negative $v$, Airy introduced the “rainbow integral”:

$$ f(z) = \int_{0}^{\infty} \cos \frac{1}{2} \pi (zt - t^3) dt \qquad (2.76) $$

where

$$ t = v/l \qquad (2.77) $$

$$ l = (3\lambda a^2 / 4h)^{1/3} \qquad (2.78) $$

$$ z = 4l(\theta - \theta_{rg})/\lambda = (12/h\pi^2)^{1/3} x^{2/3}(\theta - \theta_{rg}) \qquad (2.79) $$

The intensity is proportional to $f^2(z)$. The Airy intensity for the rainbow of $p-1$ order can be calculated by [102]

$$ I_{Airy}^{p}(\theta) = |E_{Airy}^{p}(\theta)|^2 \qquad (2.80) $$
<figure><img src="image_12.png" /><figcaption>Figure 2.6: Intensity distribution calculated by Airy theory for a sphere of refractive index $m = 1.333$ and radius $a = 50$ µm</figcaption></figure>

where $E_{Airy}^p(\theta)$ is the electric field in the vicinity of $\theta_{rg}$ and is given by

$$E_{Airy}^{p}(\theta) = x \left( \frac{2\pi \sin \alpha_{rg}}{\sin \theta_{rg}} \right)^{1/2} \frac{x^{1/6}}{h^{1/3}} t^{21}(\alpha_{rg}) \times [r^{11}(\alpha_{rg})]^{p-1} t^{12}(\alpha_{rg}) \\ \times Ai\left(\frac{-x^{2/3}\Delta}{h^{1/3}}\right) \times \exp(2\pi i L^R / \lambda) \exp\left[(ix\Delta)\left(\frac{p^2 - m^2}{p^2 - 1}\right)^{1/2}\right] \quad (2.81)$$

where $\Delta$ in Eq. (2.81) is the angular difference between observation (scattering) angle $\theta$ and rainbow angle $\theta_{rg}$

$$\Delta = \theta - \theta_{rg} \qquad (2.82)$$

$L^R$ is the optical path of the rainbow ray:

$$L^R = 2a(pm \cos \beta_{rg} - \cos \alpha_{rg}) + 2a \qquad (2.83)$$

$\beta_{rg}$ is the refractive angle corresponding to the incident angle $\alpha_{rg}$, and is given in Eq. (2.61). $Ai$ is the Airy function defined by

$$Ai(x) = \frac{1}{\pi} \int_{0}^{\infty} \cos \left( \frac{t^3}{3} + xt \right) dt \qquad (2.84)$$

$t^{21}$ is the transmission factor of the electric field from the exterior to the interior of the sphere. $t^{12}$ is the transmission factor from the interior to the exterior and $r^{11}$ is the internal reflection coefficient.

As an example, Fig. 2.6 shows the intensity near the primary rainbow angle as function of scattering angle calculated by Airy theory. Supernumerary structure can be observed clearly. Geometrical rainbow angle is also indicated for comparison.
The positions of the Airy maxima $\theta_{Airy,K}$ can also be found according to the properties of Airy functions [75]

$$|\theta_{Airy,0} - \theta_{rg}(p, m)| = 1.087376 \left(\frac{\pi^2}{12x^2} h\right)^{1/3} \qquad (2.85)$$

for main peak $K = 0$, and

$$|\theta_{Airy,K} - \theta_{rg}(p, m)| = \left[ \frac{9\pi^2}{4x^2} \left( K + \frac{1}{4} \right)^2 h \right]^{1/3} \qquad (2.86)$$

for the subsidiary supernumerary bows ($K = 1, 2, 3...$). $h$ depends on the refractive index $m$ and $x$ is the size parameter. Therefore it is clear that the angular positions of supernumerary bows depend on the particle size and its refractive index.

Besides from Eqs. (2.85) and (2.86), we deduce that

$$\theta_{rg}(p, m) = \frac{\theta_{Airy,K} - c\theta_{Airy,0}}{1 - c} \qquad (2.87)$$

where the constant $c$ is defined by

$$c = \frac{1}{1.087376} \left[ 27 \left( K + \frac{1}{4} \right)^2 \right]^{1/3} \qquad (2.88)$$

Similar formula can be derived for the determination of geometrical rainbow angle by two other Airy bows $K_1$ and $K_2$.

$$\theta_{rg}(p, m) = \frac{\theta_{Airy,K_1} - d\theta_{Airy,K_2}}{1 - d} \qquad (2.89)$$

where constant $d$ is defined by

$$d = \left( \frac{K_1 + \frac{1}{4}}{K_2 + \frac{1}{4}} \right)^{2/3} \qquad (2.90)$$

Both $c$ and $d$ are independent of the refractive index, the particle size and the order of rainbow. This means that they can be applied to determine the geometrical rainbow angles from two Airy bows of any order of rainbow for a particle of any size and material.

## 2.4 Conclusion and comments

Two typical classical methods for dealing with the light scattering by a small particle are presented.
The first one, Lorenz-Mie theory, is a rigorous solution of the Maxwell equation, which will serve us as a reference to validate the model developed in this thesis and to evaluate its precision. Especially its series version - Debye theory is particular important for the development of VCRM since it permits to decompose different modes of scattering, very suitable to compare directly with the results of VCRM.

Geometric optics approximation is an approximate method, but it is very flexible allowing to deal the scattering of non-spherical particles. It is largely used in the optical particle characterization, directly or indirectly. It is also the base of our VCRM. The Airy theory, also an approximate method, has been applied naturally since long time in the rainbow refractometry. Its precision is so satisfactory that few researchers intend to examine the approximation and calculation involved in the theory. We present this theory in details for two purposes:

* To provide a basic knowledge to understand and apply VCRM in the characterization of non-spherical particles.

* To reexamine it later in the framework of VCRM to show that the Airy theory can be further improved, and the idea of Airy serves us to deal with non-spherical particle.

We will present in the next chapter the principle of VCRM and its application to the characterization of non-spherical particles.
# Chapter 3

## VCRM for scattering by an ellipsoid

Geometrical optics is a simple and intuitive method for treating the interaction of an object with electromagnetic or light waves when the dimension of the object is much larger than the wavelength. It can be applied to objects of complex shape, which are hard or even impossible to be dealt with by rigorous theories or most numerical techniques. Many researchers have contributed to its improvement. Some take into account the forward diffraction or other particular wave effects (Airy theory for the rainbow [56, 102, 107] and Marstons model for the critical scattering [108, 109]. Others combine directly geometrical optics with the electromagnetic wave method [65, 66]. However, in these studies, interference effects of all order rays are rarely taken into account. Besides, when ray optics is extended to a three dimensional (3D) object of irregular shape, it becomes a heavy task (see [110–113] and references therein) because of the difficulties in the tracing of rays in a complex shaped particle, the calculation of local divergence factors and the phase shift due to focal lines.

To overcome these difficulties, a new model of light-particle interaction, called Vectorial Complex Ray Model (VCRM) is developed [67]. It is based on the geometrical optics but a new property – the wave front curvatures is introduced. In this model, all waves are described by vectorial complex rays. The ray direction and the wave divergence/convergence after each interaction of the wave with a dioptric surface as well as the phase shifts of each ray are determined by the vector Snell law and the wavefront equation according to the curvatures of the surfaces. The total scattered field is the superposition of the complex amplitude of all orders of the rays emergent from the object.

In this chapter, we will firstly describe Vectorial Complex Ray Model in general case, i.e. reflection and refraction of an arbitrarily shaped wave by a dioptric surface in arbitrary direction. Then VCRM is applied to predict the light scattering by an ellipsoid.
## 3.1 Vectorial Complex Ray Model

In VCRM, the wave is considered as bundles of vectorial complex rays. And each ray possesses five properties:

* direction of propagation $\hat{k}$

* polarization state $X$

* amplitude $A_{p,X}$

* phase $\phi_p$

* wavefront curvature

where the wavefront curvature is newly introduced compared with classical ray models making VCRM easy to extend to irregularly shaped 3D objects. The complex amplitude $S_{p,X}$ of an emergent ray can be then expressed as

$$S_{p,X} = A_{p,X} \exp(i\phi_p) \qquad (3.1)$$

where subscripts $X = 1$ or $2$ correspond respectively to the perpendicular and parallel polarization to the plane of incidence. $p$ is the order of the emergent ray.

The direction of the rays before and after reflection or refraction are described respectively by the normalize wave vectors $\hat{k} = k/k'$ and $\hat{k}' = k'/k'$. The two wave vectors $k$ and $k'$ satisfy the vector Snell's law:

$$(k - k') \times n = 0 \qquad (3.2)$$

with $n$ the normal vector of the dioptric surface at incident point.

Two factors play rule in the amplitude of the emergent field: divergence/convergence of the waves and reflection/refraction on the dioptry surface. In VCRM, the divergence/convergence of the wave, described by the divergence factor, is calculated directly by the wave front curvature radii. The divergence factor $D_p$ of an emergent ray after $q = p + 1$ interactions with the dioptry surface can be calculated by

$$D_p = \frac{R'_{11}R'_{21}}{R_{12}R_{22}} \cdot \frac{R'_{12}R'_{22}}{R_{13}R_{23}} \cdots \frac{R'_{1,q}R'_{2,q}}{(r + R'_{1,q})(r + R'_{2,q})} \qquad (3.3)$$

where $r$ is the distance from the emergent point to the observation point, $R_{1,q}$ and $R_{2,q}$ are the two principal wavefront curvature radii of rays before interaction. $R'_{1,q}$ and $R'_{2,q}$ are those after interaction (reflection or refraction). The curvature radii are determined by the curvature matrix of the corresponding surface. Suppose that an arbitrary wave of curvature matrix $Q$ impinges on a dioptric surface of curvature
matrix C. Then the curvature matrix $Q'$ of the wave after intersection can be deduced from the wavefront matrix equation [67]

$$ (k' - k) \cdot nC = k'\Theta'^T Q'\Theta' - k\Theta^T Q\Theta \quad (3.4) $$

where the letters with prime represent the quantities after refraction or reflection, the superscript T represents the transpose of the matrix, and Θ is the projection matrix between the base unitary vectors of the coordinates systems on the planes tangent to the wavefront ($s_1, s_2$) and to the dioptric surface ($t_1, t_2$)

$$ \begin{pmatrix} s_1 \cdot t_1 & s_1 \cdot t_2 \\ s_2 \cdot t_1 & s_2 \cdot t_2 \end{pmatrix} $$

$\epsilon_{p,X}$ stands for the relative amplitude of the X polarized ray of order p and is determined by Fresnel coefficients:

$$ \epsilon_{p,X} = \begin{cases} r_{0,X} & p = 0 \\ t_{0,X} t'_{p,X} \prod_{n=1}^{p-1} r'_{n,X} & p \ge 1 \end{cases} \quad (3.5) $$

$r_{0,X}$ and $t_{0,X}$ are respectively the Fresnel reflection and refraction coefficients for a ray impinging the particle surface from surrounding medium, and $r'_{n,X}$ and $t'_{n,X}$ are those for a ray arriving at the dioptric interface from inside of the particle. To be consistent with VCRM and simply the calculation, the Fresnel coefficients for reflected ray in Eqs. (2.38) and (2.39) can be written as function of the normal components of the wave vector:

$$ r_1 = \frac{k_n - k'_n}{k_n + k'_n} \quad (3.6) $$

$$ r_2 = \frac{m^2 k_n - k'_n}{m^2 k_n + k'_n} \quad (3.7) $$

and the amplitude ratio of transmitted waves $t_X$ can be calculated by

$$ \begin{aligned} t_1 &= \sqrt{1 - r_1^2} \\ &= \frac{2\sqrt{k_n k'_n}}{k_n + k'_n} \end{aligned} \quad (3.8) $$

$$ \begin{aligned} t_2 &= \sqrt{1 - r_2^2} \\ &= \frac{2m\sqrt{k_n k'_n}}{m^2 k_n + k'_n} \end{aligned} \quad (3.9) $$

**Attention:** It should be noted that for a non-spherical particle, the Fresnel coefficients must be calculated at each intersection point because the incident and refractive angles change each time a ray encounters the surface.
The amplitude $A_{p,X}$ of the emergent ray is product of Fresnel factor $\epsilon_X$ and the square root of the divergence factor $D_p$

$$A_{p,X} = |\epsilon_{p,X}|\sqrt{D_p} \qquad (3.10)$$

The phase of each emergent ray $\phi_p$ is composed of four parts: the phase of the incident wave $\phi_{inc}$; the phase due to the optical path $\phi_{op}$; the phase due to the focal point or focal line $\phi_{fl}$ and the phase due to the reflection $\phi_{rl,X}$. The first one is calculated directly by the phase function of the incident electromagnetic wave expression. The phase due to the optical path $\phi_{op}$ can be computed according to the optical trajectory. Each time the sign of the curvature radius changes, a phase shift $\pi/2$ is added for the phase due to the focal point or focal line $\phi_{fl}$. And the phase due to the reflection is included in the Fresnel coefficient.

Knowing the amplitude and the phase of each ray, we calculate the total scattered field by the summation of the complex amplitudes of all emergent rays and the diffraction field.

## 3.2 Scattering in a symmetric plane of an ellipsoid

Consider an ellipsoidal particle of refractive index $m$ illuminated by a plane wave propagating in the symmetric plane of the ellipsoid. Suppose that the semi-axes of the ellipsoid in $x$, $y$ and $z$ directions are respectively $a$, $b$ and $c$. And the incident plane wave propagates in the $xz$ plane and making an angle $\theta_i$ ($0 \le \theta_i \le \frac{\pi}{2}$) with the $z$ axis as shown in Fig. 3.1. Here we only focus on the scattering in $xz$ plane (symmetric plane). The problem is therefore a 2D scattering since the rays remain always in this plane. But the curvature of the particle surface in the direction perpendicular to this plane counts in the VCRM. This is a key difference of VCRM from the classical geometrical optics.

### 3.2.1 Ray tracing

The cross section of the ellipsoid in this plane is described by

$$\frac{x^2}{a^2} + \frac{z^2}{c^2} = 1 \qquad (3.11)$$

#### Incident rays

In VCRM, the plane wave is considered as bundles of vectorial complex rays. The linear equation in the slope-intercept form for each incident ray in the $xz$ plane is
<figure><img src="image_13.png" /><figcaption>Figure 3.1: Schema of an ellipsoid illuminated by a plane wave</figcaption></figure>

then

$$x = \gamma z + d \qquad (3.12)$$

where $\gamma = \tan \theta_i$ is the slope of the ray and $d$ is the x-intercept of the incident ray. Substituting Eq. (3.12) into Eq. (3.11), we can get

$$\frac{(\gamma z + d)^2}{a^2} + \frac{z^2}{c^2} = 1 \qquad (3.13)$$

To determine the illumination region, we rewrite the Eq. (3.13) in the quadratic equation in one unknown $z$

$$(c^2\gamma^2 + a^2)z^2 + 2\gamma dc^2z + (c^2d^2 - a^2c^2) = 0 \qquad (3.14)$$

The discriminant $\Delta$ of this quadratic equation is

$$\Delta = 4\gamma^2 d^2 c^4 - 4(c^2\gamma^2 + a^2)(c^2d^2 - a^2c^2) \qquad (3.15)$$

When $\Delta = 0$, Eq. (3.13) has one solution, corresponding to the rays tangent to the ellipsoid surface. The two solutions of $d$ in $\Delta = 0$ are the two $x$-intercepts of the incident rays (see Fig. 3.1)

$$d_{max,min} = \pm \sqrt{c^2\gamma^2 + a^2} \qquad (3.16)$$
The rays with x-intercept $d_{min} < d < d_{max}$ interact the ellipsoid with the coordinate

$$z = \frac{-2\gamma dc^2 \pm \sqrt{\Delta}}{2(c^2\gamma^2 + a^2)} \qquad (3.17)$$

Substitute Eq. (3.17) in Eq. (3.12), we can get the coordinate $x$ of intersection point of the ray with the ellipsoid.

**Attention:** in our code, we suppose that the rays arrive on the ellipsoid from negative axis of $z$ to positive axis and $0 \le \theta_i \le \pi/2$, so only the solution with the smaller $z$-value meets our demand.

## Interaction point

The refracted ray after the incident point will propagate and interact with the particle surface, and is split to the reflected ray and refracted ray. The internal reflected ray will continue to propagate and interact with the particle surface. Suppose that the two successive interaction points are $(z_{j-1}, x_{j-1})$ and $(z_j, x_j)$, and the slope of the ray is $\gamma_{j-1}$, then we have:

$$x_j - x_{j-1} = \gamma_{j-1}(z_j - z_{j-1}) \qquad (3.18)$$

Substituting this equation into Eq. (3.11), we can get $z_j$ of next interaction point with the ellipsoid.

$$z_j = z_{j-1} - \frac{2\gamma_{j-1}x_{j-1}c^2 + z_{j-1}a^2}{c^2\gamma_{j-1}^2 + a^2} \qquad (3.19)$$

Knowing the coordinate $(z_{j-1}, x_{j-1})$ of the $(j-1)_{th}$ interaction point, the coordinate $(z_j, x_j)$ of the next interaction point are given by Eqs. (3.18) and (3.19).

## Wave vectors

In VCRM, the direction of a ray is described by the wave vector. When a ray arrives at the particle surface, its wave vector can be expressed by

$$\mathbf{k} = k_x \mathbf{e}_x + k_z \mathbf{e}_z \qquad (3.20)$$

where $k_x$ and $k_z$ are respectively the components of the wave vector $\mathbf{k}$ in $x$ and $z$ directions. The components of the wave vector in the normal and tangent directions of the particle surface can be expressed as

$$k_n = \mathbf{k} \cdot \mathbf{n} = k_x n_x + k_z n_z \qquad (3.21)$$

$$k_{\tau} = \mathbf{k} \cdot \mathbf{\tau} = -k_z n_x + k_x n_z \qquad (3.22)$$

where $\mathbf{n}$ and $\mathbf{\tau}$ are the unit vectors, respectively, normal and tangent to the particle surface. Their calculations will be given below. According to the Snell-Descartes law
in Eq. (3.2), we can get the tangent components of the reflected wave vector $k_{\tau}^l$ and the refracted wave vector $k_{\tau}'$ are equal to that of the incident wave vector

$$k_{\tau} = k_{\tau}^{l} = k_{\tau}' \qquad (3.23)$$

The normal component of the reflected wave vectors is

$$k_{n}^{l} = -k_{n} \qquad (3.24)$$

and the reflected wave vector is then

$$\mathbf{k}^{l} = k_{\tau}^{l} \boldsymbol{\tau} + k_{n}^{l} \boldsymbol{n} = k_{\tau} \boldsymbol{\tau} - k_{n} \boldsymbol{n} \qquad (3.25)$$

It can also be given in $x$ and $z$ components:

$$\mathbf{k}^{l} = k_{x}^{l} \mathbf{e}_{x} + k_{z}^{l} \mathbf{e}_{z} \qquad (3.26)$$

where

$$k_{x}^{l} = k_{n}^{l} n_{x} + k_{\tau}^{l} \tau_{x} \qquad (3.27)$$

$$k_{z}^{l} = k_{n}^{l} n_{z} + k_{\tau}^{l} \tau_{z} \qquad (3.28)$$

For the refracted wave vector $\mathbf{k}'$, we have

$$k' = mk \qquad (3.29)$$

$$\mathbf{k}' = k_{\tau}' \boldsymbol{\tau} + k_{n}' \boldsymbol{n} \qquad (3.30)$$

where its normal component

$$k_{n}' = \sqrt{k'^{2} - (k'_{\tau})^{2}} \qquad (3.31)$$

### Normal and tangent vectors of ellipsoid

Suppose that our solution for the coordinate of the interaction point in $xz$ plane is $(z, x)$, and the outward-pointing unit normal vector $\mathbf{n}$ of the surface of ellipsoid in $xz$ plane is the normalized gradient of the surface equation $f(x, z) = \frac{x^2}{a^2} + \frac{z^2}{c^2} - 1$, i.e.

$$\mathbf{n} = \frac{\frac{2x}{a^2} \mathbf{e}_x + \frac{2z}{c^2} \mathbf{e}_z}{\sqrt{(\frac{2x}{a^2})^2 + (\frac{2z}{c^2})^2}} = \frac{c^2 x \mathbf{e}_x + a^2 z \mathbf{e}_z}{\sqrt{c^4 x^2 + a^4 z^2}} = n_x \mathbf{e}_x + n_z \mathbf{e}_z \qquad (3.32)$$

It should be noted that for the first refraction and reflection of the incident rays on the ellipsoid, the unit normal vector of the surface should be reversed, i.e. the outward-pointing normal vector is changed to the inward-pointing normal vector. Then the unit tangent vector to the surface at this point is

$$\boldsymbol{\tau} = \tau_x \mathbf{e}_x + \tau_z \mathbf{e}_z = -n_z \mathbf{e}_x + n_x \mathbf{e}_z \qquad (3.33)$$

Even if $n_z \mathbf{e}_x - n_x \mathbf{e}_z$ is also the unit tangent vector to the surface, we chose the definition in Eq. (3.33) to make sure that our reflection angle and refraction angle are in the interval $[0, \pi/2]$.
### 3.2.2 Amplitudes and phases of emergent rays

In VCRM, the complex amplitude $S_{p,x}$ is calculated step by step at each interaction of the rays with the particle. Its module $A_{p,x}$ is determined by three factors: incident wave amplitude, divergence factor and the Fresnel coefficients. For the plane wave incidence, the amplitude and the phase of incident rays are constant. Here the calculations of divergence factor and the Fresnel coefficients will be presented. The phase shifts will be discussed later.

#### Divergence factor

When the rays encounter the curved dioptric surface, it will be converged or diverged. This convergence or divergence, called the divergence factor in VCRM, is essential to calculate the amplitude of the emergent rays and the phase shift due to the focal line. In VCRM, the divergence factor $D_p$ defined in Eq. (3.3) is calculated directly by the wavefront curvature radii. To get the curvature radii of the wavefront after each interaction, the wavefront matrix in Eq. (3.4) can be written as two scalar equations due to the symmetry of the problem [67]

$$ \frac{k_n'^2}{k'R_{1,q}'} = \frac{k_n^2}{kR_{1,q}} + \frac{k_n' - k_n}{\rho_1} \qquad (3.34) $$

$$ \frac{k'}{R'_{2,q}} = \frac{k}{R_{2,q}} + \frac{k'_{n} - k_{n}}{\rho_{2}} \qquad (3.35) $$

For a plane wave, its curvature radii of wavefront are infinity. The two principal curvature radii $\rho_1$ and $\rho_2$ of the ellipsoid surface at the plane $y = 0$ are given by

$$ \rho_1 = -\frac{c^2}{a} \left[ 1 + \left( \frac{a^2}{c^2} - 1 \right) \frac{z^2}{c^2} \right]^{3/2} \qquad (3.36) $$

$$ \rho_2 = -\frac{b^2}{a} \left[ 1 + \left( \frac{a^2}{c^2} - 1 \right) \frac{z^2}{c^2} \right]^{1/2} \qquad (3.37) $$

The curvature radii can be positive or negative depending on its convergent or divergent property. Here we adopt the convention that the curvature radii are positive if the focal line/point is after the intersection point.

Therefore the divergence factor $D_0$ of the reflected wave $p = 0$ is given by

$$ D_0 = \frac{R''_{11} R''_{21}}{(r + R''_{11})(r + R''_{21})} \qquad (3.38) $$

$R''_{11}$ and $R''_{21}$ are wave front curvature radii of the reflected wave. To make a difference with the wave front curvature radii of the refracted wave, in the following equations
we use $R''_{1,q}$ and $R''_{2,q}$ to specialize the wavefront curvature radii of the reflected rays. And the divergence factor $D_1$ of the emergent ray $p = 1$ is

$$D_1 = \frac{R'_{11} R'_{21}}{R_{12} R_{22}} \cdot \frac{R'_{12} R'_{22}}{(r + R'_{12})(r + R'_{22})} \qquad (3.39)$$

where $R_{12}$ and $R_{22}$ are the curvature radii of incident wave before the first internal reflection.

$$R_{12} = R'_{11} - d \qquad (3.40)$$

$$R_{22} = R'_{21} - d \qquad (3.41)$$

$d$ is the distance between the two successive interaction points. And $R'_{12}$ and $R'_{22}$ are the curvature radii of the first order refraction ray.

### Fresnel coefficients

The factor $\epsilon_{p,X}$ of each emergent ray can be calculated by Eq. (3.5).

### Phase of emergent rays

To take into account the effect of interference, generally, the phase of each ray must be computed, which is composed of four parts:

1. the phase of the incident wave $\phi_i$

2. the phase due to the optical path $\phi_{op}$

3. the phase due to the focal point or focal line $\phi_{fl}$

4. the phase due to the reflection $\phi_{rl,X}$

The phase $\phi_i$ of the incident wave is calculated directly according to the wavefront of the incident wave expression. It is a constant for the plane wave and plays no role in the interference, so we neglect it in the calculation.

The phase $\phi_{op}$ due to the optical path is calculated according to the optical path in and out of the particle relative to the reference ray, which propagates in the vacuum and arrives at the center of the particle in the same direction as the incident ray, then goes out in the direction of the emergent ray. Therefore, the reflected ray ($p = 0$) with a shorter path than the reference rays has a positive phase shift. The refracted rays ($p = 1, 2, 3, \dots$) have a longer path and negative phase shift.

The phase due to the focal lines indicates the phase shifts due to the crossing of two adjacent rays or crossing of a ray with the central optical axis. Each time a
ray passes a focal line, the phase advances $\pi/2$ according to van de Hulst [56]. This
holds also for a focus point where two focal lines coincide, so the phase advances by
$\pi$. It is easy to calculate the phase due to focal lines for a spherical particle or a
circular cylinder, in which, for a given incident ray, the deviation angles of all orders
can be obtained analytically and the focal lines or focal points can also be counted
accordingly. However for a particle of complex form, it's impossible to obtain an
analytical expression of phase due to the focal line. In VCRM we integrate the wave
property in the ray model and then the phase shifts due to the focal lines can be
determined directly by the change of sign of the wavefront curvature. Each time the
curvature (radius) changes the sign, we add a phase shift $\pi/2$.

The phase due to reflection is included in the Fresnel coefficient. An additional
phase $\pi$, is added when the Fresnel coefficient is negative. But when the incident light
is traveling from an optically denser medium towards an optically less dense medium
and the incident angle is greater than the so-called critical angle $\alpha_c = \sin^{-1} m$, the
total reflection occurs and the energy is totally reflected. In this case, the wave number
in the second medium is less than the tangent component of wave vector in the first
medium, i.e. $k' < k_{\tau}$. Fresnel coefficient will be a complex number. We can write
$k'_n = i\sqrt{k'^2 - k'^2}$. The Eqs. (3.6) and (3.7) are turned to be:

$$
\begin{align}
r_1 &= \frac{k_n - i\sqrt{k_{\tau}^2 - k'^2}}{k_n + i\sqrt{k_{\tau}^2 - k'^2}} \nonumber \\
   &= \frac{(k_n^2 - k_{\tau}^2 + k'^2) - 2k_n\sqrt{k_{\tau}^2 - k'^2}i}{k_n^2 + k_{\tau}^2 - k'^2} \tag{3.42}
\end{align}
$$

$$
\begin{align}
r_2 &= \frac{m^2 k_n - i \sqrt{k_{\tau}^2 - k'^2}}{m^2 k_n + i \sqrt{k_{\tau}^2 - k'^2}} \nonumber \\
    &= \frac{(k_{\tau}^2 - k'^2 - m^4 k_n^2) + 2m^2 k_n \sqrt{k_{\tau}^2 - k'^2} i}{m^4 k_n^2 + k_{\tau}^2 - k'^2} \tag{3.43}
\end{align}
$$

The phase shifts of the reflected rays are then

$$
\phi_{rl,1} = \arctan \frac{-2k_n \sqrt{k_{\tau}^2 - k'^2}}{k_n^2 - k_{\tau}^2 + k'^2} \quad (3.44)
$$

$$
\phi_{rl,2} = \arctan \frac{-2m^2 k_n \sqrt{k_{\tau}^2 - k'^2}}{m^4 k_n^2 - k_{\tau}^2 + k'^2} \quad (3.45)
$$

The total phase of a ray is then

$$
\phi = \phi_i + \phi_{op} + \phi_{fl} + \phi_{rl,X} \quad (3.46)
$$
### 3.2.3 Total scattered field and intensity

#### Amplitude of each ray

Finally we can get amplitude of each emergent ray $A_{p,X}$. As defined in Eq. (3.10), the amplitude $A_{p,X}$ of each emergent ray is the product of the Fresnel factor $\epsilon_{p,X}$ and the square root of the divergence factor $D_p$. For example, the amplitude of the external reflected ray ($p=0$) is

$$ 
\begin{aligned}
A_{0,X} &= |\epsilon_{0,X}|\sqrt{D_0} \\
&= |r_{0,X}|\sqrt{\frac{R''_{11} R''_{21}}{(r + R''_{11})(r + R''_{21})}}
\end{aligned} 
\qquad (3.47) $$

and for the refracted ray of $p=1$

$$ 
\begin{aligned}
A_{1,X} &= |\epsilon_{1,X}|\sqrt{D_1} \\
&= |t_{0,X} t'_{1,X}|\sqrt{\frac{R'_{11} R'_{21}}{R_{12} R_{22}} \cdot \frac{R'_{12} R'_{22}}{(r + R'_{12})(r + R'_{22})}}
\end{aligned} 
\qquad (3.48) $$

where $t_{0,X}$ and $t'_{1,X}$ are respectively the Fresnel refraction coefficients at the first and second interaction points.

#### Total scattered field

The total scattered field at a given angle is the superposition of the complex amplitudes $S_{p,X}$ (Eq. (3.1)) of all the rays arriving in that direction

$$ S_X = \sum_{p=0}^{\infty} S_{p,X} \qquad (3.49) $$

In VCRM, the directions of emergent rays are determined by the angles and the positions of the incident rays, the orders of the emergent rays and the refractive index of the particle. Their emergent angles have irregular intervals. As an example, Fig. 3.2 shows the variation of the amplitude as function of the emergent rays. It is clear that we can't calculate the summation of the amplitude of the rays of orders $p=0$ and 3 in a given angle since we do not know the phase and the amplitude of the emergent rays at that angle. For example at 128.8° as indicated in Fig. 3.2, a careful interpolation of the amplitudes and the phase shifts (not shown here) are necessary to calculate the summation of the complex amplitude. It should be noted that several rays of the same order may arrive at the same angle. Here the rays of order $p=3$ arrive twice in the angles less than the rainbow angle (129.1° here). Usually, the
<figure><img src="image_14.png" /><figcaption>Figure 3.2: Amplitude of scattered light near the second rainbow angle for a plane wave $\lambda = 0.6328$ µm scattered by a sphere of refractive index $m = 1.333$ and radius $a = 50$ µm</figcaption></figure>

interpolation algorithm applies to the monotone interval, therefore we must split the data obtained by ray tracing into groups in each of which the validation of the angle is monotone, i.e. only increasing or only decreasing, but never the both. Then we can calculate the amplitudes and the phases of all the rays arriving at the same angle by interpolation. And the summation of them gives rise to the complex amplitude of the total field. The intensity is just the square of the module of the total complex amplitude.

The scattered wave in the far field has the character of spherical wave. The direction of scattering is defined by the deviation angle $\theta$ with respect to the direction of propagation of the incident light and an azimuth angle $\varphi$. Let $I_0$ be the intensity of the incident light, $I$ the intensity of the scattered light in a point at a large distance $r$ from the particle. Since $I$ must be proportional to $I_0$ and $r^{-2}$, we may write

$$I = \frac{I_0 i(\theta, \varphi)}{k^2 r^2} \qquad (3.50)$$

Here the square of wave number $k$ is added to make $i(\theta, \varphi)$ a dimensionless function of the direction but independent of $r$. It also depends on the orientation of the particle with respect to the incident wave and the polarization state of the incident wave. The relative values of $i$, may be plotted in a polar diagram, as a function of $\theta$ in a fixed plane (given $\varphi$) respect to the direction of incidence. This diagram is called scattering diagram of the particle. So to get the final scattering diagram, our results of intensities should be still multiplied by $k^2 r^2$.
<figure><img src="image_15.png" /><figcaption>Figure 3.3: Scattering diagram predicted by VCRM and LMT for a sphere of radius $a = 50~\mu\text{m}$ illuminated by a plane wave of wavelength $\lambda = 0.6328~\mu\text{m}$. Refractive index of particle $m = 1.333$. The results of VCRM are shifted by $10^2$ for clarity.</figcaption></figure>

### 3.2.4 Numerical results and comparisons

In this thesis, I have written a code of VCRM in C# for the scattering of a plane wave by an elliptical particle in a symmetric plane. In this subsection, we focus on the validation of this code by comparing the scattering diagrams of the spherical particles calculated by VCRM with the results of Lorenz-Mie theory (LMT), which provides rigorous solution for light scattering by spherical particles. Then in chapter 4 this code is applied to the light scattering by large spheroidal particles to show the influence of particle non-sphericity on rainbow patterns.

In the following numerical calculations, the wavelength $\lambda$ of the incident wave is fixed to $0.6328~\mu\text{m}$. The particle is embedded in vacuum and its refractive index $m$ is set to be $1.333$. The maximum order of ray is 5 because the calculation contribution of higher orders are very small [114]. And numbers of incident and scattered rays are respectively 20000 and 10000. Besides, diffraction, which is very important in the forward direction, is also taken into account.

In principle, the ray model is only valid for the particles of size much larger than the wavelength. Therefore we start our comparisons for a water droplet of radius $a = 50~\mu\text{m}$. Fig. 3.3 illuminates the scattering diagram predicted by VCRM and LMT in perpendicular and parallel polarizations. From the figures we can find that VCRM is in agreement with LMT especially for perpendicular polarization. However discrepancies exist in the vicinity of $90^\circ$ and Alexander's dark region.

To further examine the effect of particle size on the precision of ray model, Fig. 3.4 shows the scattering diagram calculated by VCRM and LMT in perpendicular polarization for a sphere of radius $a = 30~\mu\text{m}$. We can find that the general agreement
<figure><img src="image_16.png" /><figcaption>Figure 3.4: Same parameter setting as in Fig. 3.3(a) except particle radius $a = 30$ µm.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 3.5: Same parameter setting as in Fig. 3.3(a) except particle radius $a = 10$ µm.</figcaption></figure>

between VCRM and LMT is still very good. However discrepancies near rainbow angle and 90° is more remarkable compared with Fig. 3.3.

In Fig. 3.5, we further reduce radius of the spherical particle to 10 µm. The discrepancies between VCRM and LMT become more significant. Even if the profiles of their scattering diagrams are still similar, details are different in the region larger than 70°. But even in this case, the overall accuracy of VCRM is still acceptable.

The scattering properties of non-spherical particles can differ dramatically from those of spheres. Now we apply VCRM to light scattering by prolate and oblate spheroids to study the influence of particle non-sphericity on the rainbow pattern. The scattering diagrams of spheroidal particles are compared with that of a sphere.

Consider a spheroid of refractive index $m$ illuminated by a plane wave propagating in its equatorial plane. The equatorial radius of this spheroid is $a$ and the distance from center to pole along the symmetry axis is $b$. Then this spheroid can be expressed
<figure><img src="image_18.png" /><figcaption>Figure 3.6: Scattering diagrams computed by VCRM for a prolate spheroid ($b = 65$ µm and aspect ratio $κ = 1.3$) with that of a sphere (radius $a = 50$ µm). The plane wave of wavelength $λ = 0.6328$ µm propagates along $z$ axis polarized in $y$ axis. The observe plane is fixed in $xz$ plane. The results of VCRM are shifted by $10^2$ for clarity.</figcaption></figure>

as

$$ \frac{x^2 + z^2}{a^2} + \frac{y^2}{b^2} = 1 \qquad (3.51) $$

It is centered at the origin and symmetric to $y$ axis. The incident plane wave propagates in the $xz$ plane with an incident angle $θ_0$ respective to $z$ axis. In this thesis we are interested only in the scattering in this symmetric plane. We take the equatorial radius $a = 50$ µm and the radius in the transversal direction $b$ as parameter. If the aspect ratio $κ = \frac{b}{a}$ is unity, the spheroid is reduced to a sphere. If $κ > 1$, it is a prolate spheroid. If $κ < 1$ it is an oblate spheroid.

Fig. 3.6 is a comparison of the scattering diagrams of a sphere (radius $a = 50$ µm) with that of a prolate spheroid ($b = 65$ µm and $κ = 1.3$). From the figure we can find that differences exist in the vicinity of the first and second rainbow angles near 139° and 127°. Compared with the scattering diagram of a sphere, the intensity of the second order rainbow of the prolate spheroid increases and becomes very close to that of the first order rainbow. This can be explained by the convergence effect due to the surface curvature of the particle in the vertical direction. In fact, the rays of the first order rainbow experience one internal reflection while the second ones experience two internal reflections. Therefore, the focalization effect decreases faster in the second rainbow than the first rainbow as the major axis increases. In next chapter, we will quantify the relation between the intensity ratio of the first to the second order rainbow and the curvature radii in the equatorial plane. Besides, a remarkable peak appears near 116°. Its formation will be examined in the next chapter.

Fig. 3.7 is a comparison of the scattering diagram of a sphere (radius $a = 50$
<figure><img src="image_19.png" /><figcaption>Figure 3.7: Scattering diagrams computed by VCRM for an oblate spheroid ($b = 45$ µm and aspect ratio $κ = 0.9$) with that of a sphere (radius $a = 50$ µm). The plane wave of wavelength $λ = 0.6328$ µm propagates along $z$ axis polarized in $y$ axis. The observe plane is fixed in $xz$ plane. The results of VCRM are shifted by $10^2$ for clarity.</figcaption></figure>

µm) with that of an oblate spheroid ($b = 45$ µm and $κ = 0.9$). From the figure we can find is a peak near 158.2° for the scattering diagram of oblate spheroid. This is our forementioned hyperbolic umbilic (HU) catastrophe caused by the eccentricity of oblate spheroid. Nye explained that this is the interference of the two light rays in the horizontal (equatorial) plane and two skew light rays which do not lie in the horizontal plane on entering the droplet, but are confined to the horizontal direction upon exiting the droplet [83]. These phenomena can be explained easily in the framework of VCRM by the analytical convergence factor. We will discuss it in detail in chapter 4.
# Chapter 4

# Rainbow of a spheroidal particle

## 4.1 Revisit of Airy theory

Airy theory (see section 2.3) can predict the intensity distribution near rainbow angle and offers a greatly simplified computation compared to rigorous theories, especially for large particles. However as an approximate method, its precision is limited. Here we will revisit its theoretical and numerical aspects, and examine its precision, especially that related to the particle metrology, by comparison with Debye series and VCRM.

### 4.1.1 Comparison of Airy theory with Debye theory

First of all, we will compare the scattering diagrams calculated by Airy theory with those obtained by the Debye series expansion – a rigorous theory. Fig. 4.1 illustrates the scattered intensity distribution around the first and the second order rainbows predicted by the two theories for a sphere of radius $a = 50 \ \mu m$ and refractive index $m = 1.333$ illuminated by a plane wave of wavelength $\lambda = 0.6328 \ \mu m$. A significant discrepancy can be observed for both the positions and the maximum intensities in supernumerary bows. This is due to the two approximations introduced in Airy theory.

The first is in the calculation of the phase of rays. Airy deduced the cubic function to describe the variation of the phase of the emergent rays as function of the coordinate $u$ perpendicular to the geometrical rainbow ray, by using the Taylor expansion of the deviation of angle of the rays $\theta$ in the neighborhood of the geometrical rainbow angles $\theta_{rg}$. This means that the cubic function is only valid in the neighborhood of geometrical rainbow angles. However, in order to obtain the amplitude of the scattered field, Airy extended the integration on the variable $u$ from $-\infty$ to $+\infty$.

The second approximation is on the amplitude of the emergent rays. To simplify
the integration, Airy supposed further that all the emergent rays have the same constant amplitude. In fact, we know well that the amplitude of the emergent ray is certainly not constant, especially, the amplitude of the rays in geometrical rainbow angle tends to infinity.

<figure><img src="image_20.png" /><figcaption>Figure 4.1: The comparison of main peak and supernumerary bows calculated by Debye series and Airy theory for a sphere of radius $a = 50 \ \mu m$ and refractive index $m = 1.333$.</figcaption></figure>

In VCRM, instead of assuming a cubic function, the phase of each emergent ray is calculated numerically. The amplitude is also evaluated directly by Fresnel formulas and divergence factor. Therefore, VCRM should be theoretically more accurate than Airy theory. In the following subsection, we will compare the predictions of Airy theory with those of Debye series and VCRM, then quantify the errors introduced in the particle measurement.

### 4.1.2 Comparison between Airy theory and VCRM

We examine firstly the scattering patterns near rainbow angles computed by VCRM and Airy theory for the case presented in the previous subsection.

The scattering diagrams near the first and the second order rainbows computed by VCRM, Airy theory and Debye series expansion are shown respectively in Fig. 4.2(a) and Fig. 4.2(b). In the calculation of VCRM, the number of incident rays and scattered rays are set to 20000 and 10000 respectively. From Fig. 4.2(a), we find that the first order rainbow predicted by VCRM is in excellent agreement with Debye series except in the vicinity of geometrical rainbow angle, while the difference between Airy theory and Debye theory increases as the order of the supernumerary bows. In the second order rainbow, the scattered intensity of VCRM is slightly smaller than that of Debye, but their peak positions are in good agreement. We constate greater
differences between the Airy theory and the Debye theory for both the intensity and the peak position. For example, the angular position differences between Airy theory and Debye theory for the third supernumerary bow is as large as 0.5° for in the first order rainbow and near 1° in the second order rainbow.

<figure><img src="image_21.png" /><figcaption>Figure 4.2: First and second order rainbows predicted by Debye theory, VCRM and Airy theory for a sphere whose radius $a = 50 \ \mu m$ and refractive index $m = 1.333$.</figcaption></figure>

For larger particles, the discrepancy between the Airy theory and the Debye theory is smaller, but the VCRM predicts always better the scattering patterns than the Airy theory. Fig. 4.3 and Fig. 4.4 show the same case as Fig. 4.2 but the particle radius are 100 $\mu m$ and 500 $\mu m$ respectively. The discrepancies between Airy theory and Debye series, i.e. the harmonic shift of the angular position of the supernumerary bows and the overestimation of the scattered intensity, are still remarkable, especially for the second order rainbow. With increase in particle size, both VCRM and Airy theory fit Debye better in angular position. And VCRM is still more accurate than Airy theory respect to Debye series. Besides we can also observe that the supernumerary structure depends on the particle size. As we can find in Eq. (2.86), the greater the particle size, the smaller the interval between two adjacent bows.

The discrepancy of Airy theory and VCRM from Debye theory depends also on the refractive index of the particle. Fig. 4.5 and Fig. 4.6 are the same case as Fig. 4.2 except that the refractive index of the particle $m$ is equal to 1.2 and 1.5 respectively. We find that VCRM is in good agreement for the first and the second rainbow of the particle with refractive index $m = 1.5$ and the first rainbow of the particle with refractive index $m = 1.2$. The discrepancy of the Airy theory from Debye theory is always much more significant than VCRM in these three cases.

In the case of the second order rainbow for a particle of refractive index $m = 1.2$ (Fig. 4.5 (b)). The scattered intensity predicted by the Airy theory is much greater than that of the Debye theory and the difference on the positions of the supernumerary
<figure><img src="image_22.png" /><figcaption>Figure 4.3: Same parameter setting as in Fig. 4.2 except particle radius $a = 100~\mu m$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 4.4: Same parameter setting as in Fig. 4.2 except particle radius $a = 500~\mu m$.</figcaption></figure>

bows is also significant. VCRM predicts well a portion of the scattering pattern in the main peak (between 164° and 170°) and a good angular position for the first supernumerary bow. But no higher orders in the supernumerary bows positions.

We can conclude that both Airy theory and VCRM predict better scattering patterns for particles of larger size and bigger refractive index. And in most cases, VCRM agrees much better with Debye theory than the Airy theory.

### 4.1.3 Discrepancy on the particle sizing

We will now investigate the discrepancy in the particle sizing caused by using the supernumerary bow positions of Airy theory. To this end, we need to extract the
<figure><img src="image_24.png" /><figcaption>Figure 4.5: Same parameter setting as in Fig. 4.2 except refractive index $m = 1.2$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 4.6: Same parameter setting as in Fig. 4.2 except refractive index $m = 1.5$.</figcaption></figure>
angular positions of the main peak and the supernumerary bows in the scattering patterns.

The peak positions of all supernumerary bows in Airy theory can be calculated analytically by using the Eqs. (2.85) and (2.86) given in section 2.4. The peak positions of the supernumerary bows of VCRM and Debye theory are to be extracted numerically from their scattering diagrams. If the angular step in the calculation is sufficiently small, the peak positions can be obtained directly by finding the local maximum in the scattering diagram. To obtain a good precision on the peak position, we must calculate for a huge number of rays and it will be very time consuming.

An alternative method is to calculate scattering diagrams with reasonably larger number of rays and then use the interpolation method to find the peaks of main and secondary supernumerary bows. The procedure is as follows.

We give a staring point $\theta_i$, usually the geometrical rainbow angles and calculate the intensity difference at two successive points, i.e. $\delta_1 = I(\theta_{i-1}) - I(\theta_i)$ and $\delta_2 = I(\theta_i) - I(\theta_{i+1})$. Then we recalculate these two differences by increasing (for the first order rainbow) or decreasing (for the second order rainbow) the scattering angle $i$ until the two differences have an opposite sign, $\theta_i$ will be the approximate peak position.

Then we suppose that the intensity around the local peak can be expressed as a quadratic function of scattering angle:

$$I(\theta) = a'\theta^2 + b'\theta + c' \qquad (4.1)$$

By using the standard routine for curve fitting in C# and the 7 points around $i$, i.e. $(\theta_{i-3}, \dots, \theta_{i+3})$, we can find the coefficients $a'$, $b'$ and $c'$. The angular position of peak $\theta_{peak}$ is then given by

$$\theta_{peak} = -\frac{b'}{2a'} \qquad (4.2)$$

The above procedure is applicable to all the supernumerary bows calculated by Debye theory and VCRM except the main peak of VCRM for which a special method must be used.

By observing the structure of the main peak in VCRM scattering patterns, we can note that there is no peak in the main bow and the intensity tends to infinity at geometrical rainbow angle, this is to be remedied by including the diffraction effect. But to be able to compare with Airy theory we intend to find its “possible peak” position which is located around the point where the second derivative is zero.

To find this point, we suppose that the scattered intensity is a cubic function of the scattering angle:

$$I(\theta) = a'\theta^3 + b'\theta^2 + c'\theta + d' \qquad (4.3)$$
Knowing the fact that the second derivative changes the sign around the point whose second derivative is zero, we calculate firstly numerically the difference of the first derivative of the curve at a given angle $\theta_j$

$$ \Delta I'(\theta_j) = \frac{I(\theta_{j+1}) - I(\theta_j)}{\theta_{j+1} - \theta_j} - \frac{I(\theta_j) - I(\theta_{j-1})}{\theta_j - \theta_{j-1}} \quad (4.4) $$

which has the same sign as the second derivative at this point $\theta_j$. Then we compare the value of $\Delta I'(\theta_j)$ at two successive points, if the signs are different then the pont to be found is around these two points $j - 1$ and $j$. Similar as for the other bows described above we choose 6 points from $j - 3$ to $j + 2$ and use the standard curve fitting routine in C# to find the coefficients in Eq. (4.3), The angular position of the main peak $\theta_{mp}$ can be calculated by

$$ \theta_{mp} = -\frac{b'}{3a'} \quad (4.5) $$

On the other hand, we can also deduce the VCRM main peak position by using the subsidiary supernumerary bow positions of VCRM found above and the relation between the main peak and subsidiary supernumerary bows given in the Airy theory. From Eqs. (2.87) and (2.88), I get the relation between angular position of main peak $\theta_{Airy,0}$ and those of supernumerary bows $\theta_{Airy,K}$ ($K = 1, 2, 3...$) in Airy theory :

$$ \theta_{VCRM,0} = \frac{\theta_{VCRM,K} - \theta_{go}(p, m)(1 - c)}{c} \quad (4.6) $$

where $\theta_{go}(p, m)$ is the geometrical optics rainbow angle. Assume that this angular relation also works for VCRM, then we can deduce the main peak position $\theta_{VCRM,0}$ in VCRM from positions of its supernumerary bows $\theta_{VCRM,K}$.

We have compared the main peak position of VCRM obtained by the two methods described above and found that the main peaks found directly by the second derivative of scattering patterns of VCRM have the same trend as Debye theory but the difference is greater than those found using Airy formula. So in what follows, we use the results calculated with Airy theory.

Fig. 4.7 shows the angular positions of the main peak $K = 0$ as well as subsidiary supernumerary bows of orders $K = 1, 2$ calculated by Debye, VCRM and Airy as function of particle radius for a sphere of refractive index $m = 1.333$ illuminated by a plane wave of wavelength $\lambda = 0.6328 \ \mu m$. We find that all supernumerary bows tend to the geometrical optics rainbow angle for large particle and the results of Airy theory and VCRM are in rather good agreement with Debye theory. The peak positions of VCRM fit always better with those of Debye than Airy theory.

To illustrate the discrepancy of the angular positions of Airy theory and VCRM from those of Debye theory, we show in Fig. 4.8 the respective angular differences
<figure><img src="image_26.png" /><figcaption>Figure 4.7: Angular positions of main peak and subsidiary supernumerary bows predicted by Airy theory, Debye and VCRM as function of particle radius. The rainbow angle of GO is also added for comparison.</figcaption></figure>

calculated from the results in the figure above. It is clear that the angular differences between VCRM and Debye is much smaller than those between Airy theory and Debye theory, especially in the first order rainbow.

Fig. 4.9 shows the peak positions differences of the Airy theory and VCRM from Debye theory as function of the refractive index $m$ for particle radius $a = 50 \ \mu m$. These differences decreases with the refractive index.

<figure><img src="image_27.png" /><figcaption>Figure 4.8: Angular differences of main peak and supernumerary bows between VCRM and Debye, and those between Airy and Debye, as function of particle radius.</figcaption></figure>

The discrepancy in the angular positions of supernumerary bows introduces an error in the particle sizing. For example, the Airy theory predicts bigger angle po-
<figure><img src="image_28.png" /><figcaption>Figure 4.9: Angular differences of main peak and supernumerary bows between VCRM and Debye, and those between Airy and Debye, as function of refractive index for a droplet whose radius $a = 50 \ \mu m$.</figcaption></figure>

sitions of bows than the Debye theory in the first order rainbow, which correspond to a smaller particle. In principle, this error in particle sizing can be evaluated by looking for the scattering diagram of Debye theory which corresponds to the angular peak positions of Airy theory or VCRM. This concerns an inverse problem and relatively difficulty to realized. In order to estimate roughly the error introduced due to the discrepancy of the peak position of the Airy theory and VCRM from the Debye theory, we use in the following the analytical relation between the peak positions and particle radius in the Airy theory.

We rewrite the size parameter in Eq. (2.85) for the main peak ($K = 0$) as an explicit function of particle radius:

$$|\theta_0 - \theta_{rg}(p, m)| = 1.087376 \left( \frac{\lambda^2 h}{48a^2} \right)^{1/3} \qquad (4.7)$$

For the supernumerary bows of order $K > 0$, we rewrite Eq. (2.86) as function of particle radius:

$$|\theta_K - \theta_{rg}(p, m)| = \left[ \frac{9\lambda^2}{16a^2} (K + \frac{1}{4})^2 h \right]^{1/3} \qquad (4.8)$$

$h$ can be calculated by the order of rainbow $p$ and refractive index $m$ according to Eq. (2.74).

For first order rainbow ($p = 2$), from Fig. 2.6, we know that the angular positions of Airy Maxima is larger than rainbow angle in the geometrical optic. Therefore the differentiation of Eqs. (4.7) and (4.8) with respect to the particle radius yields

$$\Delta\theta_0 = -\frac{1.087376}{3} \left( \frac{\lambda^2 h_{p=2}}{6} \right)^{1/3} a^{-5/3} \Delta a \qquad (4.9)$$
for Airy main peak, and for other supernumerary bows ($K > 0$)

$$ \Delta\theta_K = -\frac{2}{3} \left[ \frac{9\lambda^2}{16} (K + \frac{1}{4})^2 h_{p=2} \right]^{1/3} a^{-5/3} \Delta a \quad (4.10) $$

where

$$ h_{p=2} = \frac{9(4 - m^2)^{1/2}}{4(m^2 - 1)^{3/2}} \quad (4.11) $$

Similarly, for second order rainbow ($p = 3$), we have

$$ \Delta\theta_0 = \frac{1.087376}{3} \left( \frac{\lambda^2 h_{p=3}}{6} \right)^{1/3} a^{-5/3} \Delta a \quad (4.12) $$

for Airy main peak, and for other supernumerary bows ($K > 0$)

$$ \Delta\theta_K = \frac{2}{3} \left[ \frac{9\lambda^2}{16} (K + \frac{1}{4})^2 h_{p=3} \right]^{1/3} a^{-5/3} \Delta a \quad (4.13) $$

where

$$ h_{p=3} = \frac{64(9 - m^2)^{1/2}}{9(m^2 - 1)^{3/2}} \quad (4.14) $$

From these equations we can find that the discrepancy of the angular position of peak will introduces an error $\Delta a$ in the particle radius measurement.

Fig. 4.10 shows the error introduced in the particle sizing by the discrepancy of peak positions differences of the Airy theory and VCRM from Debye theory as function of the particle radius for the refractive index $m = 1.33$. We find that in the first order rainbow, the relative error in the particle sizing is almost constant and equal to about 5% for the Airy theory and 0.5% for VCRM. While in the second order rainbow, it is the absolute error is almost constant. The particle radius found by Airy theory from the first and the second supernumerary bow positions $K = 1, 2$, is about 30 µm larger than the given size. The particle radius determined by the second VCRM supernumerary bow position $K = 2$ is in good agreement with given size and that found from the first VCRM supernumerary bow position is smaller than the given size. The particle size found from the main peak is larger than the given one for both Airy theory and VCRM, the relative error is about 200% for Airy theory and 150% for VCRM for a particle of 50 µm. The relative errors are 60% and 40% for a particle of 500 µm.

## 4.2 Rainbows of spheroidal particles

Scattering diagrams of non-spherical particles differ dramatically from those of spheres. We cite in the following the research on two typical non-spherical particles, one can be
<figure><img src="image_29.png" /><figcaption>Figure 4.10: Size error introduced by Airy theory and VCRM respect to Debye series, as function of particle radius.</figcaption></figure>

considered as a prolate spheroidal particle and the other is an acoustically levitated oblate spheroidal particle.

In an experimental study of rainbow formed by pendant drops in 1978, Sassen observed that the relative intensity of different orders of rainbow depends on the deformation degree of drops and noted [92] “near-spherical particles can generate not only exceptionally well-defined primary, secondary, and associated supernumerary rainbows, but unexpectedly strong higher order bows as well”. He intended to explain the phenomena concerning the first and second order rainbows, “perhaps due to the generally reduce near-backscattering efficiency of the pendant drops, the relative strength of the secondary bow in comparison with that of the primary bow typically increases with increasing drop size, often resulting in secondary bows of equal or greater strength for the largest drops studied” [92] (page 1088, paragraph 2 in right).

Later in 1998, Langley and Marston, in their study of rainbow by slightly oblate drops, noted that “Sassen found that large pendant drops showed increased rainbow scattering, especially for the higher-order bows, because of decreased curvature in the equatorial plane, and suggested that large raindrops having reduced curvature at their equators would also show an increased Airy brightness. The mechanism described in our study, however, shows that vertically focused rays outside the equatorial plane may contribute significantly to rainbow brightness for smaller drops that are oblate but have their greatest curvatures at the equator” [93] (page 1525, paragraph 2 in right, line 7-18).

We would firstly comment that the particles in their study are quite different. The drop levitated of Langley and Marston is oblate spheroidal particle while the pendant drop of Sassen can be considered as a prolate spheroid. Their convergence effect can
be very different. Secondly, Langley and Marston talked about the size of particle, but the characteristics of rainbow of particles much larger than the wavelength of incident wave depend much more on the shape of the particle than in their size. We note both of them intend to conjecture the mechanism according to their observations, however, none of their explanations is convincing.

In the framework of VCRM, the amplitude of the scattered field, the phase and the direction of all the emergent rays in the symmetric plane of a spheroidal particle can be calculated analytically. So we can predict the position of the supernumerary bows easily. In this section, we will give a detailed description of this method and intend to quantify the relation of the intensity ratio of the first order rainbow to the second order rainbow with the curvature radius of the particle surface in the direction perpendicular to the scattering plane, i.e. equatorial plane.

### 4.2.1 Amplitude of scattered field near rainbow angle

#### 1. Wavefront curvature radii of wavefront at each interaction

Consider a spheroid of refractive index $m$, centered at $O$ with $y$ as its symmetric axis, illuminated by a plane wave propagating in its equatorial plane, i.e. perpendicular to $y$ axis. The radius of this spheroid in this plane is $a$ and its half-axis in $y$ direction is $b$. Then this spheroid is described in a Cartesian coordinate system $Oxyz$ by Eq. (3.51).

The direction of the incident plane wave is along the $z$ axis. Here we are interested only by the scattering in this symmetric plane. By the symmetry of the problem, we know that all incident rays in this plane remain always in this plane. The two principle curvature radii can be evaluated step by step according to Eqs. (3.34) and (3.35). To facilitate the analysis in a familiar manner, we rewrite the two equations in the function of the incident angle $\alpha$ and the refraction angle $\beta$:

$$ \frac{k' \cos^2 \beta}{R'_{1,q}} = \frac{k \cos^2 \alpha}{R_{1,q}} + \frac{k' \cos \beta - k \cos \alpha}{\rho_1} \quad (4.15) $$

$$ \frac{k'}{R'_{2,q}} = \frac{k}{R_{2,q}} + \frac{k' \cos \beta - k \cos \alpha}{\rho_2} \quad (4.16) $$

Eq. (4.15) gives the relation of the curvature radii of wavefronts with the local surface curvature radius of the particle in the equatorial $xz$ plane while Eq. (4.16) is that for the direction perpendicular to the incident plane. Both $\alpha$ and $\beta$ are constant for the case under study. $R_{1,q}$ and $R_{2,q}$ ($q = 1, 2, \dots, p+1$) are the two principal wavefront curvature radii of rays before the $q^{th}$ interaction. $R'_{1,q}$ and $R'_{2,q}$ are those after the interaction (reflected or refracted). To ease the discussion, we use, from now on, $R'_{1,q}$ and $R'_{2,q}$ to designate the wavefront curvature radii of refracted rays, and $R''_{1,q}$ and $R''_{2,q}$ for the reflected rays.
The curvature radii of the particle surface and the wavefront may be positive or negative depending on the convex or concave property. In our convention, the curvature radii are positive if the surface is concave, i.e. the focal line/center is after the intersection point. Therefore, at the first interaction point, the two main curvature radii of the spheroid at equatorial plane (y = 0) are given by:

$$ \rho_1 = a \qquad (4.17) $$

$$ \rho_2 = \frac{b^2}{a} \qquad (4.18) $$

The wavefront curvature radii $R_{11}$ and $R_{21}$ of the plane wave are infinite, we can deduce the curvature radii $R''_{11}$ and $R''_{21}$ of reflected rays by Eqs. (4.15) and (4.16)

$$ R''_{11} = \frac{a \cos \alpha}{2} \qquad (4.19) $$

$$ R''_{21} = \frac{b^2}{2a \cos \alpha} \qquad (4.20) $$

The positive sign indicates that the reflected wave is divergent in the two directions: perpendicular and parallel to the incident plane, i.e. the curvature center of reflected wave is on the other side of emergent wave.

The curvature radii $R'_{11}$ and $R'_{21}$ of rays refracted into spheroid are

$$ R'_{11} = \frac{ma \cos^2 \beta}{m \cos \beta - \cos \alpha} = \frac{ma^2 \cos^2 \beta}{h} \qquad (4.21) $$

$$ R'_{21} = \frac{mb^2}{a(m \cos \beta - \cos \alpha)} = \frac{mb^2}{h} \qquad (4.22) $$

where $h = a(m \cos \beta - \cos \alpha)$. Here the positive values of $R'_{11}$ and $R'_{21}$ correspond well the fact that the wave refracted into the particle is convergent also in the two directions.

The distance between two successive interaction points in the particle is $d = 2a \cos \beta$, therefore the wavefront curvature radii of the ray arriving at the second interaction point are given by

$$ \begin{aligned} R_{12} &= R'_{11} - d \\ &= -\frac{a \cos \beta (m \cos \beta - 2 \cos \alpha)}{m \cos \beta - \cos \alpha} \\ &= -\frac{a \cos \beta (h - a \cos \alpha)}{h} \end{aligned} \qquad (4.23) $$

$$ \begin{aligned} R_{22} &= R'_{21} - d \\ &= \frac{mb^2}{h} - 2a \cos \beta \\ &= -\frac{2ah \cos^2 \beta - mb^2}{h} \end{aligned} \qquad (4.24) $$
Then by using again Eqs. (4.15) and (4.16), we can obtain the wavefront curvature radii of emergent ray of order $p = 1$

$$ 
\begin{aligned}
R'_{12} &= -\frac{a \cos \alpha (m \cos \beta - 2 \cos \alpha)}{2(m \cos \beta - \cos \alpha)} \\
&= -\frac{a \cos \alpha (h - a \cos \alpha)}{2h}
\end{aligned}
\qquad (4.25) $$

$$ 
\begin{aligned}
R'_{22} &= -\frac{b^2(2ah \cos \beta - mb^2)}{2h(mb^2 - ah \cos \beta)} \\
&= \frac{b^2}{2h} - \frac{ab^2 \cos \beta}{2(mb^2 - ah \cos \beta)}
\end{aligned}
\qquad (4.26) $$

And for the first internal reflected ray, its wavefront curvature radii $R''_{12}$ and $R''_{22}$ are calculated by

$$ 
\begin{aligned}
R''_{12} &= \frac{a \cos \beta (m \cos \beta - 2 \cos \alpha)}{m \cos \beta - 3 \cos \alpha} \\
&= a \cos \beta + \frac{a \cos \alpha \cos \beta}{m \cos \beta - 3 \cos \alpha}
\end{aligned}
\qquad (4.27) $$

$$ 
\begin{aligned}
R''_{22} &= \frac{(2ah \cos \beta - mb^2)b^2}{2a \cos \beta(2ah \cos \beta - mb^2) - b^2h} \\
&= \frac{Hb^2}{2aH \cos \beta - b^2h}
\end{aligned}
\qquad (4.28) $$

where $H = 2ah \cos \beta - mb^2 = 2a^2m \cos^2 \beta - 2a^2 \cos \alpha \cos \beta - mb^2$.

Then using the same principle, the wavefront curvature radii of the ray arriving at the third interaction point are given by

$$ 
\begin{aligned}
R_{13} &= R''_{12} - d \\
&= \frac{a \cos \alpha \cos \beta}{m \cos \beta - 3 \cos \alpha} - a \cos \beta \\
&= -\frac{a \cos \beta (m \cos \beta - 4 \cos \alpha)}{m \cos \beta - 3 \cos \alpha}
\end{aligned}
\qquad (4.29) $$

$$ 
\begin{aligned}
R_{23} &= R''_{22} - d \\
&= -\frac{(4a^2 \cos^2 \beta - b^2)(2ah \cos \beta - mb^2) - 2ab^2h \cos \beta}{2a \cos \beta(2ah \cos \beta - mb^2) - b^2h} \\
&= -\frac{4a^2 \cos^2 \beta H - 2ab^2h \cos \beta - Hb^2}{2aH \cos \beta - b^2h}
\end{aligned}
\qquad (4.30) $$
The wave curvature radii of emergent ray of order $p = 2$

$$ 
\begin{aligned}
R'_{13} &= -\frac{a \cos \alpha (m \cos \beta - 4 \cos \alpha)}{2(m \cos \beta - 2 \cos \alpha)} \\
&= -\frac{a \cos \alpha (h - 3a \cos \alpha)}{2(h - a \cos \alpha)} \\
&= -\frac{a^2 \cos \alpha (m \cos \beta - 4 \cos \alpha)}{2(h - a \cos \alpha)}
\end{aligned} 
\quad (4.31) $$

$$ 
\begin{aligned}
R'_{23} &= -\frac{b^2[2a \cos \beta (2aH \cos \beta - b^2h) - Hb^2]}{2H(hb^2 - aH \cos \beta)} \\
&= -\frac{b^2[4a^2 \cos^2 \beta H - 2ab^2h \cos \beta - Hb^2]}{2H(hb^2 - aH \cos \beta)}
\end{aligned} 
\quad (4.32) $$

The wavefront curvature radii of the ray after second internal reflection can be calculated by

$$ R''_{13} = \frac{a \cos \beta (m \cos \beta - 4 \cos \alpha)}{m \cos \beta - 5 \cos \alpha} = a \cos \beta + \frac{a \cos \beta \cos \alpha}{m \cos \beta - 5 \cos \alpha} \quad (4.33) $$

$$ R''_{23} = \frac{b^2(4a^2H \cos^2 \beta - 2ab^2h \cos \beta - Hb^2)}{4a \cos \beta (2a^2H \cos^2 \beta - ab^2h \cos \beta - Hb^2) + hb^4} \quad (4.34) $$

The wavefront curvature radii of the ray arriving at the fourth interaction point are given by

$$ 
\begin{aligned}
R_{14} &= R''_{13} - d \\
&= \frac{a \cos \alpha \cos \beta}{m \cos \beta - 5 \cos \alpha} - a \cos \beta
\end{aligned} 
\quad (4.35) $$

$$ 
\begin{aligned}
R_{24} &= R''_{23} - d \\
&= \frac{b^2(4a^2H \cos^2 \beta - 2ab^2h \cos \beta - Hb^2)}{4a \cos \beta (2a^2H \cos^2 \beta - ab^2h \cos \beta - Hb^2) + hb^4} - 2a \cos \beta
\end{aligned} 
\quad (4.36) $$

The wavefront curvature radii of emergent ray of order $p = 3$

$$ R'_{14} = -\frac{a \cos \alpha (h - 5a \cos \alpha)}{2(h - 2a \cos \alpha)} \quad (4.37) $$

$$ R'_{24} = -\frac{b^2 [2(2a^2H \cos^2 \beta - ab^2h \cos \beta - Hb^2)(4a^2 \cos^2 \beta - b^2) + m^6]}{(2a^2H \cos^2 \beta - ab^2h \cos \beta - Hb^2) [4mab^2 \cos \beta - 2h(4a^2 \cos^2 \beta - b^2)]} \quad (4.38) $$

Similarly, we can get the wavefront curvature radii of other order rays.

## 2. Divergence factor of emergent rays for an spheroid

The divergence factor $D_p$ of an emergent ray after $q = p + 1$ time interactions with the dioptry surface is given in VCRM by Eq. (3.3). The divergence factor for the
reflected ray of order $p = 0$ is

$$D_0 = \frac{R''_{11} R''_{21}}{r^2} = \frac{b^2}{4r^2} \qquad (4.39)$$

The divergence factor for the emergent ray of order $p = 1$ is

$$ \begin{aligned} D_1 &= \frac{R'_{11} R'_{21}}{R_{12} R_{22}} \cdot \frac{R'_{12} R'_{22}}{r^2} \\ &= \frac{m^2 b^4 \cos \alpha \cos \beta}{4r^2 (m \cos \beta - \cos \alpha)^2 (mb^2 - ma^2 \cos^2 \beta + a^2 \cos \alpha \cos \beta)} \end{aligned} \qquad (4.40) $$

The divergence factor for the emergent ray of order $p = 2$ is

$$ \begin{aligned} D_2 &= D_1 \cdot \frac{R'_{13} R'_{23}}{R_{13} R_{23}} \\ &= \frac{m^2 b^6 \cos \alpha \cos \beta}{4r^2 (m \cos \beta - 2 \cos \alpha)(2ma^2 \cos^2 \beta - 2a^2 \cos \alpha \cos \beta - mb^2)} \times \\ & \quad \frac{1}{(2mb^2 \cos \beta - b^2 \cos \alpha - 2ma^2 \cos^3 \beta + 2a^2 \cos \alpha \cos^2 \beta)} \end{aligned} \qquad (4.41) $$

And the divergence factor for the emergent ray $p = 3$ can be calculated by

$$D_3 = D_2 \frac{R'_{14} R'_{24}}{R_{14} R_{24}} \qquad (4.42)$$

The divergence factor for emergent ray of higher order can be calculated by following recurrence formula

$$D_p = D_{p-1} \frac{R'_{1,p+1} R'_{2,p+1}}{R_{1,p+1} R_{2,p+1}} \qquad (4.43)$$

### 3. Divergence factor of emergent rays for a sphere

To check the divergence factor of a spheroid obtained above, we will apply them to a simple case where $a = b$, i.e. a sphere, for which we know the analytical formula. In the case, Eq. (4.39) is simplified to

$$D_{s,p=0} = \frac{a^2}{4r^2} \qquad (4.44)$$

The divergence factor for the emergent ray of order $p = 1$ is

$$D_{s,p=1} = \frac{m^2 a^2 \cos \alpha \cos \beta}{4r^2 (m \cos \beta - \cos \alpha)^2 (m \sin^2 \beta + \cos \alpha \cos \beta)} \qquad (4.45)$$
and that for emergent ray $p = 2$ is

$$D_{s,p=2} = \frac{m^2 a^2 \cos \alpha \cos \beta}{4r^2(m \cos \beta - 2 \cos \alpha)(2m \cos^2 \beta - 2 \cos \alpha \cos \beta - m)} \times \frac{1}{(2m \sin^2 \beta \cos \beta - \cos \alpha + 2 \cos \alpha \cos^2 \beta)} \quad (4.46)$$

On the other hand, the divergence factor $\mathcal{D}$ of order $p$ in GO is given by Eq. (2.51), recalled here for convention:

$$\mathcal{D}_p = \frac{\sin(2\alpha)}{4 \sin \theta \left( p \frac{\cos \alpha}{\sqrt{m^2 - \sin^2 \alpha}} - 1 \right)} \quad (4.47)$$

where $\theta$ is the scattering angle, which is related to the deviation angle $\theta_p$ by

$$\theta_p = q\theta - p\pi + 2n\pi \quad (4.48)$$

where $q = \pm 1$, $n$ is an integer so that $\theta$ is between $[0, \pi]$. $\theta_p$ is determined by the incident angle and the refraction angle by

$$\theta_p = 2p\beta - 2\alpha + \pi \quad (4.49)$$

By introducing Eqs. (4.48) and (4.49) into Eq. (4.47), the divergence factor of GO can be expressed as the incident angle and the refraction angle, which must be in agreement with those obtained by VCRM. Before the theoretical verification, it worth to note that the divergence factor of GO is independent of particle radius while that of VCRM includes the particle dimension (the two half axes $a$ and $b$ of the spheroid). When the latter is simplified to spherical case, the divergence factor of GO differs from that of VCRM by a factor $a^2/r^2$.

For the reflected ray of order $p = 0$, we have $\theta = \pi - 2\alpha$, so we get

$$\mathcal{D}_{p=0} = \frac{1}{4} \quad (4.50)$$

This is just $\frac{r^2}{a^2} \times D_{s,p=0}$, as it should be.

For the emergent ray of order $p = 1$, $\theta = 2(\beta - \alpha)$, we have

$$\begin{aligned} \mathcal{D}_{p=1} &= \frac{\sin(2\alpha)}{4 \sin \theta \left( \frac{\cos \alpha}{m \cos \beta} - 1 \right)} \\ &= \frac{\sin(2\alpha)}{4 \sin[2(\beta - \alpha)]} \cdot \frac{m \cos \beta}{\cos \alpha - m \cos \beta} \\ &= \frac{m^2 \cos \beta \cos \alpha}{4(\cos \alpha - m \cos \beta)^2 (\cos \beta \cos \alpha + m \sin^2 \beta)} \end{aligned} \quad (4.51)$$
This is equal to $\frac{r^2}{a^2} \times D_{s,p=1}$ as it should be.

For the emergent ray of order $p = 2$, $\theta = \pi - 2(2\beta - \alpha)$, then

$$ 
\begin{aligned}
\mathcal{D}_{p=2} &= \frac{\sin(2\alpha)}{4 \sin \theta \left(\frac{2 \cos \alpha}{m \cos \beta} - 1\right)} \\
&= \frac{\sin(2\alpha)}{4 \sin[2(2\beta - \alpha)]} \cdot \frac{m \cos \beta}{2 \cos \alpha - m \cos \beta} \\
&= \frac{m \sin \alpha \cos \alpha \cos \beta}{4(2 \cos \alpha - m \cos \beta) \sin(2\beta - \alpha) \cos(2\beta - \alpha)}
\end{aligned}
\quad (4.52) $$

This is also equal to $\frac{r^2}{a^2} \times D_{s,p=2}$.

From Eqs. (4.50)-(4.52) above, we can find that, for order $p = 0, 1, 2$, after multiplied a factor of $r^2/a^2$, our divergence factors derived in VCRM are identical to those defined in GO for a sphere. In principle, this verification can be done for high order rays but the mathematical calculation is very tedious. However, the numerical check can be made easily for any order.

We compile in Table 4.1 the divergence factors of orders $p = 3, 4$ calculated by VCRM for a sphere of refractive index $m = 1.33$ and radius $a = 50 \ \mu m$ according to

$$ \mathcal{D}_{s,p=3} = \mathcal{D}_{p=2} \cdot \frac{R''_{13} R''_{23}}{R'_{13} R'_{23}} \cdot \frac{R'_{14} R'_{24}}{R_{14} R_{24}} \quad (4.53) $$

and

$$ \mathcal{D}_{s,p=4} = \mathcal{D}_{s,p=3} \cdot \frac{R''_{14} R''_{24}}{R'_{14} R'_{24}} \cdot \frac{R'_{15} R'_{25}}{R_{15} R_{25}} \quad (4.54) $$

$\mathcal{D}_{p=2}$ is given as a starting point. The divergence factors calculated by GO is exactly the same as those of VCRM, so only one $\mathcal{D}_{s,p}$ value is given.

### 4.2.2 Peak intensity of supernumerary bows

In the scattering of light by a spheroidal particle in the equatorial plane, the phase of each ray is the same as that in the scattering by a spherical particle and the rainbows as well as the supernumerary bows are located in the same angles as those of a sphere. But the intensities of the bows of a spheroidal particle are different from those of a sphere due to the convergence caused by the curvature of the particle surface in the direction perpendicular to the equatorial plane. This curvature is constant for all rays, therefore we can calculate the peak intensity of each bows analytically. This section is devoted to the exploitation of the method.

In VCRM, the complex amplitude $S_{p,X}$ of an emergent ray is

$$ S_{p,X} = A_{p,X} \exp(i\Phi_p) \quad (4.55) $$
<table><caption>Table 4.1: Wavefront curvature radii and divergence factor of high order rays calculated VCRM for a sphere of refractive index $m = 1.33$ and radius $a = 50$ µm. The divergence factors are exactly the same as those obtained by Eq. (4.47) of GO.</caption><thead><tr><th>Incident Angle α (deg.)</th><th>α = 0</th><th>α = 30</th><th>α = 45</th><th>α = 60</th><th>α = 90</th></tr></thead><tbody><tr><td>Q<sub>p=2</sub></td><td>0.985</td><td>1.125</td><td>1.573</td><td>-34.384</td><td>-3.403</td></tr><tr><td>R'<sub>13</sub>R'<sub>23</sub></td><td>9926</td><td>8672</td><td>8349</td><td>-95094</td><td>0.001</td></tr><tr><td>R''<sub>13</sub>R''<sub>23</sub></td><td>1323</td><td>1441</td><td>1648</td><td>2052</td><td>3602</td></tr><tr><td>R<sub>14</sub>R<sub>24</sub></td><td>4048</td><td>2935</td><td>1622</td><td>-278</td><td>-1428</td></tr><tr><td>R'<sub>14</sub>R'<sub>24</sub></td><td>4887</td><td>3217</td><td>1516</td><td>-200</td><td>-9.042</td></tr><tr><td>Q<sub>s,p=3</sub></td><td>0.159</td><td>0.205</td><td>0.290</td><td>0.536</td><td>-9.443</td></tr><tr><td>R''<sub>14</sub>R'<sub>24</sub></td><td>1696</td><td>2254</td><td>299175</td><td>142</td><td>667</td></tr><tr><td>R<sub>15</sub>R<sub>25</sub></td><td>3460</td><td>1824</td><td>-430699</td><td>3244</td><td>1507</td></tr><tr><td>R'<sub>15</sub>R'<sub>25</sub></td><td>3900</td><td>1699</td><td>-5508</td><td>4444</td><td>4.0×10<sup>-4</sup></td></tr><tr><td>Q<sub>s,p=4</sub></td><td>0.062</td><td>0.134</td><td>0.732</td><td>-0.522</td><td>1.7×10<sup>-7</sup></td></tr></tbody></table>

where subscripts X designates the state of polarization, perpendicular or parallel to the plane of incidence. $\Phi_p$ is the phase of the $p^{th}$ order ray, and $A_{p,X}$ is its amplitude given by

$$A_{p,X} = E_i \sqrt{|D_p|} \epsilon_{p,X} \quad (4.56)$$

$\epsilon_{p,X}$ stands for the relative amplitude of the emergent ray and is determined by Fresnel coefficients.

In general, the total amplitude of scattered field at a given angle is the superposition of the complex rays of all orders arriving at that direction. However, if the aspect ratio of the spheroid is near unity (slightly deformed sphere), the low order rainbows as well as their supernumerary bows are separated far enough so that their intensity can be considered as formed by the two rays of the same order. Suppose that the incident angles of the two rays are respectively $\alpha_1$ and $\alpha_2$. From the formation mechanism of Airy structure we know that the phase difference between these two rays is multiple of $2\pi$. Therefore, their resultant amplitude at that point is equal to the sum of the amplitudes of the two rays:

$$I_{p,X}(\theta) = |A_{p,X}(\alpha_1) + A_{p,X}(\alpha_2)|^2 \quad (4.57)$$

Combined with expressions of divergence in subsection 4.2.1, this intensity can be calculated easily.

Then the intensity ratio $\chi_{p,X}^{p'}$ of two rainbows of orders $p'$ and $p$ can be calculated by

$$\begin{aligned} \chi_{p,X}^{p'} &= \frac{I_{p',X}(\theta')}{I_{p,X}(\theta)} \\ &= \left| \frac{A_{p',X}(\alpha'_1) + A_{p',X}(\alpha'_2)}{A_{p,X}(\alpha_1) + A_{p,X}(\alpha_2)} \right|^2 \end{aligned} \quad (4.58)$$
where $\theta$ and $\theta'$ are the main peak positions of the two rainbows. $\alpha_1$, $\alpha_2$ and $\alpha'_1$, $\alpha'_2$ are the incident angles of the corresponding rays.

### 4.2.3 Incident angles of rays arriving at bow peaks

In order to calculate the amplitudes of the two emergent rays arriving at the super-numerary bow peak according to the vertical and horizontal curvature of the particle surface, we need their incident angles. In VCRM these angles cannot be calculated directly. From subsection 4.1.2 we know that the angular positions of main peaks of rainbow predicted by VCRM and Airy theory are very close, so we can determine the incident angles of the rays arriving at the bow peak by using the analytical relation in the Airy theory.

The main peak position $\theta_{Airy,0}$ and the peak positions of the subsidiary supernu-merary bows ($K \ge 1$) are given by Eqs. (2.85) and (2.86) respectively. The rainbow angle in the geometrical optics is calculated by Eq. (2.64). For a given spheroidal particle ($m$ and $a$ are known) and given order $p$, the peak positions of the main and subsidiary bows can be solved from Eqs. (2.85) and (2.86). Then the corresponding incident angles $\alpha_1$ and $\alpha_2$ can be calculated by using Eqs. (4.48) and (4.49), as well as the Snell's law.

Here we take an example of a plane wave scattered by a spheroid of symmetrical semi-axis $a = 50 \ \mu m$ and refractive index $m = 1.33$. The angle position of the first order rainbow ($p = 2$) can be calculated by Eq. (2.64)

$$ \theta_{rg}(2, 1.33) = 137.5^{\circ} \qquad (4.59) $$

The main peak position is

$$ \theta_{Airy,0} = 139.0^{\circ} \qquad (4.60) $$

Since we know that the scattering angle is related to the incident angle $\alpha$ and the refracted angle $\beta$ by $\theta = \pi - 2(2\beta - \alpha)$, and $\beta = \arcsin(\sin \alpha/m)$, the two incident angles $\alpha_1$ and $\alpha_2$ corresponding to this scattering angle can be solved by following equation

$$ \theta_{Airy,0} = 2 \left[ 2 \arcsin \left( \frac{\sin \alpha}{m} \right) - \alpha \right] \qquad (4.61) $$

which give rise to $\alpha_1 = 50.85^{\circ}$ and $\alpha_2 = 67.45^{\circ}$.

This procedure can also be explained graphically. The solid curves in Fig. 4.11 shows main peak positions of different orders of rays calculated by Eqs. (4.48) and (4.49) as function of incident angle. The angular position of each main peak is indi-cated by a horizontal line. From the figure, we can find that there are two incident angles corresponding to the scattering angle for each main peak.
<figure><img src="image_30.png" /><figcaption>Figure 4.11: Scattering angle of rays $p = 2, 3, 4, 5$ as function of incident angle, and angular positions of main supernumerary peaks for a spheroid of refractive index $m = 1.33$) and symmetric semi-axis $a = 50~\mu m$.</figcaption></figure>

The main peak position as well as the two incident angles depend also on the size of the particle. Figs. 4.12-4.15 show respectively the scattering angle of rays $p = 2, 3, 4$ and $5$ as function of incident angle. The rainbow angles of geometric optics and the angular positions of main peaks of different order rainbows for spheroid of $m = 1.33$ with different semi-axis $a = 20, 50, 100$ and $1500~\mu m$ are indicated with dashed straight lines. For convenience of calculation, these data are also given in Tab. 4.2. It is evident that the difference between the two incident angles decreases as function of the particle size. This can be understood by the fact that to have the same phase difference the two incident rays must be farther for a smaller particle. Besides, the incident angles of each main peak can be got easily from figures which will be useful for our study of relative intensity as function of vertical and horizontal curvature radii in the next subsection.

### 4.2.4 Rainbow intensity ratio of a prolate spheroid

For the scattering in the symmetric plane of a spheroid, the rainbow angles as well as the angular positions of the supernumerary bows are the same as a spherical of the same size. But their intensities are affected by the curvature radius of the spheroid in the direction perpendicular to the scattering plane, which is related directly to the aspect ratio of the spheroid. So it is possible to calculate the intensity ratio between
<figure><img src="image_31.png" /><figcaption>Figure 4.12: Scattering angle of the first rainbow ($p = 2$) as function of incident angle, as well as its geometric optics rainbow angle and the angular positions of the main peak for spheroid/sphere of $m = 1.33$, $a = 20, 50, 500$ and $1500$ µm. The incident wavelength is $0.6328$ µm.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Figure 4.13: Same as Fig. 4.12 but for the second rainbow ($p = 3$).</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>Figure 4.14: Same as Fig. 4.12 but for the third rainbow ($p = 4$).</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure 4.15: Same as Fig. 4.12 but for the second rainbow ($p = 5$).</figcaption></figure>

different orders of rainbow.

When the refractive index and the two axes *a* and *b* of a spheroidal particle are given, the positions of Airy peaks can be calculated by Eqs. (2.85) and (2.86). Then
<table><caption>Table 4.2: Comparison of rainbow angle calculated by GO and angular position of Airy maxima by Airy theory for a water sphere ($m = 1.33$) of radius $a = 20, 50, 500$ and $1500$ µm. The incident wavelength is $0.6328$ µm.</caption><tr><td></td><td><i>m</i> = 1.33</td><td><i>p</i> = 2</td><td><i>p</i> = 3</td><td><i>p</i> = 4</td><td><i>p</i> = 5</td></tr><tr><td>Geometric Optics</td><td><i>&theta;</i><sub>go</sub>(<i>p</i>, <i>m</i>)</td><td>137.5</td><td>129.9</td><td>42.8</td><td>42.3</td></tr><tr><td rowspan="4">Airy Theory</td><td><i>a</i> = 1500</td><td>137.6</td><td>129.6</td><td>42.4</td><td>42.8</td></tr><tr><td><i>a</i> = 500</td><td>137.8</td><td>129.3</td><td>42.0</td><td>43.4</td></tr><tr><td><i>a</i> = 50</td><td>139.0</td><td>127.1</td><td>38.9</td><td>47.4</td></tr><tr><td><i>a</i> = 20</td><td>140.3</td><td>124.7</td><td>35.5</td><td>51.7</td></tr></table>

the corresponding incident angles $\alpha_1$ and $\alpha_2$ are determined for each order of rainbow according to the method described above. On this base, the divergence factor, the amplitude of each emergent ray as well as the intensity ratio of two orders of rainbow $p'$ and $p$ noted by $\chi_{p,X}^{p'}$ can be calculated according to the procedure described in Subsection 4.2.2.

It is worth to point out that the two incident angles corresponding to a given rainbow main peak is determined by the Airy theory, the intensity calculated analytically with these angles $I_{p,X}$ may differ from the scattering intensity $I_{p,X}^{VCRM}$ calculated numerically by the code VCRM. To investigate this discrepancy, we compare in Tab. 4.3 the intensities calculated by the two methods for a sphere. We find that the intensities calculated with the analytical expressions are in good agreement with those by the code VCRM. In the following, we will examine the intensity ratio of different orders of rainbow as function of the aspect ratio of a prolate spheroid.

<table><caption>Table 4.3: Comparison of the intensity calculated by the code VCRM and the analytical expressions for a sphere of refractive index $m = 1.33$ and radius $a = 50$ µm, and intensity ratio of Airy main peak of order $p = 3, 4$ and $5$ to order $2$. The incident wavelength is $\lambda = 0.6328$ µm and the polarization is perpendicular $X = 1$.</caption><thead><tr><th>p</th><th>θ<sub>Airy,0</sub></th><th>I<sup>VCRM</sup><sub>p,1</sub></th><th>α<sub>1</sub> [deg.]</th><th>α<sub>2</sub> [deg.]</th><th>I<sub>p,1</sub></th><th>I<sup>VCRM</sup><sub>p,1</sub>/I<sub>p,1</sub></th><th>χ<sup>p</sup><sub>2,1</sub></th></tr></thead><tbody><tr><td>2</td><td>139.0</td><td>148936</td><td>50.9</td><td>67.5</td><td>146433</td><td>1.02</td><td>-</td></tr><tr><td>3</td><td>127.1</td><td>13740</td><td>64.0</td><td>79.4</td><td>14137</td><td>0.97</td><td>0.097</td></tr><tr><td>4</td><td>38.9</td><td>3715</td><td>69.2</td><td>84.2</td><td>3853</td><td>0.96</td><td>0.026</td></tr><tr><td>5</td><td>47.4</td><td>758</td><td>72.1</td><td>87.1</td><td>756</td><td>1.00</td><td>0.005</td></tr></tbody></table>

Consider a spheroid illuminated by a plane wave of wavelength $\lambda = 0.6328$ µm and polarization is perpendicular to the symmetric axis (i.e. the scattering plane).

First of all, we reexamine in Fig. 4.16 the scattering diagrams of a sphere ($a = b = 50$ µm) with that of a spheroid ($a = 50$ µm and $b = 65$ µm. The refractive index of the particle is $m = 1.333$. We find that the intensity of second order of a spheroid increases and becomes very close to that of first order rainbow. Besides, a remarkable peak appears near $116^\circ$. It is due to the convergence effect of the rays of order $p = 3$.
<figure><img src="image_35.png" /><figcaption>Figure 4.16: Comparison of scattering diagrams for a sphere and a spheroid illuminated by a plane wave of $\lambda = 0.6328 \ \mu m$ in perpendicular polarization.</figcaption></figure>

For further study of influences of aspect ratio on intensity distribution, Fig. 4.17 illustrates first and second order rainbow near the rainbow angle when $b = 50, 55, 60$ and $65 \ \mu m$, i.e. aspect ratio $\kappa = b/a = 1.0, 1.1, 1.2$ and $1.3$. From the figure we find that, compared with rainbow of a sphere, when the transversal semi-axis $b$ increases, the intensity of the second order rainbow also increases while that of the first order rainbow decreases, leading that the relative intensity of first to second order rainbow decreases. When $\kappa = 1.3$, the second order rainbow is even in the same order as the first order rainbow.

From the above analysis, we see that the intensity ratio between the different orders of rainbow is sensible to the aspect ratio of the prolate spheroid. So we can use the intensity ratio to measure the deformation of a particle. To quantify this relation, we compile in Tab. 4.4 the intensity ratio of the second to fifth order rainbows ($p = 3$ to $5$) to that of the first order rainbow ($p = 2$) for the 4 aspect ratios ($\kappa=1, 1.1, 1.2, 1.3$). The two incident angles $\alpha_1$ and $\alpha_2$, the corresponding amplitudes and the curvature radii of the emergent rays, as well as the intensity of the main peak are calculated with the method described in the previous section. We find that the ratio of the second order rainbow to the first order rainbow increases as function of the aspect ratio, while the relation between the intensity ratio and the aspect ratio is not monotone for the higher order rainbows.

It is clear that relative intensity of different order rainbows depends on the deformation of the particle due to the curvature of the reflected and refracted waves. VCRM provides a powerful tool to this study.
<table><caption>Table 4.4: Intensity ratio of (p − 1)th order rainbow (main peak K = 0) with that of the first order rainbow (p = 2) of a prolate particle of a = 50 µm, m = 1.333 with aspect ratio κ = b/a = 1, 1.1, 1.2 and 1.3.</caption><thead><tr><th>p</th><th>κ</th><th>α (°)</th><th>R'<sub>1,q</sub> (μm)</th><th>R'<sub>2,q</sub> (μm)</th><th>√|D<sub>p</sub>|</th><th>ε<sub>p,1</sub></th><th>A<sub>p,1</sub></th><th>I<sub>p,1</sub></th><th>χ<sup>p</sup><sub>2,1</sub></th></tr></thead><tbody><tr><td rowspan="6">p=2</td><td>1</td><td>50.90<br>67.50</td><td>125.50<br>-28.70</td><td>76.19<br>72.28</td><td>74.56<br>58.03</td><td>0.24<br>0.35</td><td>18.24<br>20.30</td><td>1485.3</td><td>1.000</td></tr><tr><td>1.1</td><td>50.90<br>67.50</td><td>125.50<br>-28.70</td><td>58.50<br>51.34</td><td>70.74<br>53.46</td><td>0.24<br>0.35</td><td>17.30<br>18.70</td><td>1296.5</td><td>1.000</td></tr><tr><td>1.2</td><td>50.90<br>67.50</td><td>125.50<br>-28.70</td><td>46.63<br>38.49</td><td>69.93<br>52.03</td><td>0.24<br>0.35</td><td>17.10<br>18.20</td><td>1246.6</td><td>1.000</td></tr><tr><td>1.3</td><td>50.90<br>67.50</td><td>125.50<br>-28.70</td><td>37.35<br>28.86</td><td>70.63<br>52.05</td><td>0.24<br>0.35</td><td>17.27<br>18.21</td><td>1259.1</td><td>1.000</td></tr><tr><td>1</td><td>64.00<br>79.40</td><td>54.01<br>-2.78</td><td>-12.06<br>-27.98</td><td>42.52<br>27.16</td><td>0.12<br>0.25</td><td>5.30<br>6.68</td><td>143.4</td><td>0.097</td></tr><tr><td>1.1</td><td>64.00<br>79.40</td><td>54.01<br>-2.78</td><td>-68.97<br>-85.81</td><td>56.78<br>35.78</td><td>0.12<br>0.25</td><td>7.08<br>8.79</td><td>251.8</td><td>0.194</td></tr><tr><td rowspan="6">p=3</td><td>1.2</td><td>64.00<br>79.40</td><td>54.01<br>-2.78</td><td>-304.16<br>-328.34</td><td>102.27<br>63.81</td><td>0.12<br>0.25</td><td>12.75<br>15.68</td><td>808.1</td><td>0.648</td></tr><tr><td>1.3</td><td>64.00<br>79.40</td><td>54.01<br>-2.78</td><td>460.04<br>470.99</td><td>119.73<br>74.81</td><td>0.12<br>0.25</td><td>14.92<br>18.38</td><td>1109.3</td><td>0.881</td></tr><tr><td>1</td><td>69.20<br>84.20</td><td>35.46<br>0.39</td><td>75.75<br>66.57</td><td>36.31<br>19.20</td><td>0.07<br>0.18</td><td>2.70<br>3.55</td><td>39.1</td><td>0.026</td></tr><tr><td>1.1</td><td>69.20<br>84.20</td><td>35.46<br>0.39</td><td>30.73<br>20.16</td><td>30.19<br>15.84</td><td>0.07<br>0.18</td><td>2.24<br>2.93</td><td>26.8</td><td>0.021</td></tr><tr><td>1.2</td><td>69.20<br>84.20</td><td>35.46<br>0.39</td><td>5.14<br>-5.48</td><td>30.59<br>15.96</td><td>0.07<br>0.18</td><td>2.27<br>2.95</td><td>27.3</td><td>0.022</td></tr><tr><td>1.3</td><td>69.20<br>84.20</td><td>35.46<br>0.39</td><td>-21.16<br>-31.77</td><td>33.94<br>17.65</td><td>0.07<br>0.18</td><td>2.52<br>3.26</td><td>33.5</td><td>0.027</td></tr><tr><td rowspan="6">p=4</td><td>1</td><td>72.10<br>87.05</td><td>27.12<br>0.75</td><td>-28.43<br>-43.29</td><td>27.57<br>11.11</td><td>0.05<br>0.13</td><td>1.32<br>1.45</td><td>7.7</td><td>0.005</td></tr><tr><td>1.1</td><td>72.10<br>87.05</td><td>27.12<br>0.75</td><td>-726.16<br>-746.17</td><td>93.34<br>36.95</td><td>0.05<br>0.13</td><td>4.46<br>4.83</td><td>86.4</td><td>0.067</td></tr><tr><td>1.2</td><td>72.10<br>87.05</td><td>27.12<br>0.75</td><td>101.27<br>93.66</td><td>34.82<br>13.91</td><td>0.05<br>0.13</td><td>1.66<br>1.82</td><td>12.1</td><td>0.010</td></tr><tr><td>1.3</td><td>72.10<br>87.05</td><td>27.12<br>0.75</td><td>45.39<br>36.51</td><td>27.62<br>10.99</td><td>0.05<br>0.13</td><td>1.32<br>1.44</td><td>7.6</td><td>0.006</td></tr><tr><td>1</td><td>72.10<br>87.05</td><td>27.12<br>0.75</td><td>-28.43<br>-43.29</td><td>27.57<br>11.11</td><td>0.05<br>0.13</td><td>1.32<br>1.45</td><td>7.7</td><td>0.005</td></tr><tr><td>1.1</td><td>72.10<br>87.05</td><td>27.12<br>0.75</td><td>-726.16<br>-746.17</td><td>93.34<br>36.95</td><td>0.05<br>0.13</td><td>4.46<br>4.83</td><td>86.4</td><td>0.067</td></tr></tbody></table>
<figure><img src="image_36.png" /><figcaption>Figure 4.17: First and second order rainbow of an spheroidal water droplet with a = 50 µm and b = 50, 55, 60, 65 µm (i.e. κ=1.0, 1.1, 1.2, and 1.3) illuminated by a plane wave of λ = 0.6328 µm in perpendicular polarization.</figcaption></figure>

### 4.2.5 Hyperbolic umbilic foci of an oblate particle

When an oblate droplet is illuminated by a plane wave in the direction perpendicular to the symmetrical axis, the slight eccentricity causes a bending of the rainbow in this plane. Decreasing the aspect ratio $κ = b/a$ leads to caustics arising from the interference of two light rays in the equatorial plane and two skew light rays which do not lie in the horizontal plane en entering the drop but confined to the horizontal direction upon exiting the droplet. As the ratio is further decreased, the four light rays come out with the same angle, i.e. at the rainbow angle. The resulting pattern is that of a hyperbolic umbilic (HU) focal section. Nye [83] gave the critical aspect ratio for such hyperbolic umbilic section as:

$$κ = [4(m^2 - 1)/3m^2]^{1/2} \qquad (4.62)$$

It is approximately 0.764 for the refractive index $m = 1.333$.

These phenomena can also be explained analytically in the framework of VCRM. As an example, we reproduce in Fig. 4.18 the figure in the paper of Onofri et al [1]. The properties of the particles are given in the figure and the wavelength of the incident wave is 532.130 µm. Apart from the first rainbow in the angle around 153° we observe another peak hyperbolic-umbilic diffraction catastrophe (HUDC) whose position depends on the aspect ratio and refractive index. In the framework of
VCRM, both the rainbow position and the second peak can be explained by the divergence factor of a spheroidal particle given in Eq. (4.41).

<figure><img src="image_37.png" /><figcaption>Figure 4.18: Comparison of VCRM and experimental normalized equatorial scattering diagrams of a Di-Ethyl-Hexyl-Sebacat droplet [1]</figcaption></figure>

The first denominator

$$T_1 = m \cos \beta - 2 \cos \alpha$$

is independent of the aspect ratio. When it tends to zero, the divergence factor, therefore the intensity of the scattered light, tends to infinity, that is the case of the first rainbow.

The second denominator

$$T_2 = 2ma^2 \cos^2 \beta - 2a^2 \cos \alpha \cos \beta - mb^2$$

depends on the aspect ratio. When it tends to zero, the divergence tends to infinity. We can evaluate the position of the peak by this term. The solution of $T_2 = 0$ for the three cases in the figure is compiled in Tab 4.5.

The third term

$$T_3 = (2mb^2 \cos \beta - b^2 \cos \alpha - 2ma^2 \cos^3 \beta + 2a^2 \cos \alpha \cos^2 \beta)$$

is positive in the three case presented in the figure. Therefore it does not generate a complementary peak.
<table><caption>Table 4.5: Variation of the HUDC peak as for an oblate particle of different refractive index and aspect ratios.</caption><thead><tr><td>case</td><td>refractive index</td><td>aspect ratio</td><td>α(°)</td><td>θ(°)</td></tr></thead><tbody><tr><td>1</td><td>1.4505</td><td>0.8938</td><td>69.69</td><td>158.25</td></tr><tr><td>2</td><td>1.4465</td><td>0.9608</td><td>82.05</td><td>171.26</td></tr><tr><td>3</td><td>1.4446</td><td>0.9816</td><td>85.14</td><td>175.84</td></tr></tbody></table>

## 4.3 Conclusion

In this chapter, we applied VCRM first to revisit the Airy theory and then to investigate the properties of the rainbows formed by a prolate and an oblate spheroidal particles. We have shown that the VCRM predict better the intensity and positions of the supernumerary bows of rainbow of a spherical particle than the Airy theory. The precision in the particle sizing by the supernumerary bows of VCRM is ten times better than Airy theory for the first order rainbow. However the VCRM is to be improved for the prediction of the main peak position.

VCRM permits also to calculate analytically the positions and the intensity of different order of rainbow of a spheroidal particle when it is illuminated in the symmetrical plane. The peak due to the focal of an oblate spheroidal particle can also be evaluated analytical by using the divergence factor of VCRM.

It is shown that the intensity ratio of different orders of rainbows is very sensible to the aspect ratio of a prolate particle. This property can be applied to characterize the deformation of particles.
# Chapter 5

## Experimental study of scattering
by a pendant drop

A rainbow measurement system by light scattering of pendant drops is set up in the laboratory. In this chapter, the experimental investigation of light scattering by such particles in the vicinity of rainbow angle is shown. By comparisons of experimental results with the predictions of VCRM, the validity of VCRM for prolate spheroidal droplets is confirmed. The influence of particle ellipticity on the intensity distribution of the scattered light is discussed.

Because the current numerical implementation of the VCRM addresses only the scattering of a plane wave in the equatorial plane of an ellipsoidal particle, our study is focused on the symmetric plane of pendant drops, more specifically, in the angular region near the rainbow angle.

### 5.1 Experimental setup

#### 5.1.1 Consideration of the measurement system

Simple but effective, the rainbow patterns of pendant water drops are full of information and suitable to validate VCRM and study the influences of particle ellipticity on the intensity distribution of scattered light. To do this, the experiment system requires high quality rainbows for pendant droplets of different shapes. Therefore following aspects are critical and considered in the design of the measurement system:

* Rainbows of high quality: a stable rainbow of high quality is very important in our experiment. To achieve it, a high quality laser, a well-designed optical path,
<figure><img src="image_38.png" /><figcaption>Figure 5.1: Experimental setup</figcaption></figure>

as well as stable and homogeneous water droplets are indispensable. Meanwhile,
the noise from environment should be reduced.

* Flexibility of system: to study the influence of particle's non-sphericity on the rainbow pattern, like the relation between the intensity ratio of the first to the second order rainbows and the deformation degree of drops, the system should have the ability to adjust the shape of droplets flexibly. Besides, the acquisition system of scattering pattern needs to be synchronized with the imaging camera of the droplet.

5.1.2 Composition of the measurement system

To meet the above demand, the following rainbow measurement system is set up in
the laboratory. The photo of the system and the schema of the configuration are
shown in Figs. 5.1 and 5.2. It is composed of 4 parts:

1. generator of droplets;

2. optical system to form and adjust a laser beam illuminating the droplets;

3. acquisition system to detect the scattered light and record the image of the droplet

4. motorized rotation platform.

The generator of droplets is made up of a water tank under pressure of 1 to 10 bar,
a rotameter (Model 1355, Brooks Instrument) and circular capillary tubes of inner
<figure><img src="image_39.png" /><figcaption>Figure 5.2: Schema of experimental set-up</figcaption></figure>

diameters vary from 0.6 to 0.95 mm. The water is stored in the tank of 6 L and flows through a pipe to the rotameter, which is connected to the capillary tube.

A beam of wavelength $\lambda = 0.6328$ µm is emitted by a high quality He-Ne laser (JDSU He-Ne Laser, Model 1145, Mode Quality > 95% and non-polarized). Its beam diameter at exit $w_{exit} = 0.7$ mm. A polarizer and a half-wave plate are used for shifting the polarization direction of the beam emitted by the laser. Then the beam is reflected by two mirrors which are designed to facilitate the adjustment of the beam. The laser beam is expanded by an expander of expansion ratio 20, so that the beam is large enough compared to the droplet. After the expander, a diaphragm is used to adjust the position of the droplet illuminated by the beam. A 4F optical system is set up after the diaphragm. This system consists of 2 two bi-convex lenses of same focal length and can be used for Fourier transform image processing. Finally the laser interacts with the suspended droplet or liquid jet.

The acquisition system consists of two camera, one to register the scattering pattern and the other to record the image of the droplet. The light scattered by the droplet is acquired by camera 1, which is a 1000 x1000 pixels, 14-bit CCD camera (C9100-02, Hamamatsu). It is mounted on a rail fixed to the rotation platform. A white screen is put before camera 1 in order to assure a sufficient large field of view. The distance $d$ between the screen and the droplet is 60 mm. The acquisition camera can also be replaced by our linear CCD (Spyder 3 SG-14, DALSA) in some cases, for example, acquisition of light scattered by a liquid jet, for which, fast acquisition is needed. And our linear CCD, because of its high responsivity (2064 DN (nJ/cm²))
at 0 dB gain), high scan frequency (68 kHz) and high dynamic range (1400 : 1), is
very suitable for fast acquisition. The image of the droplet is registered by camera 2,
a 2048x2048 pixels, 10-bit CCD camera (Jai TM-4200CL). And the imaging system
was equipped with a long distance microscope allowing ~ 2.5 X magnification of the
droplet images. To calibrate the magnification factor of the lens system, one records
an image of a grid positioned at the exit of capillary tube. Subsequently, the relation
between the pixel number and the droplet size is found. The scale of the resulting
spatial intensity data was 460 pixels/mm, and the resolution of the complete imaging
system was on the order of 5 µm. The two cameras are synchronized by BNC Pulse
Generator.

The rotation platform consists of two motorized rotation stages (RSA200, Zolix). They are driven by a stepper motor controller (SC300, Zolix). The controller can be connected to a PC for multi-axis motion control and parameters setting with the software provided. It can also be programed by the communication protocol. In addition, the default parameter settings saved within the controller allow remote operations without computer. The two rotation stages rotate on same axis and can provide continuous rotation of arc-second resolution (0.00125°). A 5-Axis stage with the capillary tube is mounted on the upper rotation stage. The 5-Axis stage can provide five adjustments, three translation and two rotation. These five independent degrees of freedom provide 13 mm of horizontal translation, 6 mm of vertical translation, and 10° of pitch and yaw adjustment. By the adjustment of 5-Axis stage, we make sure that the capillary tube is located on the axis of rotation stage and suspended drops / liquid jets are illuminated by laser in the direction perpendicular to their axis. In this way, the upper rotation stage can be used to adjust the incident direction of beam by rotating the capillary tube. The lower stage is connected with the rail where camera 1 is mounted. Thus we can measure the scattered light in all directions.

5.1.3 Acquisition and image processing software

**Acquisition software for camera** Two different softwares HoKaWo and Dynam-
icStudio are used for the image acquisition of camera 1 and 2 respectively. They
include all basic functions required to acquire and process images. Both acquisition
softwares support the external trigger mode, making it possible to synchronize the
two cameras by BNC Pulse Generator. For the linear CCD, based on the develop-
ment platform of DALSA, an image acquisition and camera control software has been
developed. This homemade software has a user-friendly interface (Fig. 5.3) and all
the basic functions like camera configuration and image acquisition. The intensity
diagram can also be exported for further treatment.
<figure><img src="image_40.png" /><figcaption>Figure 5.3: Interface of the image acquisition and camera control software for our linear CCD</figcaption></figure>

**Image processing software** The image of scattering pattern is processed by ImageJ, a well-known image processing software. It was designed with an open architecture that provides extensibility via Java plugins and recordable macros. Custom acquisition, image analysis and processing plugins can be developed using ImageJ's built-in editor and a Java compiler. In our experiment, three macros of image processing are developed to meet our demand. Their purposes and principles will be shown in section 5.3. The details of the macros will be presented in Appendix B.

## 5.2 Experimental observations

Fig. 5.4 illustrates four images of scattering diagrams around the first and the second rainbows. Image 5.4(a) is a typical image of scattering diagram of a drop almost spherical. The rainbows are similar to what we observe in the nature. Two visible bows correspond to the primary (right) and the secondary (left) order rainbows. The supernumerary bows, well known as Airy structure, are clearly visible in both the first and the second rainbows. Each rainbow consist of a main peak and the subsidiary bows. For a sphere, Airy structure has a concave shape bending over toward the backscattering, leaving a constant angular width, named as Alexanders dark band. When pendant drop becomes bigger and is elongated in the vertical direction, correspondingly, the bows are deformed. The rainbow patterns in images 5.4(b) and 5.4(c), notably the second order on the left, are more curved. In image 5.4(d), the second order rainbow is completely twisted. The forming mechanism of these twisted structure will be discussed lately. Besides by simple observation we find that the intensities of
the two orders of rainbow in image 5.4(a) and image 5.4(b) are about in the same
order. However when the drop is deformed, the intensity of the first rainbow becomes
very weak compared to the second one in image 5.4(c) and 5.4(d). This phenomena
agrees with the observation of Sassen[92] (section 1.3.3) and attracts our attention. In
this thesis, we are mainly concerned about the relation between the rainbow intensity
and the curvature radii of droplets in the equatorial plane.

5.3 Image processing

In previous experiment, scattering pattern near the rainbow angle, as well as images
of drops, are recorded by two synchronized cameras. For the quantitative comparisons
with VCRM, the image processing is needed.

Firstly, the noises from environment, especially the light reflected by the tube
exist in the scattering pattern recorded in camera 1. They can not be neglected.
Therefore, before the comparison, we should subtract the noise from the images of
scattering pattern.

Secondly, we need to extract two principal curvature radii of the droplet in the equatorial plane for simulation. In VCRM, the two principal curvature radii of the droplet in the equatorial plane are needed for the calculation of the scattering diagram in this plane.

Finally, it is needed to get the scattering diagram in equatorial plane from images
of scattering pattern. The existing software, ImageJ for example, can get intensity
distribution of an indicated region from images, but when the shape of drop changes,
its equatorial plane moves vertically. We need to locate this plane in the image.
Besides, for all the image processing procedures above, batch processing function is
also essential.

To meet these demands, three plugins of ImageJ are developed in Macro language which is a scripting language built into ImageJ and allows to automate a series of ImageJ commands. Programs written in macros, Macros for short in this thesis, can be used to perform sequences of actions.

The first Macro is for batch background noise reduction (see Appendix B.1). The
principle of this macro is that we define the image without droplet as a background,
and then subtract it from the particle images. The existing ImageJ command can
realize the calculation between two images, but it is limit to the number of images and
cannot meet our demand of batch processing. So a macro of images batch processing
is developed here, to automate the subtraction of background noise from the images
of droplets.

The second Macro is developed to extract the two principal curvature radii of the
<figure><img src="image_41.png" /><figcaption>Figure 5.4: Images of pendant drops and their scattering pattern near rainbow angle. The region between two yellow line is the equatorial region selected. The yellow and red curve along the rainbow pattern are respectively the outline extracted and curve fitted for the main peak of rainbow.</figcaption></figure>
drop in the equatorial plane (see Appendix B.2). This Macro will make binary the image to separate the droplet from the background; extract the contour of droplet near equatorial plane by Moving Average Method; curve fit the contour; locate the equatorial plane and calculate two principal curvature radii in this plane from the approximating function.

The third Macro is for the calculation of the intensity ratio of the two orders rainbows (see Appendix B.3). When the equatorial plane is found, we can export the intensity distribution at this plane, and calculate the intensity ratio of the first to the second order rainbow. The following cases should be considered:

**location of the equatorial plane** When the drop is nearly spherical like in image B.7(a), both the first and the second rainbow have a concave shape and symmetry to the equatorial plane. Therefore once their main peaks are extracted, we can fit each of them with a quadratic function respectively and locate the equatorial plane by calculating vertexes of these quadratic functions. However, as observed before, when the pendant drop becomes bigger, the scattering patterns changes. The second order rainbows in images B.7(b), B.7(c) and B.7(d) are deformed and no longer symmetrical. On the contrary, the first order rainbow remains symmetrical relative to the equatorial plane. In this case, we determine the equatorial by the first order rainbow.

**Selection of a region near equator for measure of scattering diagram** Once the equatorial plane in an image of scattering diagram is located, the intensity distribution in this plane can be got. Image 5.5(a) is an image of 2D scattering diagram, the yellow line is the calculated equatorial plane. The intensity distribution in the equatorial plane is shown in image 5.5(b). From image 5.5(b), we can find that noise cannot be ignored. The intensity ratio of the first to the second order rainbow cannot be calculated. Therefore a region of wide w=40 pixels are selected near the equatorial plane to calculate the mean value of the intensity in image 5.6(a). After such treatment the intensity distribution of image 5.6(b) becomes rational, both first and second order rainbows are identifiable for further calculations.

**Conversion of intensity** In experiment, the main peak value $G_{value}$ exported from scattering diagram is the gray level. It is proportional to its illuminance. Furthermore, all the elements of the screen are not perpendicular to the light rays arriving on it and the difference of their distance from the drop is not negligible. Thus for comparison with VCRM, we define a dimensionless factor E referring Fig. 5.2 by

$$E = \frac{G_{value}}{\cos \alpha} \cdot r^2 = \frac{G_{value}}{h} \cdot r^3 \qquad (5.1)$$

This factor is theoretically proportional to the distance free intensity defined in Eqs. (2.24) and (2.25).
<figure><img src="image_42.png" /><figcaption>Figure 5.5: image of scattering diagram, and its intensity distribution in equatorial plane got from experiment</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>Figure 5.6: image of scattering diagram, and its intensity distribution in a selected window (wide =40 pixels) near equatorial plane</figcaption></figure>
<table><caption>Table 5.1: Comparison of the rainbow intensities of experimental measurement with the numerical results simulated by VCRM.</caption><thead><tr><th>Image Number</th><th>Minor and major axes (μm)</th><th>Experiment E<sub>max1</sub>/E<sub>max2</sub></th><th>Theory I<sub>max1</sub>/I<sub>max2</sub></th><th>Theory/Experiment</th></tr></thead><tbody><tr><td>1</td><td>(807, 849)</td><td>1.65</td><td>2.16</td><td>1.31</td></tr><tr><td>2</td><td>(1067, 1149)</td><td>1.63</td><td>1.98</td><td>1.21</td></tr><tr><td>3</td><td>(1195, 1325)</td><td>1.44</td><td>1.79</td><td>1.24</td></tr><tr><td>4</td><td>(1292, 1418)</td><td>1.39</td><td>1.89</td><td>1.36</td></tr><tr><td>5</td><td>(1362, 1548)</td><td>1.29</td><td>1.68</td><td>1.30</td></tr><tr><td>6</td><td>(1393, 1651)</td><td>0.76</td><td>1.17</td><td>1.53</td></tr><tr><td>7</td><td>(1431, 1765)</td><td>0.74</td><td>0.95</td><td>1.28</td></tr><tr><td>8</td><td>(1453, 1864)</td><td>0.66</td><td>0.91</td><td>1.38</td></tr><tr><td>9</td><td>(1476, 1798)</td><td>0.62</td><td>0.82</td><td>1.32</td></tr></tbody></table>

## 5.4 Results and conclusions

In this section, we will compare experimental results with our predictions of VCRM. Because we are interested only in the scattering intensity distribution in the equatorial plane, only the geometry of the droplet in this plane plays its role. We need only the local curvature radii in the equatorial (horizontal) plane $a_1$ and in the vertical plane $a_2$. For better understanding, the pendant drop can be considered as a prolate spheroid of semi-minor axis equal to $a_1$ and the semi-major axis $b$ equal to $b = \sqrt{a_1a_2}$. The minor and major axes extracted from a series of 9 images are given in Tab 5.1. Based on them, our code of VCRM developed in C# can calculate the full scattering diagram in the equatorial plane of pendant drops. The intensity ratio of the first to the second order rainbows are then extracted from the scattering diagram calculated. The experimental intensity ratios of first to second order rainbows are extracted from the corresponding images of scattering diagram recorded by camera 1. From Tab 5.1 we can find firstly that as indicated above, in contrary to the case of spherical particle, the intensity of the first order rainbow may be weaker than that of the second order rainbow for pendant drops. Besides, it is found that, even the theoretical ratios are around 30% larger than the ratios obtained from the experiment, their trends are in good agreement.

Fig. 5.7 shows the intensity ratios of the first to the second order rainbow obtained as function of aspect ratio $b/a_1$. The intensity ratio of an equivalent sphere is also added for comparison. We know that, for a sphere, the intensity of the first rainbow is always stronger than the second one because the loss of energy in the reflection and in the divergence of the wave on the particle surface. The theoretical simulation in this figure shows that the intensity ratio of the first to the second order rainbow of equivalent spheres in the range of size 800 to 1500 μm, is almost constant and equal to about 2.5. While the intensity ratio of the pendant drops decreases as the aspect ratio. The simulations by VCRM for prolate spheroids are in agreement with the experimental results. The decrease of the intensity ratio as function of the aspect
<figure><img src="image_44.png" /><figcaption>Figure 5.7: Intensity ratio of first to second order rainbow as function of aspect ratio.</figcaption></figure>

ratio can be explained by the vertical convergence. In fact, the rays of the first rainbow experience one internal reflection while the second ones experience two internal reflections. Therefore, the focalization effect decreases faster in the second rainbow than the first rainbow as the aspect ratio increases.

From the above studies, we find that the scattering properties of non-spherical particles can differ dramatically from those of spheres. The intensity ratio of the first to the second order rainbows is very sensitive to the deformation. It is therefore can be used as a parameter to measure the deformation of the deformation of the droplets. Meanwhile, in this application of VCRM in the simulation and interpretation of scattering patterns of naturally generated pendant drops, we have shown that this model can predict correctly complex structure in the scattering diagrams of non-spherical particles. The results demonstrate the large potential of the model in the characterization of non-spherical particles. Besides, we observe in the experiment many interesting phenomena near the rainbow angle and in forward direction (not shown in this thesis). However by the limitation of the numerical simulation with 3D VCRM, our study is limited to the equatorial plane of the drops. For further studies of the influence of particle non-sphericity on the rainbow pattern, 3D code of VCRM is essential.
98CHAPTER 5. EXPERIMENTAL STUDY OF SCATTERING BY A PENDANT DROP
# Chapter 6

# Conclusions and Perspectives

This thesis is devoted to the extension of Vectorial Complex Ray Model (VCRM) and its application to the metrology of non-spherical particles. In this chapter we will give the conclusions of presented work and some perspectives in further studies.

## 6.1 Conclusions

In various research domains, such as the environment control, the biochemistry, the fluid mechanics and the combustion field, we need to measure the characteristics of the particles, like the particle size, shape, velocity, temperature etc. Among different kinds of measurement techniques, optical metrology is largely employed because of being accurate, reliable and non-intrusive. A large effort has been devoted to develop powerful optical techniques. However, in most optical techniques, the shape of the particles is considered to be spherical, while the scattering properties of non-spherical particles can differ dramatically from those of spheres. This is mainly due to a lack of theoretical model to describe the interaction of light with large non-spherical particles. Therefore characterization of non-spherical particles is still a veritable challenge in the development of optical metrology in very large scope of applications.

Various theoretical and numerical methods have been developed for non-spherical particle, such as spheroid, ellipsoid, or elliptical cylinder to take into account the shape of the particles. However, the sizes or shapes of the particles are often limited because of the numerical difficulty. The Vectorial Complex Ray Model (VCRM) recently developed in the laboratory permits to deal with the scattering of large irregularly shaped particles. It is based on the geometrical optics but a new property – the wave front curvatures is introduced. This makes VCRM possible to calculate the divergence factor of a single ray bundle and is easy to extend to irregularly shaped 3D objects. VCRM has been proved to achieve good agreement with the rigorous theories. In this thesis, I devoted myself to the further extension of VCRM and its application to large
non-spherical particles.

After a brief recall of the fundamental models and concepts in light scattering, the general principles of the VCRM are presented. The model is then applied to the scattering of light in the symmetric plane of an ellipsoid. The scattering diagrams of spherical particles by VCRM are compared with those of Lorenz-Mie theory (LMT) – a rigorous theory. Good agreement can be observed between VCRM and LMT for spherical particles of diameter larger than tens of wavelength. The effect of particle size on the precision of VCRM is examined. In addition, the scattering diagrams of spheroidal particles are presented to show the influence of particle’s non-sphericity on the rainbow pattern. For the scattering diagram of a prolate spheroid, compared with that of a sphere, its intensity of the second order rainbow increases and becomes very close to that of first order rainbow. For the scattering of light by an oblate spheroid, hyperbolic umbilic diffraction catastrophe occurs.

VCRM is also applied to revisit the Airy theory and then to investigate the prop-
erties of the rainbows formed by a prolate and an oblate spheroidal particles. The
comparisons show that VCRM predict better the intensity and positions of different
supernumerary bows of rainbow of a spherical particle than the Airy theory. The
precision in the particle sizing by the supernumerary bows of VCRM is better than
Airy theory for the first order rainbow. However the VCRM is to be improved for the
prediction of the main peak position.

The positions and the intensities of different orders of rainbow of spheroidal parti-
cles illuminated by a plane wave in the symmetrical plane are calculated analytically
in the framework of VCRM. The relation between the intensity ratio of different orders
of rainbows and the curvature radius of a prolate spheroid in the direction perpendic-
ular to the equatorial plane is quantified. The results show that the intensity ratio of
different order rainbow is very sensible to the aspect ratio of a prolate particle. This
property can be applied to characterize the deformation of particles. The peak due
to the focal effect of an oblate spheroidal particle is evaluated analytical by using the
divergence factor of VCRM.

Besides, a rainbow measurement system for light scattering by pendant drops is set
up in the laboratory. The experimental investigation of light scattering by a pendant
drop in the vicinity of rainbow angle is shown. By comparisons of experimental results
with predictions of VCRM, the validity of VCRM for prolate spheroidal droplets is
confirmed. The influence of particle's ellipticity on intensity distribution of scattered
light is discussed. We conclude that the focalization effect decreases faster in the
second order rainbow than the first rainbow as the vertical curvature radius increases,
leading that the relative intensity of first to second order rainbow decreases.