Inherent optical properties of the coccolithophore: Emiliania huxleyi

: A realistic nonspherical model for Emiliania huxleyi (EHUX) is built, based on electron micrographs of coccolithophore cells. The Inherent Optical Properties (IOP) of the EHUX are then calculated numerically by using the discrete dipole approximation. The coccolithophore model includes a near-spherical core with the refractive index of 1 . 04 + m i j , and a carbonate shell formed by smaller coccoliths with refractive index of 1 . 2 + m i j , where m i = 0 or 0 . 01 and j 2 = − 1. The reported IOP are the Mueller scattering matrix, backscattering probability, and depolarization ratio. Our calculation shows that the Mueller matrices of coccolithophores show different angular dependence from those of coccoliths.


Introduction
Coccolithophores, or coccolithophorids, are found in large numbers in global ocean waters.They are unicellular and photosynthetic organism, and provide important food sources for aquatic environments.An interesting feature of coccolithophores is that they generate calcium carbonate plates called coccoliths.The production of coccoliths plays important roles in the carbon cycle through their CaCO 3 production which responds to CO 2 partial pressure change [1,2].Emiliania huxleyi (EHUX) is the most abundant species of coccolithophores.They can produce massive blooms that have large impacts on the environment and fisheries.On the global scale, it is indispensable to use ocean color sensors such as MODIS and VIIRS to study the effects of coccolithophore primary production on the carbon cycle and to monitor EHUX bloom events [3].The CALIOP satellite lidar system provides valuable information on ocean waters [4,5].The combination of passive and active systems will generate a plethora of multidimensional information on coccolithophores and other oceanic particles.
To assist the interpretation of both passive ocean color satellite images and active lidar signals, it is necessary to have a knowledge of the Inherent Optical Properties (IOP) of these particles.A short list of these properties includes absorption and scattering coefficients, beam attenuation coefficient, scattering matrix, and a few other parameters used in different applications [6].Both theoretical simulations and experimental measurements have been used to obtain the IOP for oceanic particles [7].Experimental measurements of ocean waters are mostly focused on absorption and extinction coefficients [8], and sometimes the volume scattering function (the 11 element of the scattering matrix) [9][10][11].The NASA Ocean Biology Processing Group maintains a large collection of such measurements along with other bio-optical data [12].It is worth noting that the Radiance in a Dynamic Ocean (RaDyO) program sponsored by the Office of Naval Research has involved simultaneous measurements of many IOP and Apparent Optical Properties (AOP) in order to understand underwater visibility and imaging of objects across the air-sea interface [13].Complete Mueller matrix measurements for ocean water and a few phytoplankton species have been reported [14,15].Shapiro et al. has measured the S 14 Mueller matrix elements of single, immobilized dinoflagellates as well as on suspensions of the species Prorocentrum micans to study the scattering of circularly polarized light [16].Recently Svensen et al. reported measurement of Mueller matrices for a variety of phytoplankton cultures [17].Volten et al. have reported the S 11 and S 12 Mueller matrix elements of phytoplankton and silt [18].
Microscopic images of coccolithophores and other particles show that overall these particles are not spherical.Nevertheless spherical models are used to simulate the IOP for oceanic particles due to the significant difficulties of simulating light scattering by large nonspherical particles [19].Limited cases have been reported on simulating the IOP of nonspherical marine particles, mainly by Gordon and coauthors [20][21][22][23][24][25].Notably, the simulated backscattering cross sections of detached coccoliths have been compared with the estimation derived from insitu radiance/irradiance profiles and detached coccolith concentration [24].The agreement is quite reasonable.In this paper we report a simulation of the IOP for whole EHUX cells, which to the best of our knowledge have not previously appeared in the literature.A realistic model has been built based on Scanning Electron Micrograph (SEM) images of the coccolith and coccolithophore (see photos by Jeremy Young, Figs. 1 and 2 in [21]).The simulation of light scattering by these nonspherical particles is done using the Amsterdam Discrete Dipole Ap-proximation code (ADDA) [26].The Mueller matrices of both coccoliths and coccolithophores are simulated and impacts of size and morphology are studied.

Coccolith model in a local coordinate system
Our model is based on the SEM of single coccolithophore cells (see Figs. 1 and 2 in [21]).The coccolith model is similar to the wagon wheel model presented in [22], but with curvature taken into account.Figure 1(a) shows a three dimensional view of the coccolith model.To define the shape of the coccolith, we use a local Cartesian coordinate system and any point in the space is denoted as r = (x, y, z).The coccolith model has a symmetry axis (the broken blue line in Fig. 1(a)), which passes through the origin O of this local coordinate system.The symmetry axis is coincident with the z-axis of the local coordinate system.A coccolith with its axis pointing in an arbitrary direction will be obtained by rotating the original coccolith using Euler's rotation theorem in the following section.In this paper we also use a polar coordinate system within which a point is denoted as r = (r, θ , φ ), where r = x 2 + y 2 + z 2 , x = r sin θ cos φ , y = r sin θ sin φ , z = r cos θ .The length of the projection of r on the x − y plane is ρ xy = x 2 + y 2 .The whole Cartesian space is denoted as V so that any point r belongs to V , i.e., r ∈ V .The space region occupied by the coccolith model shown in Fig. 1(a) is denoted by V c , which consists of three parts (regions): an inner spherical cap V I, an outer spherical cap V II, and a hollow cylinder connecting the two caps V III.The spherical caps are parts of concentric spherical shells, with the center of the spherical shells defined as the coordinate origin O. Figure 1(b) shows a sketch of a cross section of the coccolith in a half plane containing the origin and axis.The radius and thickness of the inner cap V I are r 0 and d c , respectively.The hollow cylinder has an inner radius of r 1 .The thickness of the cylinder is also d c .The height of the cylinder V III is d h , which is defined as the distance between the two centers where two caps intersect with the cylinder in Fig. 1(b).The thickness of the outer spherical cap is assumed to be the same as the inner cap.The projection of the inner or outer spherical cap in the x-y plane is shown in Fig. 1(c).It is fairly obvious that this projection is very similar to the pin-wheel model by Gordon (see Fig. 3 in [21]).However, it should be noted that our model consists of spherical caps with which curvature is built in.In addition, we use two rings, one inner and one outer, around the spherical pinwheel to make the model closer to the SEM of the detached coccolith.
Assuming the number of the pinwheel is N = 40, and N = 2 * N , we have the following definitions: The space V I that the inner spherical cap occupies can be defined in the following way: where ∩ and ∪ are used to show the intersection and union of two space volumes, respectively.In a similar way the outer spherical cap is defined as: The hollow cylinder region, V III, is defined as: The detached coccolith is then: This coccolith model is defined in the local coordinate system.In the next section, we will obtain the coccolith model in the global coordinate system within which its symmetry axis may point in an arbitrary direction.

Coccolith model in a global coordinate system
The global coordinate system shares the same origin with the local coordinate system.For any point in the space, its coordinates can be expressed in both the local and global coordinate systems.We denote the local and global coordinates by r l = (x l , y l , z l ) T and r g = (x g , y g , z g ) T , respectively, where the superscript T stands for the transpose of the column vector.They can be related with a matrix according to Euler's rotation theorem: If r l is a point within V c given in Eq. ( 5), r g will be the same point expressed in the global coordinate system.To find the expression of M, we assume that rc = (sin θ c cos φ c , sin θ c cos φ c , cos θ c ) is the unit vector in the global coordinate system pointing in the direction of the coccolith axis.The local x-axis could be defined as x l = rc × ẑg /|r c × ẑg | in the global coordinate system, where ẑg = (0, 0, 1) T is the unit vector of the z-axis of the global coordinate system; and × denotes the vector cross product.The local y axis is then ŷ l = rc × x l , because rc is the local z-axis expressed in the global coordinate system.The transformation matrix is written as: because this form of M can transform (1, 0, 0) T , (0, 1, 0) T , and (0, 0, 1) T into x l , ŷ l , and ẑ l , respectively.The transformation from the global to local coordinate system is then: where M −1 is the inverse matrix of M.
There are two steps to determine if an arbitrary point r g in the global coordinate system is occupied by a coccolith: 1. transform r g to r l using Eq. ( 8) with M specified by θ c and φ c ; 2. use Eq. ( 5) to find whether the point is in one of the coccolith domains.

Coccolithophore model
Two models, symmetrical and asymmetrical coccolithophores, are built to study the sensitivities of IOP on the microscopic morphologies.Figures 2(a The refractive index m of the spherical core of the symmetrical coccolithophore is given by: where m i is 0 and 0.01 in our study.These values are relative to pure water and consistent with [18].The value m i = 0.01 is used to study the impacts of absorption on the IOP of coccoliths and coccolithophores.
To form a carbonate shell, we use a set of coccoliths with different values of θ c and φ c .They are evenly distributed on and cover the surface of the spherical core with the radius of r 0 .The number of θ c values is: as we intend to cover the zenith angle range from 0 to π with coccoliths whose half angle span is θ 3 .The NINT () operator takes the integer part of the enclosed real number.The θ c values are: For each value of θ c , the number of φ c values are given by: The φ c is: For each pair of θ c and φ c values, the coccolith domain is defined by the two step procedure described in Sec 2.2.If a point is in the coccolith domains, the refractive index is specified as: which is defined relative to water [21].There are cases where a point in the space could belong to multiple coccoliths.This does not matter because the same refractive index value is used for all coccoliths.We just assign the refractive index of 1.2 to the point if this happens.Figure 2(a) shows regions that different coccoliths entangle with each other, manifesting this fact.By a scalar value for the refractive index, we have ignored the birefringence of the carbonate shell.Using the Rayleigh-Gans scattering approximation, Gordon has demonstrated that the relative difference between a randomly oriented birefringence disk and an isotropic disk is within 20% for backscattering cross section [23].Based on this argument, birefringence has been neglected for coccolith simulation in Gordon et al. [24].Implementing birefringence for coccolithophores is quite complicated and needs significant additional efforts, which we would defer to a future study.The asymmetrical coccolithophore is designed to mimic a coccolithophore more realistically, as observed in many SEM images.Starting from the symmetrical model, we put another layer of coccolith shell surrounding the symmetrical model shown in Fig. 2(a).To do this, r 0 = r 0 + d h + 2d c is used to replace r 0 from Eq. ( 1) to (5).In addition, the condition of η > 0.5 (15) is used to determine whether a coccolith in the extra outer layer should be kept before we put it in the model, where η is a random number uniformly distributed in the interval of [0, 1].This will effectively reduce the number of coccoliths by half in the additional layer of the carbonate shell.After the extra layer of coccoliths is added, we further compress the z coordinate values of all the points in the global space by 80%: z = 0.8z.The coefficient of 0.8 is determined by the ratio of the smallest and largest dimensions of the SEM shown in Fig. 1 of [21].Figure 2(b) shows one realization of the asymmetrical coccolithophore model.
In this work we fixed the ratios between different parameters: r 1 /r 0 = 0.18, r 2 /r 0 = 0.46, d c /r 0 = 0.07, and d h /r 0 = 0.18.The size of the particle is then primarily determined by r 0 .We have used r 0 = 1.5, 1.8, 2.1, 2.4, and 2.7μm in the simulation of light scattering by coccolithophores to study the effect of different sizes on the IOP.The maximum diameters of coccolithophores corresponding to these r 0 values are roughly: D = 2(r 0 + d h + d c ) = 3.75, 4.50, 5.25, 6.00 and 6.75μm.The values of r 1 , r 2 , d c , and d h could be set to vary independently if future research finds better statistical representations for them.For the case of coccoliths, we have studied the cases of r 0 = 3.0, 3.3, 3.6, 3.9, 4.2, and 4.5μm, in additional to those used for coccolithophore cases.The diameters of the corresponding coccolith projections (see Fig. 1(c)) range from roughly 1.63μm (for r 0 = 1.5μm) to 4.73μm for (r 0 = 4.5μm).Note that the values of r 0 > 3.0μm are probably higher than a typical value suggested by SEM images, which are used here for theoretical considerations.For light scattering by coccolithophores, r 0 > 3.0μm requires too much computational resource which are ruled out in this study.
In order to compare IOP for different particles, we use the volume equivalent size parameters x = 2πr v /λ to denote the particle size, where λ is the wavelength; r v = [3V e /(4π)] 1/3 is the volume equivalent radius; and V e is the total geometric volume occupied by the particles.The usage of size parameters is motivated and justified by the scale invariance rule, i.e., the scattering properties of a particle depend on its relative size to the wavelength, not on its absolute size (see Sec. 5.8.2 of [27]), if the morphology is the same.It is evident that for particles with structures as diverse as detached coccoliths, symmetrical and asymmetrical coccolithophores, the volume equivalent size parameter is only one of many parameters which influence the particles interaction with the electromagnetic field.However, as its definition is valid for all the structures, it can be used to obtain a simplified but unified representation of the different particles.

Numerical methods and IOP
There are several methods which can simulate light scattering by nonspherical particles.These include the Discrete Dipole Approximation [26,28], Finite Difference Time Domain (FDTD) [29][30][31][32], T-matrix [33], and some other methods summarized in [27].The T-matrix method works best for particles with certain symmetrical properties.We chose the Amsterdam Discrete Dipole Approximation code (ADDA) [26] to calculate the IOP of coccoliths and coccolithophores because of its accuracy, versatility, and faster convergence for soft particles rather than the FDTD method [34].
To simulate the IOP of the coccolithophore with the ADDA code, the three dimensional space enclosing the scattering particle is divided into small grids with the grid size of λ /12 along each dimension, where λ = 0.532μm.Each grid is equivalent to a dipole and the coordinates of its center are used to pass the test described above to judge if it is occupied by the core, the carbonate shell, or empty.The refractive index is given to the dipole according to the region where it belongs.The coordinates and refractive indices of the dipoles are written into the input files of the ADDA code.Random orientation averaging is applied to the IOP calculated by the ADDA.The ADDA uses three Euler angles to define orientations of a scattering particle.Orientation averaging over one Euler angle is equivalent to rotating the scattering plane, which can be done with a single internal field calculation.Averaging over the other two Euler angles is performed by Romberg integration in an adaptive iterative mode, i.e., only a necessary number of orientation angles have to be evaluated.The termination of orientation averaging is based on estimations of relative errors of scattering cross sections.We have used the criteria of EPS=0.001 to perform the random orientation averaging, where EPS is the maximum relative error of the scattering cross sections.In general, the numbers of orientations are around 160, 120, and 220, for each size of detached coccolith, symmetrical and asymmetrical coccolithophores, respectively.Readers are referred to [26] and the reference within for details of random orientation averaging used in the ADDA code.
The ADDA package calculates the scattering matrix, which is the Mueller matrix for light scattering by a single particle [35]: where I, Q,U,V are the Stokes parameters; the subscripts i and s denote the incident and scattered light, respectively; k = 2π/λ .For a mirror symmetric particle, the scattering matrix is block diagonal after random orientation averaging, i.e., S 13 = S 14 = S 23 = S 24 = S 31 = S 32 = S 41 = S 42 = 0 [27].Strictly speaking the coccolith and coccolithophore model described in this paper do not have mirror symmetry.However, it is somewhat close to a mirror symmetric particle.The calculation result confirms that the off block diagonal elements are virtually zero (S i j /S 11 < 10 −7 ) for the detached coccoliths, and symmetrical coccolithophores.The off block diagonal elements S i j /S 11 of asymmetrical coccolithophores are smaller than 0.1 for all the scattering angles.They will not be presented in the following.In addition, we have confirmed that the simulated scattering matrix elements satisfy the following inequality for all scattering angles and all particles [36]: In the radiative transfer community, the phase function is often used to separate scattering cross section C scat from S 11 , which is defined as: (18) so that 4π pdΩ = 1 (19) Note that in the atmospheric radiative transfer community, the term phase function is used to denote p = p * 4π.We have adopted the convention in Eq. ( 19) to be consistent with the ocean optics community [6].It is also convenient to define and study the scattering efficiency in the light scattering theory: where r v is the volume equivalent radius in this application.Moreover, the term "reduced scattering matrix" is often used to denote the scattering matrix divided by S 11 .
The scattering matrix contains all the light scattering information by a particle.In reality, different applications require different forms of the information.Two important pieces of information is the backscattering probability B bp and the depolarization ratio δ L : where Θ is the scattering angle.The backscattering probability B bp has important applications in ocean optics [6,7], while δ L is important for particle categorization in lidar applications [37].Both B bp and δ L are studied in the next section.

Results
In this section we present the simulation results for detached coccoliths, symmetrical, and asymmetrical coccolithophores.Figures 3(a) -3(d) show the phase function and the reduced scattering matrix elements as a function of scattering angle for detached coccoliths.Two imaginary refractive indices m i = 0 and m i = 0.01 are used.In addition to simulation results, the reduced matrix elements of the average ocean water from Voss and Fry [14] are also plotted for comparison.Overall the differences between m = 0 and m = 0.01 are small, especially for small sizes of coccoliths.The reduced scattering matrix elements do show similar features to those of the average ocean water.However, they are different in some ways., consistent with the measurements by Voss and Fry.However, there is no exact match, which is not surprising because the measurement in [14] was done for clear ocean waters.Furthermore, the coccolith is just one of the possible components in natural water and there are a plethora of different organisms which would have to be considered such as diatoms, dinoflagellates, zooplankton, just to mention a few.Another observation is that the reduced scattering matrix elements are closer to the Rayleigh scattering matrix for smaller particles, as expected.However, they show larger deviations as size increases.Figures 4 and 5 show the same as Figs.3(a) -3(d) but for symmetrical and asymmetrical coccolithophores, respectively.Symmetrical coccolithophores show rapid resonant features, similar to spherical particles.Asymmetrical particles are smoother as a function of scattering angle, which is more realistic.Note that one could also achieve a smooth scattering matrix for symmetrical coccolithophores by taking account of the particle size distribution.The influence of refractive index is somewhat larger in these coccolithophores than that in coccoliths, mainly due to larger sizes for coccolithophores shown in Figs. 4 and 5.The impacts are mostly in the backscattering directions of Θ > π/2.Different from coccoliths, the reduced scattering matrix elements for coccolithophores show sharp peaks around backscatter.This suggests that it is difficult to extrapolate measurements towards 180 • .A similar feature has been observed in the scattering matrix elements of large aerosol particles (see Fig. 3 in [38]).It is interesting to observe the smallest values of S 22 /S 11 shift to larger scattering angles (Θ around 160 • ) for coccolithophores.This may mean that the average ocean water scattering matrix contains contributions largely from smaller particles.The contribution of large particles could be more significant for coccolithophore blooms, however.Svensen et al. have reported the measurement of the scattering matrix for EHUX (see Fig. 6 in [17]).The simulation of S 12 and S 33 is similar to what they have shown.However, it is hard to do a more quantitative comparison because the numerical values of the measurement are not available to us.In addition, they have assumed that EHUX is a spherical particle so (S 33 + S 44 )/2 is used to represent S 33 , which does not apply to our nonspherical model.Figure 6 shows the difference between S 33 /S 11 and S 44 /S 11 for symmetrical coccolithophores.In the forward scattering direction, the difference is small.However, they show large differences in the backscattering region, which is also an indicator of nonsphericity, in addition to the S 22 element.Besides the scattering matrix, it is informative to study the scattering efficiency, backscattering, and lidar depolarization ratio for coccoliths and coccolithophores.Figures 7(a) -7(c) show these three quantities as a function of the equivalent volume size parameter.It is shown that Q scat , B bp , and δ L all decrease as the imaginary refractive index increases.Spherical particles also have similar trend for Q scat [39].It looks like particles with more nonsphericity scatter light more efficiently than those particles of less nonsphericity with the same equivalent volume size parameter.Our reasoning is based on the ranking of the three particles "sym", "asym", and "detached", from the least to highest nonsphericity.The backscattering probability first decreases as x increases if particles are small (x < 10).Then it increases as x increases.After x > 35, the rate of increase decreases a bit.Overall the B bp for detached coccoliths is ranging from 0.005 to 0.025.Gordon et al. have shown that the backscattering probabilities are around 0.04 for their modified "fishing reel" model [24].Though the difference are noticeable but in the same order of magnitude mainly due to subtle different assumptions used in constructing the model.Furthermore, Gordon et al. have found an upper limit of the backscattering cross section estimated from in situ measurements, which is 0.123 ± 0.039μm 2 per coccolith at 500 nm [24].Our simulation shows that the backscattering cross section (not plotted) ranges from 0.04μm 2 to 0.07μm 2 for coccoliths with diameters smaller than 2.83μm, which is consistent with [24].The linear depolarization ratio shows irregular variation as a function of x.Asymmetrical coccolithophores do show larger δ L than symmetrical coccolithophores in general.

Conclusion
A realistic nonspherical model for EHUX is built and the optical properties are calculated using the ADDA code.First the coccolith model is built and a number of (14-21) coccoliths are grouped, rotated, and rearranged around a nearly spherical core to form a representation of the coccolithophore cell.The index of refraction is set to 1.2 for the coccoliths and 1.04 for the spherical core.Two models are built for coccolithophores: symmetrical and asymmetrical models.The model is then transformed into dipole representations which can be used as inputs to the Amsterdam Discrete Dipole Approximation code (ADDA) to calculate the complete set of IOP for this cell.Our simulation results show that the Mueller matrices of coccolithophores show different angular dependence from those of coccoliths.The scattering matrix elements are smoother as a function of scattering angle for the asymmetrical model than those of the symmetrical model.The most interesting element is S 22 , which indicates the nonsphericity of the particles.The smallest value of S 22 /S 11 of the detached coccolith is at around 90 • to 120 • , which is similar to Voss and Fry [14].This position shifts to larger scattering angles for the coccolithophore models.It is notable that our simulation shows S 33 is not equal to S 44 for coccolithophores, which is another indicator of nonsphericity.A particle with larger imaginary refractive index shows smaller scattering efficiency, backscattering probability, and depolarization ratio.

Fig. 1 .
Fig. 1.(a) The detached coccolith model with its axis shown.(b) a sketch of a cross section of the coccolith in a half plane containing the axis.(c) the projection of the inner or outer spherical cap in the x − y plane.
) and 2(b) show a realization of the symmetrical and asymmetrical coccolithophore models, respectively.The symmetrical model consists of a spherical core surrounding by a carbonate shell formed by individual coccoliths.The asymmetrical model is derived by deforming the symmetrical model with a random process.For clarity, Figs.2(a) and 2(b) show only the carbonate shell part while the spherical core is hidden.
12 S 13 S 14 S 21 S 22 S 23 S 24 S 31 S 32 S 33 S 34 S 41 S 42 S 43 S 44 The two elements S 12 /S 11 and S 33 /S 11 are closer to those of Voss.The differences in S 22 /S 11 are larger.Note that S 22 /S 11 is always 1 for spherical particles so the deviation of S 22 /S 11 from 1 is an indicator of nonsphericity.Our coccolith model is highly nonspherical which is the reason that S 22 /S 11 is smaller than 1 except for Θ = 0.The S 22 /S 11 values for different sizes of coccoliths all shows smallest values around 90 • to 120