Angle‐dependent light scattering in tissue phantoms for the case of thin bone layers with predominant forward scattering

The cochlea forms a key element of the human auditory system in the temporal bone. Damage to the cochlea continues to produce significant impairment for sensory reception of environmental stimuli. To improve this impairment, the optical cochlear implant forms a new research approach. A prerequisite for this method is to understand how light propagation, as well as scattering, reflection, and absorption, takes place within the cochlea. We offer a method to study the light distribution in the human cochlea through phantom materials which have the objective to mimic the optical behavior of bone and Monte‐Carlo simulations. The calculation of an angular distribution after scattering requires a phase function. Often approximate functions like Henyey–Greenstein, two‐term Henyey–Greenstein or Legendre polynomial decompositions are used as phase function. An alternative is to exactly calculate a Mie distribution for each scattering event. This method provides a better fit to the data measured in this work.


| INTRODUCTION
Those affected by sensorineural hearing loss can have some of their hearing restored with the help of an electrical cochlear implant (eCI) [1][2][3].This involves transmitting acoustic environmental stimuli to the brain by means of electrical stimulation of the spiral ganglion neurons (SGNs).However eCIs are limited in their function by a wide current spread resulting in a low resolution of the spectral coding of auditory information [3,4].A new method to avoid this limitation is an optical cochlear implant (oCI) [2,5,6].Here, instead of electrodes, light emitters are placed along the tonotopic axis of the cochlea to stimulate the SGNs, which were made light sensitive through Optogenetics [7].Light promises to have a higher spatial resolution than electrical excitation, thus increasing the spectral selectivity of artificial sound encoding [8,9].If this concept proves to be feasible in clinical applications, it could facilitate speech recognition in background noise as well as the understanding of music and melodies [2,10].For further development of the oCIs and to demonstrate the higher spatial resolution of light, we need to understand how light is scattered in the human cochlea.Considering the size and accessibility of the auditory system, the use of tissue-simulating objects to mimic the properties of human tissues can help to understand scattering and absorption behavior of the human cochlea tissue [11][12][13][14][15][16][17].These so-called "phantoms" have to be tuned, to match the desired tissue properties, and can then be used for a number of applications like initially testing system designs or comparing performance between systems [18].For our research we want a phantom that mimics the optical properties within the human cochlea and is suitable for experimental and simulative investigation [18][19][20].In the following, we will discuss a Monte-Carlo approach to simulate the spread of light of micro emitters through tissue phantoms and contrast it with an experimental method to understand the propagation of light through thin bone layers.These phantoms consist of a filler matrix and a proportion of homogeneous, spherical, polymer scatterers.For our simulative approach, we combine a wave optical method, to calculate the Mie-scattering [21,22] behavior on a single polymer sphere, with the ray optical approach of the Monte-Carlo simulation [23,24].In this way, we calculate the angle distribution of one scatterer which we utilize as the phase function for scattering events in the Monte-Carlo simulation.The phase function describes the spatial distribution of scattered light for a single scattering event.In both approaches, we used a red-light source which resembles the actual emitters that are used in oCIs.Sources of that wavelength range have the advantage of better neural stimulation with less tissue scattering and a lower risk of phototoxicity [4].Calculating the exact angle distribution for each scattering event has proven to yield a higher agreement with measurements of phantoms than a simulation of those layers by using the Henyey-Greenstein formula as approximated phase function.

| Tissue phantoms
To create a phantom that behaves close to the actual tissue when exposed to radiation, it is important to know the main optical characteristics of the tissue that interacts with light [18].For that we want to match the scattering coefficient μ s λ ð Þ, the absorption coefficient μ a λ ð Þ, and the anisotropy factor g λ ð Þ (see Section 2. 3) [18], which represent the intensity-weighted mean value of the cosine of the scattering angle θ and indicate the directionality of the scattering.Considering the fact that the shape of the phase function is usually not known, the angle scattering distribution is often characterized by this single parameter, the anisotropy factor g [11,17].Since light within the cochlea will interact mainly with the bone tissue of this structure, it is important to find a phantom that approximates the corresponding scattering characteristics.For this purpose, it has proven successful to regard the bone of the cochlea as a layer of spherical, homogeneous scatterers.The known optical parameters of bone are listed in Table 1.Polymer spheres [25], which are placed in an epoxy [26] matrix with a packing fraction η, are used as phantom material to mimic the optical properties of bone according to the parameters of Table 1 (see Refs. [18,27,28]).The optical parameters of those polymer scatterers are calculated as described in Section 2.2 and are listed in Table 2.As an idealization, we neglect the absorption in the phantoms (μ a λ ð Þ ¼ 0) (see Table 2).In Figure 1, we simulated a bone layer in intercellular fluid and the resulting scattering angle distribution with the optical parameters for μ a and μ s of Pifferi et al. [15] and a second simulation with an idealized μ a ¼ 0:0 cm À1 .One can see, that such small absorption coefficients have almost no impact on the scattering behaviour, so the absorption is neglectable here.The handling of the anisotropy factor g is explained in Section 2.3.By adjusting the packing fraction it is possible to fit the scattering properties (μ s ) for phantoms (Table 2) to the values of bone given in Table 1.This is essential to calculate the scattering coefficient and the mean free path for the Monte-Carlo simulation (see Section 2.2).

| Scattering coefficient
Starting from the scattering efficiency of one polymer scatterer in epoxy, see Appendix A.1, we calculate the effective cross section σ.Combined with a density of scatterers with a volume packing fraction η and particle radius R, we calculate the scattering coefficient and the mean free path between two scattering events Together with the angle-dependent intensity distribution of unpolarized light for one scatterer (see Appendix A.1), which is used as the phase function for the Monte-Carlo simulation we can define a tissue layer with μ s , l s and with no need for an anisotropy factor g, since the shape of the phase function is defined through S U θ ð Þ.Alternatively, as a phase function, the Henyey-Greenstein function (HGS) is often used as an approximation for the shape of an actual phase function for the scattering of radiation by a particle.

| Henyey-Greenstein function
In the case of using the HGS, the anisotropy factor g can be calculated exactly from the intensity distribution S U θ ð Þ for a single scattering event A g-value of g ¼ 1 equals total forward scattering, g ¼ 0 means isotropic scattering and g ¼ À1 total backward scattering.We describe scattering of a single, homogeneous sphere with the angle-dependent intensity distribution This gives us, together with the previously explained Mie distribution, see Appendix (Appendix A.1), two possible phase functions to define a scattering event in the Monte-Carlo simulation.The first one uses the Mie distribution of light scattering on one sphere, that is, Equation (4).The later uses the HGS ( 6) with the addition of the anisotropy factor g calculated with (4) and ( 5).In Figure 2, one can see the one-term and two-term HGS compared with the scattering pattern of a single polymer sphere (S U θ ð Þ).Since the two-term HGS overfits the exact calculated angle distribution more than the one-term HGS, we will use the one-term HGS as phase function in order to make a closer comparison between the two phase functions.

| Monte-Carlo simulation
To describe a layer of tissue with multiple scattering events, we can combine the previously introduced waveoptical description for a single scattering event with a Simulated scattering angle distribution of a 100 μm bone layer for a wave, unpolarized to the plane of measurement with the known parameters of Pifferi et al. [15] (orange curve) and with an idealized μ a ¼ 0:0 cm À1 (blue curve).
Monte-Carlo simulation.These simulations for the light propagation require a free path length of the photons, scattering angles, as well as reflection or transmission at the boundaries of the layer.The probability that, in the Monte-Carlo simulation, a photon is scattered within a distance smaller than s is expressed by the probability distribution function P S < s ð Þ. S is a random variable representing the distance that a photon travels between two scattering events.s is a value in the definition range of S, which is greater than zero [29].
P s ð Þ is the probability distribution function.In contrast f s ð Þ is the probability density function, which expresses the probability that a photon is scattered while traveling between s and s þ ds.The integral of each density function f s ð Þ is equal to the probability distribution function P s ð Þ [29], A random number a 0, 1 ½ with a ¼ P s ð Þ is assigned to a random variable s, such as the path length of a photon.Equation ( 8) is used to obtain a series of values for a, which specify a series of values s.The histogram of the s-values corresponds to the probability density function f s ð Þ.For the described problem, a steady-state simulation will be used.The aim is to find the intensity distribution after scattering processes.
To define a tissue layer in the Monte-Carlo simulation with the Mie distribution as phase function (Equation 4), we first simulate a scattering event of an incoming plane wave on one polymer sphere in epoxy.In Figure 3, we show the angle distribution p θ ð Þ for scattering events on different-sized spheres, ranging from 50 to 3000 nm in diameter.With increasing size of the scatterers, the polarization of the source becomes less important.Furthermore we see strong backward scattering for small spheres (50, 100 nm) and increasing forward scattering with increased sphere size.Since we want to mimic the scattering behavior in bone tissue, we aim for a diameter of the scatterers of 1000 nm.With this diameter the scattering distribution shows the Mie scattering regime and does not change when we look at different polarizations (see Figure 3).Therefore a scattering event in the Monte-Carlo simulation will result in an angle distribution, which corresponds to the distribution on a single sphere with d ¼ 1000 nm as shown in Figure 3.This distribution is calculated with Equation (4) which results out of the Equations (A19) (see Appendix A.1).The second step is to define a layer containing scatterers with radius r and calculating the packing fraction dependent mean free path between two scattering events and the scattering Scattering angle distribution on polymer spheres with a diameter between 50 and 3000 nm for a wave (λ ¼ 635 nm), perpendicular, parallel, and unpolarized to the plane of measurement.
F I G U R E 2 One-term Henyey-Greenstein function (HGS; blue curve) and two-term HGS (orange curve) compared with the Mie scattering S U θ ð Þ of one polymer sphere (green curve).
coefficient (see Section 2.2).Thus, the layer is defined with those parameters for the scatterer with the radius R. If we use Henyey-Greenstein as a phase function, the Monte-Carlo simulation additionally uses the computed anisotropy parameter g for a scattering event.In Figure 4, we can see how the mean free path depends on the size of the scatterers.Considering the density, l s is proportional to 1/R in the range of Rayleigh scattering.In this range the scattering efficiency will decrease and l s will increase for small radii compared with the wavelength (R ( λ), because the cross section of the Rayleigh scattering is.
with σ Th as the Thomson effective cross section.As the radius increase, the scattering efficiency and μ s increases, which in turn reduces the mean free path length due to l s ¼ 1=μ s .Since the human cochlea consists of very thin bone layers, we aim to generate a strong forward scattering via small particles (see Table 2) and thin layers of polymer and epoxy mixtures.This will result in very few scattering events per photon which will contribute to the expected strong forward scattering of thin bone layers.We therefore use spheres, which have the same order of magnitude in diameter as the used wavelength of the light source.The light source for each simulation run is defined as a beam consisting of 10 7 photons with a wavelength of λ ¼ 635 nm.We assume the layer as infinitely extending in the x and y directions, therefore the starting position of the photons is of no importance.The layer thickness is varied between 100 and 500 μm.Like previously mentioned, the spheres are placed in the medium with a packing fraction η.We will use η ¼ 0:00625 and η ¼ 0:01015 for a low number of scattering events.With those packing fractions, we calculate the values for μ s and g given in Table 2.
The location and angle of every photon is saved when it reaches a border of the simulation box.This results in an angle distribution p θ ð Þ which corresponds to the scattering intensity I θ ð Þ with the scattering angle θ of the photon beam through the layer.To be able to make comparisons to measurements, we normalized the angle distribution to the maximum of the values of the distribution.

| Experimental approach
For sample preparation, the polymer spheres were mixed with a packing fraction of η ¼ 0:00625 or η ¼ 0:01015 to match the scattering parameters expected of bone from Table 1 with the epoxy medium for 5 min and then applied to the first glass plate with m ¼ 1:515, x 1 ¼ 18 mm (see Figure 5).The mixture was then hardened in vacuum for 6 h.After that, the second glass plate was placed on the sample followed by completely drying the sample in vacuum.This procedure prevented the formation of air bubbles during both mixing and placement of the glass plates (see Refs. [18,27]).The layer thickness of those samples is controlled via spacers which are manufactured according to the DIN 988 [30] standard and is, like in the simulation, varied between 100 and 500 m.A schematic structure of a sample is shown in Figure 6.
The experimental setup shown in Figure 7, containing the described samples, the laser diode CPS635S [31] which has a wavelength of 635 nm, an elliptical, collimated output beam (3.8 mm Â 12 mm) and a beam divergence (full-angle, 1/e 2 beam width) of 1.6 mrad, and the light sensors TSL2591 [32]  Carlo simulation does not include any glass layers but is compared with our measurements we included a Snellius correction in the simulation, with α and β being the incoming and outgoing angle and m matrix the refractive index of the surrounding epoxy medium.We can utilize the unpolarized transmission coefficient T u as a threshold as a photon will get reflected with a predetermined probability.For that we generate a random number z 0, 1 ½ .If z is smaller than T u the corresponding photon will get transmitted.If not, the photon will get reflected.With this method, we can follow the reflection processes of each photon until it has actually left the layer.This compares favorably to a correction of the experimental data where we can only correct the first time a photon tries to exit the layer.

| COMPARISON
As shown in Figure 8, we first modeled a scattering event on a single polymer scatterer.We use the calculations from Appendix A.1, which are implemented in a Mie scattering module (PyMieScatt) [34].The angle distribution shows the typical Mie scattering behaviour with a large amount of forward scattering, which is ideal for a layer of spherical scatterers.Therefore, the angle distribution is used as a phase function to define a scattering event in the Monte-Carlo simulation.Together with the scattering coefficient μ s ¼ 33:57 cm À1 , the packing-fraction-dependent mean free path l s ¼ 297:80μm, we obtain the angle distribution for η ¼ 0:01015 and d ¼ 500 μm shown in Figure 9 (orange curve).For these parameters, an average of 3:51 scattering events occur and 18:5% of the photons remain unscattered.Compared with the scattering intensity of a single sphere, this results in a wider distribution, which is still mainly forward scattered.For the comparison in Figure 9, we excluded unscattered photons which make F I G U R E 6 Schematic structure of a sample.
up a significant portion of the angle distribution.The blue curve in Figure 9 is the same distribution as shown in Figure 8 and the small blue curve is the backscattered light from said distribution.The same layer simulated with Henyey-Greenstein as phase function (see Figure 10 and Equation ( 6)) and an additional anisotropy factor of g ¼ 0:95 (brown curve in Figure 10) results in the green curve.This distribution shows a more peaked pattern with less backscattering.If we compare theory and experiment, we see the angle distributions are shown in Figure 12.Since the polarization of the incident light has no influence on the angular distribution in the prevailing size ratio of wavelength and sphere diameter, the unpolarized portion of the light is compared.In the measurements, the unscattered light is included showing a sharp distribution for small angles.The Monte-Carlo simulation that uses the Mie distribution from one sphere as scattering function (blue curve) shows a strong agreement with the experimental data (green curve) for larger angles.The second Monte-Carlo simulation which uses the HGS as scattering function (orange curve) produces lower values than the measured data.All three approaches show a strong forward-directed scattering behavior with a large proportion of unscattered light (compare Figure 9).Since we broadened the unscattered light of the simulations with a Gaussian fit to a zero measurement of the laser through an epoxy layer (see Figure 11), comparisons of the curve form between simulations and measurements in the range of 0 ∘ À 14 ∘ are not meaningful.Nevertheless, the proportion of unscattered light can be detected and compared in this area.
The discrepancy between the simulation methods at small angles decreases with reducing the layer thickness.
In Figure 13, we consider a thinner layer, which is d ¼ 200 μm thick with a packing fraction of η ¼ 0:01015, μ s ¼ 33:57cm À1 , and l s ¼ 297,80μm.Here, the agreement between both simulation methods is very high for small angles.For larger angles, the Monte-Carlo simulation with the Mie distribution also has a higher agreement with the experimental data for this thinner layer than with the HGS.For layers with few scattering events simulations with an exactly calculated distribution as phase function seems to better represent experimental data for large angles than simulations that use the Henyey-Greenstein as phase function.

| SUMMARY AND CONCLUSIONS
We investigated the scattering properties of human tissue by using polymer-epoxy layers as phantoms.On the one hand, this was done with a Monte-Carlo simulation with which we showed that for a description of scattering in tissue, the HGS is not necessary, but the scattering distribution of a scatterer as a phase function provides an accurate description of the behavior of light in a tissue layer.
Together with the mean free path l s and its corresponding scattering coefficient μ s , which are calculated from one scattering event, we are able to define tissue layers without using the Henyey-Greenstein, or another approximated phase function.This is supported by an experimental investigation in which the agreement between the simulative method and the corresponding measurements of the simulated layers is shown.For layers with few scattering events, the Mie distribution as phase functions shows higher agreement with the experimental data than the HGS for large angles.Therefore calculating the exact phase function for every scattering event has proven to provide a more realistic scattering behaviour in tissue layers than the often used HGS which is one of any function that meets the condition of normalization 2π and anisotropy 2π To find an analytic phase function that is more consistent with the exact Mie distribution, one could search for a function that fits the Mie angle distribution and therefore meets the said conditions.This has the advantage that the phase function is not just an assumption but is based on a calculated distribution.In further work, we will provide this kind of function together with an analytical model, thereby calculating the angle distributions and transmissions after any number of scattering events.Furthermore, we want to verify the scattering pattern, and therefore applicability, of used phantoms with the scattering pattern of real cochlea bone tissue and compare the phase functions and their corresponding angle distributions to the measurements of real cochlea bone.To improve the presented approach and make it come closer to real cochlea bone one should consider multiple sphere sizes in the layers.For this purpose, we will first produce appropriate phantoms in subsequent work and perform corresponding measurements and simulations based on the results of the measurements of real cochlea bone.In addition, we will investigate the scalability of phantom layers to create structures that have the shape and optical properties of a human cochlea.This will allow further investigation of the scattering behaviour of an emitter array in the human cochlea.

APPENDIX A A.1 | MIE SCATTERING THEORY
The interaction of an electromagnetic wave with spherical, homogeneous scatterers is the subject of Lorenz-Mie theory [35].If one considers a spherical scatterer in isolation, taking this theory into account, scattering parameters for a single scattering event within a layer can be calculated on this sphere.The incident light can be described as a plane wave polarized in the x-direction which gets scattered into a wave Hankel functions of the first kind (superscript (3)).The indices e and o mean even/odd (symmetrical/antisymmetrical).The scattered waves overlap with the part of the plane wave that remains unscattered.The coefficients are the amplitudes of the partial waves that make up the scattered wave [35].They consist of the Ricatti-Bessel functions ζ x ð Þ and ψ x ð Þ which depend on the ratio of particle circumference to wavelength as well as the complex refractive index and are derived from the Hankel-Bessel functions.From the coefficients a n , b n , we calculate the polarizationdependent scattering intensities and the scattering efficiency U stands for "unpolarized."The factors π and τ are calculated as follows and VEML7700 [33].The sensors were mounted on an arm which rotates around the sample in a range of 1 ∘ 230 ∘ .Measurements of the scattering intensity I θ ð Þ with the scattering angle θ are done in steps in air behind the layer.Since the Monte-F I G U R E 4 Double logarithmic plotted mean free path of photons versus radius R of the polymer spheres in the simulated layer for η ¼ 0:01015 and λ ¼ 635 nm.F I G U R E 5 Dimensions of borosilicate glass and spacers for phantom samples with x 1 ¼ 18mm, x 2 ¼ 25mm, and d ¼ 0:1 mm, d ¼ 0:2 mm, d ¼ 0:3 mm, or d ¼ 0:5 mm.

F
I G U R E 7 Schematic structure of the experimental setup.F I G U R E 8 Scattering intensity I θ ð Þ for one polymer sphere in epoxy resulting from an incoming wave which is unpolarized to the direction of propagation.F I G U R E 9 Logarithmically plotted scattering intensity I θ ð Þ with Mie scattering as phase function (orange curve) and with Henyey-Greenstein function (HGS) as phase function (green curve) of a simulated d ¼ 500 μm polymer-epoxy layer and for one polymer sphere in epoxy (blue curve).

F I G U R E 1 0F I G U R E 1 2
Henyey-Greenstein function for multiple gvalues.Logarithmically plotted measured scattering intensity (green curve) of a d ¼ 500 μm polymer-epoxy layer compared with the Monte-Carlo simulation with Mie distribution (blue curve) and the Monte-Carlo simulation with Henyey-Greenstein function (HGS) for the corresponding parameter set.F I G U R E 1 1 Logarithmically plotted measured scattering intensity of an d ¼ 100 μm epoxy layer.F I G U R 1 Logarithmically plotted measured scattering intensity (green curve) of a d ¼ 200 μm polymer-epoxy layer compared with the Monte-Carlo simulation with Mie distribution (blue curve) and the Monte-Carlo simulation with Henyey-Greenstein function (HGS) for the corresponding parameter set.
Values of optical parameters of polymer spheres at a wavelength of λ ¼ 635nm and a packing fraction η in epoxy.