Characterization of ejecta in shock experiments with multiple light scattering

Upon impact, the free surface of a solid metal may eject a cloud of fast and fine particles. Photon Doppler Velocimetry (PDV) is one of the optical diagnostics used to characterize these ejecta. Although the technique provides a direct way to estimate the particle velocities in the single scattering regime, it has been shown that multiple scattering cannot be neglected in real ejecta. Here we derive a model for PDV measurements starting from first principles of wave scattering. We establish rigorously the relationship between the specific intensity and the measured signal, as well as the radiative transport equation (RTE) that describes the evolution of the specific intensity upon scattering and absorption in a dynamic ejecta, including the effects of inelastic scattering and inhomogenities in the optical properties. We also establish rigorously the connection between the Monte-Carlo scheme used for numerical simulations and the solution to the RTE. Using numerical simulations, we demonstrate the crucial role of multiple scattering and inhomogeneities in the particle density and size-velocity distribution. These results could substantially impact the analysis of ejecta by PDV.


I. INTRODUCTION
The extreme conditions of shock-compression experiments make it possible to observe unique physical phenomena.Although they have specific applications, for example in inertial fusion [1], from a purely physical point of view, these regimes are worth investigating in their own right.Under the correct impact conditions, shockwaves inside a material cause it to partially melt and its free surface to eject a cloud of particles called ejecta.It has been shown that ejecta creation is mainly due to surface imperfections of machined materials [2,3].Since then, ejecta formation has been intensively studied [4] with the goal of predicting both the total amount of ejected mass and the size-velocity distribution of the ejecta.This effort was twofold.
On the ejecta modelling side, there have been advances in theory [1,5] showing that ejecta are a limiting case of Richtmyer-Meshkov instabilities.Numerically, continuum simulations [6] and more recently molecular dynamics simulations [7][8][9] have made it possible to determine expected size-velocity distributions.Parallel to these advances, a large number of diagnostics were developed and used in experiments to obtain as much information as possible on the ejecta [4,10].Among other optical diagnostics, Photon Doppler Velocimetry (PDV) has played a key role thanks to its ability to simultaneously measure the velocities of several targets, and to the trade-off between a fine temporal resolution and a broad range of accessible velocities [11,12].However, current analysis of PDV signals relies on the single scattering hypothesis (light collected is assumed to be scattered only once in the medium), which provides a direct relationship between the measured Doppler shifts and the velocity of the scatterers.Since in practice the thickness of the ejecta can largely exceed the photon scattering mean-free path, accounting for multiple scattering is of utmost importance for a quantitative analysis of the experiments.
In the last decade, several models have been developed to account for multiple scattering [13][14][15][16][17][18][19][20].A key point is to assume that the PDV spectrograms correspond to the specific intensity detected by the probe, and that the specific intensity obeys the Radiative Transfer Equation (RTE) [21].An advantage of the RTE is its ability to deal with complex particle size and velocity distributions, while being solvable numerically in realistic geometries.The use of the RTE raises several basic issues that remains to be clarified.Although a phenomenological approach is widely used to introduce the specific intensity and the RTE, a rigorous framework based on the wave equation and a statistical description of the scattering medium is available [22,23].This theoretical framework should allow one to establish rigorously the relationship between the specific intensity solution to the RTE and the spectrogram deduced from the PDV signal.Moreover, the RTE that is used is usually modified to account for Doppler shifts and inhomogeneities in the scattering properties of the medium.This generalized form of the RTE also deserves to be derived in more than an intuitive way.Finally, Monte-Carlo simulations are widely used [14,18,20], but the fact that they provide solutions to the (generalized) RTE is often overlooked.Proving that the Monte-Carlo approach amounts to solving the RTE would fill a gap in the chain connecting simulations to experiments.
The purpose of this work is to answer these questions, in order to provide a complete theory supporting the analysis of PDV experiments in ejecta.To proceed, we first establish a rigorous connection between the PDV signal and the specific intensity scattered by the ejecta and collected by the probe.Then, we derive the RTE from the wave equation, taking into account all important features of shock ejecta such as the size and velocity distributions of the scatterers, and heterogeneities in the number density.Finally, we establish a random walk scheme that naturally supports the connection between Monte-Carlo simulations and the solution to the RTE.The model is used to demonstrate unambiguously the importance of multiple scattering in shock ejecta.
The paper is organized as follows.In Sec.II, we describe a typical shock experiment including the PDV setup used to characterize the ejecta.We also explain how a PDV spectrogram can estimate accurately the particle velocity distribution under the assumption of single scattering of the probe light.Section III is dedicated to the analysis of the multiple scattering regime, and the link between the spectrograms computed from PDV measurements and the specific intensity of the collected light.The derivation of the RTE describing multiple scattering of light scattering in ejecta is reported in Sec.IV.Finally, the justification of the Monte-Carlo model is presented in Sec.V, together with some applications to the analysis of real spectrograms.

II. SHOCK EXPERIMENTS AND PHOTON DOPPLER VELOCIMETRY IN THE SINGLE SCATTERING REGIME
In this section, we recall some standard results for the characterization of shock experiments using PDV in the single scattering regime.In a typical shock experiment, a shockwave is released by impact in a material.Under this extreme excitation, the material partially melts and the shockwave interacting with the irregularities of the free surface generates microjets of matter.These micro-jets extend before fragmenting and forming a cloud of fast and fine particles, referred to as the ejecta [see Fig. 1 (a)].
Ejecta are usually characterized optically using interferometric techniques such as a PDV [11,12].The usual setup is a Michelson interferometer, in which the reference beam interferes with the scattered light coming from a measuring arm, as represented schematically in Fig. 1 (b).Since the particles and the free surface are moving, the scattered field is slightly Doppler shifted in frequency compared to the light in the reference arm.This creates a beating signal on a photodiode detector, that allows one to directly measure the Doppler spectrum of the scattered field.In practice, a time-frequency analysis of the signal is performed, resulting in a spectrogram as that shown in Fig. 1 (c).The usual representation directly transforms frequency shifts into particle velocities, under the fundamental assumption of single scattering.As a result, a vertical line of the spectrogram is understood as the velocity distribution of the particles at a given time.
Let us briefly summarize the main steps leading to a the-oretical description of the spectrogram in the single scattering regime.The measured signal corresponds to the intensity , where Ē(t) is the analytical signal associated to the real-valued field E(t) describing the light received on the detector.Focusing on the time dependence, the analytic field is defined as where E(ω) is the Fourier transform of E(t).We use a scalar description of the light field, leaving aside polarization effects.Considering that the illuminating field is monochromatic at frequency ω 0 , the associated analytical signal is Ē0 (t) = A 0 exp(−iω 0 t).The analytical signal associated to the scattered field can be written where N(t) is the number of particles, A j (t) is the amplitude of the field scattered by particle number j, and δ ω j (t) is the corresponding Doppler shift.These quantities evolve on a time scale T c , which corresponds to the characteristic time of changes in the statistical properties of the cloud of particles.The field also contains the time scale δ T = 2π/δ ω, where δ ω is the typical Doppler shift.We note that these time scales are expected to be much larger than the period T 0 = 2π/ω 0 of the incident light field.The intensity received by the photodiode is where the superscript * denotes the complex conjugate.We have dropped the term | Ēs (t)| 2 since the amplitude of the scat-tered field is much smaller than the amplitude of the illuminating field.The first term |A 0 | 2 is a DC component that can be ignored, and the second term defines the real-valued useful detected signal Using Eq. ( 2), we immediately find that the detected signal is In practice, a Short Term Fourier Transform (STFT) is performed on the signal, defining the spectrogram as where w(t) is a time window function centered at t = 0. Various shapes [24] can be used for w(t), the key parameter being its duration T w that is chosen to clearly separate the time scales T c and δ T characterizing the ejecta dynamics.To this aim, the inequalities T 0 ≪ T d ≪ δ T ≪ T w ≪ T c need to be satisfied, where T d is the time scale associated to the digitizer bandwidth.Typical orders of magnitudes for usual experiments are T 0 = 10 −15 s, T d = 10 −10 s, δ T = 10 −9 s, T w = 10 −8 s and T c = 10 −6 s.For illustration purposes, we consider here a rectangular function for w(t), with unit amplitude and width T w .This choice is not restrictive and allows us to describe the main features of PDV spectrograms.In these conditions, from Eqs. ( 5) and ( 6), we can show that where sinc(x) = sin(x)/x.In addition, since T w is large compared to T 0 and δ T , the sinc functions can be approximated by Dirac delta functions, leading to In the single scattering regime, a mapping between the frequency shifts and the particle velocities is obtained by writing that for a scattering event occuring on particle j with velocity v j (t), the Doppler shift [25] is In this expression k s and k i are the scattered and incoming wavevectors respectively, such that |k s | = |k i | = k 0 = ω 0 /c = 2π/λ , with λ the incident wavelength.In the configuration presented in Fig. 1 (b), where the backscattered light is detected in direction k s = −k i , the Doppler shift is simply δ ω j (t) = 4πv j (t)/λ .Inserting this expression into Eq.( 8), we immediately see that in this configuration we only have access to the absolute value v j (t) of the particle velocities.The analysis above relies on the assumption of single scattering.In a statistically homogeneous cloud with longitudinal size L, this assumption holds as long as the optical thickness b = L/ℓ s ≪ 1, where ℓ s is the scattering mean free path.In an inhomogeneous cloud, the criterion becomes b = dz ℓ s (z) ≪ 1 (10) where z denotes the coordinate along the direction perpendicular to the free surface.In practice, many situations of interest rather correspond to b > 1, and even b ≫ 1 [10,14,26].For example, we will describe in Sec.V a situation with b = 42, corresponding to a real experiment.Obviously, analyzing experiments in the deep multiple scattering regime with the formalism described above is irrelevant, and leads to artefacts such as apparent velocities smaller than the free surface velocity [14], as seen in Fig. 1 (c).Establishing a rigorous treatment of spectrograms in the multiple scattering regime is the purpose of the next sections.

III. SPECTROGRAM AND SPECIFIC INTENSITY
The RTE appears as a natural tool to describe light scattering in ejecta in the presence of multiple scattering.Since the RTE is a transport equation for the specific intensity, we will start by establishing a rigorous relationship between the specific intensity and the spectrogram.The key step in defining the spectrogram from the measured PDV signal in Eq. ( 6) is a time analysis, and we will focus on the time dependence and omit the space dependence of the field in the notations in the beginning of this section.
For a rectangular time window function with width T w and unit amplitude, the spectrogram in Eq. ( 6) becomes In statistical optics, the specific intensity is defined as the time and space Wigner transform of the field [22,23,27,28].Focusing on the time dependence, the specific intensity of the scattered field is where the notation • • • denotes a statistical average over an ensemble of realizations of the scattering medium.Since we no longer assume a single scattering regime, we cannot use Eqs.( 2) and ( 5) anymore.Nonetheless, the definition given for I (t) in Eq. ( 4) remains valid.To link both expressions, we start by introducing the Fourier transform of I (t) in Eq. ( 11), and rewrite the spectrogram in the form The frequency Ω associated to the time variable t encodes the slow variations of I on the scale of T c .Since T w ≪ T c we have ΩT w ≪ 1, and we can simplify the expression which be-comes Transforming back to the time domain for I (t), we find that The gating window having a width T w large compared to δ T and T 0 and small compared to T c , the second integral can be replaced by a statistical average over the configurations of the disordered cloud by invoking ergodicity.This leads to Since δ T ≪ T w , the sinc 2 function can be replaced by a Dirac delta function, which leads to (17) Making use of Eq. ( 4) and keeping only the terms corresponding to frequencies on the order of δ ω, the spectrogram becomes with I s (t, ω) the specific intensity defined in Eq. ( 12).This expression connects the spectrogram to the specific intensity.
The full expression also exhibits a dependence on space variables (that were not used in the above derivation).Indeed, the specific intensity I s (r, u,t, ω) is defined by the relation where k R is the real part of the effective wavevector in the scattering medium and u is a unit vector.The meaning of this definition, and the fact that I s (r, u,t, ω) is the solution to the RTE, will be clarified in Sec.IV.In the large wavelength limit, it can be shown that I s is always positive and can be interpreted as the radiative flux at position r, in direction u, at time t and at frequency ω [29].From this expression, the full spectrogram finally takes the form Here G is the etendue of detection of the photodiode (surface of the photodiode and angular aperture), du means integration over the solid angle, and n is the normal to the detector surface.In this expression, the specific intensity for the scattered field can be replaced by the specific intensity for the full field I(r, u,t, ω 0 ± ω) since the incident field E 0 does not play any role for detection in the backward direction.It is given by Equation ( 20) establishes a rigorous connection between the spectrogram and the specific intensity that is often used [13,15,17,19].In particular, we note that the spectrogram at ω is the sum of the Doppler contributions at both +ω and −ω.In the next section, we derive the transport equation satisfied by the specific intensity in a scattering and absorbing medium, accounting for both inelastic scattering and inhomogeneities in the optical properties.

IV. INELASTIC AND INHOMOGENEOUS RADIATIVE TRANSFER EQUATION
The RTE is a convenient tool to describe multiple scattering of light, due to its ability to handle complex geometries, and different transport regimes from single to deep multiple scattering.Historically, the RTE was established based on a phenomenological energy balance [21].Its derivation from the wave equation (first principles) is also well-known, and requires the definition of the specific intensity in statistical optics [Eq.( 21)], together with a theory of multiple scattering [22,23].In order to account for the motion of scatterers in an ejecta, and for spatial and temporal statistical heterogeneities in its optical properties, a generalized form of RTE is needed.The purpose of this section is to derive this generalized RTE from first principles, following the approach in Ref. 30 to account for Doppler shifts, and drawing inspiration from Ref. 31 to handle statistical inhomogeneities in the medium (space and time-dependent particle number density, size and velocity distribution).

A. Bethe-Salpeter equation and statistical quasi-homogeneity
The specific intensity, introduced in Eq. ( 21), is based on the field-field correlation function, which is known to obey the Bethe-Salpeter equation.The latter is the starting point to derive the generalized RTE.It can be derived from first principles using the multiple scattering theory [23,28,32].We simply recall its formulation here.In the multiple scattering regime, light depolarizes after propagation over a distance on the order of the scattering mean-free path [33].For this reason, we consider the scalar approximation.Thus the Bethe-Salpeter is given by We note that the field-field correlation function is written for the real field E here while the specific intensity involved the analytic field.We will rigorously make the switch to Ē at the appropriate moment.Also, Eq. ( 22) is written in the time domain, a formulation which is required to take into account the quasi-homogeneous aspect of the statistical properties of the medium.In this equation, G is the average Green's function.
Physically, G(r, r ′ ,t,t ′ ) corresponds to the average electric field generated at point r and time t by a point source at point r ′ emitting an infinitely short pulse at time t ′ .We note that G is not translationally invariant in a statistically inhomogeneous medium.Γ is known as the irreducible vertex, that contains the contributions of all scattering sequences connecting the field-field correlations functions at different points and times.The Bethe-Salpeter equation can therefore be seen as a propagation equation for the field-field correlation function.
To derive the quasi-homogeneous RTE, we first assume that the statistical properties of the medium (typically the number density of particles, their size or their velocity distribution) are expected to vary slowly in space and time compared to all other characteristic lengths and times.We perform a change of variables first in the average Green function that reads where r 0 = (r + r ′ )/2 and t 0 = (t + t ′ )/2 respectively.This change of variable is standard in coherence theory [29].In the following, r 0 and t 0 will be kept fixed in particular when performing Fourier transforms thanks to the slow space and time evolution of the statistical properties of the medium.We perform a similar operation for the intensity vertex that becomes The Bethe-Salpeter equation can be rewritten in the form This transformation in G and Γ is central in the so-called quasi-homogeneous approximation.Indeed, at fixed r 0 and t 0 , this amounts to considering the statistical properties of the medium as homogeneous.Taking the spatio-temporal Fourier transform of Eq. ( 25), leads to where we have dropped the product of average fields in the right-hand side that would contribute as a source term.The next step of the derivation consists in performing a second change of variables with . Thus, Eq. ( 26) becomes where (28) To go further, the average Green function Ĝ needs to be explicitly given.It is driven by the Dyson equation that reads [23,28,34,35] The self energy Σ is the operator that accounts for all the possible scattering sequences connecting the average field at different points and times.G 0 is the Green function for the homogeneous medium without the scattering particles (i.e., air).It corresponds to the electric field solution of the wave equation with a Dirac delta source term, and an outgoing radiation condition.It is given by We perform again the change of variables defining a central position r 0 and a central time To solve Eq. ( 31), we take its Fourier transform and knowing that We now replace the average Green function by its expression given by Eq. ( 32) in Eq. (27).Making use of the relation (AB At this point, we have only performed an approximation of quasi-homogeneity of the statistical properties of the medium and additional approximations are required in order to obtain the RTE.

B. Weak extinction limit and radiative transfer equation
We now perform the large-scale approximation |q| ≪ {|k|, |k ′ |} which is also known as the radiative transfer limit [23,36,37].It assumes that the field-field correlation function varies on length scales 2π/|q| much larger than the effective wavelength in the medium 2π/k R , the expression of which will be given below.Since the typical distance over which the field-field correlation function varies is the extinction mean-free path ℓ e , the large-scale limit holds as soon as k R ℓ e ≫ 1.This regime of weak extinction is relevant to ejecta in shock experiments.A large scale asymptotics for the time dependence is also relevant, and we assume |Ω| ≪ {|ω|, |ω ′ |}.Assuming non-resonant scatterers and weak non-locality, the q and Ω dependencies can be dropped in the average Green function Ĝ , the self-energy Σ and the intensity vertex Γ in Eq. (33).However, in order to obtain the correct expression for the transport velocity, it is crucial to keep the first term of the Taylor series in Ω for Re Σ.We get The weak extinction regime also amounts to assuming | Σ| ≪ k 2 0 .Thus the average Green function Ĝ in Eq. ( 32) is very peaked around |k| = k 0 [23].In this case, the self-energy Σ can be evaluated on-shell at |k| = k 0 .By making use of the identity lim where PV is the Cauchy principal value, we find that Im Ĝ (r 0 , k,t 0 , ω) = πδ k 2 − k 2 0 − Re Σ (r 0 , k 0 ,t 0 , ω) (36) where we have assumed that the disorder is statistically isotropic which means that Σ depends only on k 0 .The relation above fixes the modulus of the wavevector k = k R (r 0 ,t 0 , ω), with k R (r 0 ,t 0 , ω) = k 2 0 + Re Σ (r 0 , k 0 ,t 0 , ω).We now introduce the group velocity as 1/v g (r 0 ,t 0 , ω) = dk R (r 0 ,t 0 , ω)/dω and one can see that We note that for non-resonant scatterers the phase velocity v p equals the group velocity v g which also equals the transport (or energy) velocity v E .We thus have where k R = n R ω/c, n R being the real part of the effective refractive index of the medium.This leads to At this stage, it is important to note that the field entering the definition of the specific intensity [see Eq. ( 21)] is the analytic field Ē, whereas the field entering the expression of our current field-field correlation F [see Eq. ( 28)] is the real field E. Since ω and ω ′ remain close to the incident frequency ±ω 0 thanks to a weak Doppler shift at each scattering event and thanks to the assumption |Ω| ≪ {|ω|, |ω ′ |}, the field-field correlation F has two peaks in the frequency domains (one for positive frequencies and the other for negative ones) that do not overlap.Thus we can easily keep only positive frequencies which means that the real field can be replaced by the analytic field Ē in the F function.Using the definition of the specific intensity given by Eq. ( 21) and taking the inverse Fourier transform of Eq. ( 38) with respect to q and Ω, we finally get where du ′ means integration over the solid angle.In the above equation, r 0 and t 0 are slow variables compared to r and t, as a consequence of the quasi-homogeneous approximation.The equation is not modified if r 0 and t 0 are replaced by r and t, respectively.The expression of the RTE in the quasihomogeneous approximation and in the presence of inelastic scattering finally takes the form In this equation the extinction mean-free path and the scattering mean-free path are respectively given by 1 ℓ e (r,t, ω) while the phase function is given by We note that with these definitions the phase function is normalized as p(r, u, u ′ ,t, ω, ω ′ )du ′ dω ′ /(2π) = 1.
Equation (40), as the standard RTE, can be understood as an energy balance.The first two terms in the left-hand side correspond to the spatio-temporal evolution of the specific intensity, which itself can be seen as a local (in space and time) and directional radiative flux at frequency ω.The third term describes losses by extinction (absorption and scattering).The right hand side corresponds to gain by scattering.The phase function describes the fraction of the incoming power along direction u ′ and at frequency ω ′ that is scattered along direction u at frequency ω, with the change of frequency corresponding to the Doppler shift acquired during the inelastic scattering process.

C. Practical expressions of mean-free paths and phase function
In order to derive explicit expressions for the mean-free paths and the phase function, we need to start from the expression for the t-matrix (or scattering matrix) of an individual scatterer.In Fourier space, and for a particle with velocity v, it can be written as [30] with t(a, k, k ′ , ω) the t-matrix of a static particle with radius a.
The Dirac delta function accounts for the Doppler shift.Given the velocities involved, the condition |k • v| ≪ |ω| is satisfied, which implies that the shift can be neglected in the t-matrix, the latter being a slowly varying function of frequency.
The self-energy is a rather complex object and we need some approximations to go further.In the weak-extinction regime, we can use the lowest order approximation to Σ, know as the independent scattering approximation (ISA).This leads to where P(r i ,t, a i , v i ) is the probability density at time t of finding particle i at position r i , with a radius a i and a velocity v i .Note that this expression is very general, but in order to have the particle number density ρ(r i ,t) appearing explicitly we performed the following splitting.We assume that P(r i ,t, a i , v i ) = ρ(r i ,t)g(r i ,t, a i , v i )/N(t), where g(r i ,t, a i , v i ) is now the probability density at time t and position r i of finding particle i with a radius a i and a velocity v i .From this expression, using again the quasi-homogeneous approximation the self-energy in Eq. ( 45) can be rewritten as Taking the Fourier transform of Eq. ( 46) and making use of Eq. ( 44), we finally get Σ(r 0 , k R u,t 0 , ω) = ρ(r 0 ,t 0 ) t(a, k R u, k R u, ω)h(r 0 ,t 0 , a)da.

(47)
In this expression h(r 0 ,t 0 , a) = g(r 0 ,t 0 , a, v)dv is the probability density of having a particle with radius a, knowing that it lies at position r 0 at time t 0 .We note that the Doppler shift does not enter the expression of the self-energy.This can be easily understood: The self-energy defines the average field which propagates in the forward direction, as can be seen in the (k, k) dependence of the t-matrix in Eq. ( 47), for which there is no Doppler shift.Finally, the expression of the extinction mean-free path is 1 ℓ e (r,t, ω) = ρ(r,t) σ e (a, ω)h(r,t, a)da (48 where σ e (a, ω) = Imt(a, k R u, k R u, ω)/k R is the extinction cross-section of a particle with radius a.
To compute the phase function, we need to express the intensity vertex.The derivation is identical to that performed for the self-energy.In the same weak-scattering limit, the intensity vertex Γ is given by In particular, from Eqs. ( 49) and ( 44), we obtain In this expression, the Doppler shift appears explicitly.We also note that the square modulus of the t-matrix is proportional to the differential scattering cross-section of a particle through the relation [23] Inserting the above equation in Eq. ( 43) allows us to obtain the explicit expression of the phase function Since we consider spherical particles for the sake of simplicity, the differential scattering cross-section and the phase function only depends on u • u ′ , or equivalently on the angle between u ′ and u.

V. NUMERICAL SIMULATIONS
The quasi-homogeneous and inelastic RTE given by Eq. ( 40) cannot be solved analytically in real geometries such as that corresponding to a shock ejecta.Monte-Carlo algorithms have become powerful and convenient tools to solve the RTE numerically, especially thanks to their simplicity of implementation and to the ever-increasing available computer resources [38].In this section we provide a justification of the correspondence between the output of a Monte-Carlo simulation and the solution to the generalized RTE.

A. Integral form of the RTE
We start by writing the solution to the RTE in the form of an integral equation, suitable to be solved by a Monte-Carlo algorithm.To proceed, we note that the time derivative in Eq. ( 40) can be neglected since the illumination in the experiment is monochromatic (steady-state) and the transit time for light withing the scattering cloud remains short compared to the evolution time scale of the cloud T c .In this case, Eq. ( 40) simplifies into We see that the time variable t now plays the role of a parameter.In practice, this means that a steady-state Monte-Carlo simulation will be performed for each time t in the evolution of the spectrogram.Considering the right-hand side in Eq. ( 52) as a source term for the RTE, denoted by S(r, u,t, ω) in the following, we first define the Green function I 0 as the solution to with the boundary condition I 0 (r, r 0 , u,t, ω) → 0 when |r| → ∞.Defining r = r • u and r ⊥ = r − (r • u)u, it is easy to show that where the splitting δ (r − r 0 ) = δ (r ⊥ − r 0,⊥ )δ (r − r 0, ) has been used for the Dirac delta function, and H is the Heaviside step function.The solution to the RTE Eq. ( 52) can now be written which directly leads to This is the integral form of the quasi-homogeneous RTE, from which a Monte-Carlo scheme can be rigorously introduced.

B. Monte-Carlo scheme
The Monte-Carlo scheme can be seen as a random walk process for energy quanta (behaving as classical particles), where each step is characterized by statistical distributions in step length, scattering direction, and frequency.Equation (56) allows us to write explicitly the probability density l of having a step of length s, in the form where σ s (a, ω) being the scattering cross section of a spherical particle with radius a.We note that the probability density l is defined through the scattering mean-free path ℓ s only.Absorption is taken into account through a weighting factor (for each energy quantum in the random walk) given by where ℓ −1 a = ℓ −1 e − ℓ −1 s is the absorption mean-free path.Similarly, the probability density of having scattered direction and frequency (u, ω) for incident direction and frequency which is simply the phase function p introduced previously.
In practice, the optical properties of individual spherical particles, such as the extinction cross-section σ e and the differ- ential scattering cross-section dσ s /du, are computed numerically using the Mie theory [39], with an average over the two possible incident polarization states for unpolarized light.In the shock ejecta geometry, we also have to take into account the free surface (which corresponds to the remaining bulk portion of the metallic plate that was initially shocked).This is done by considering specular reflection at the interface to define a reflected direction, followed by the use of the Doppler shift formula to define the reflected frequency.

C. Implementation
In the practical implementation of the Monte-Carlo simulation, we discretize space into several layers or cells, depending on space variations of the statistical properties of the particle in the cloud.In a given layer or cell, the scattering and absorption mean-free paths ℓ s and ℓ a , and the phase function p, are taken to be constant.They are computed thanks to the routine given in Ref. 40.The probe geometry is mimicked by starting each random walk from the center of the probe, with a direction along the optical axis and towards the scattering medium.The collection by the probe is represented by an acceptance cone with angle θ p and the necessity for the energy quanta to hit a pupil with size φ p centered at the probe location.
To draw a distance s between two consecutive scattering events, we first normalize all lengths in each layer or cell by the corresponding scattering mean-free path ℓ s .Then, we randomly draw a distance d in an exponential distribution with unit parameter, and rescale the distances with ℓ s as we propagate from layer to layer or from cell to cell until we reach the end of the segment.To draw scattered direction and frequency, we first draw a velocity v and a particle radius a following the distribution g(r,t, a, v).Then, we draw a scattered direction u from the differential scattering cross-section, and finally we compute the outgoing frequency ω ′ using the Dirac delta function.This process is certainly not optimized in terms of computation time, but remains low cost in term of computer memory.Thus we are able to compute the detected flux G I(r, u,t, ω 0 + ω)u • ndudr which, once symmetrised on ω, gives the spectrogram at time t thanks to Eq. (20).Repeated for different times spaced on the order of T c , we reconstruct the full spectrogram S(t, ω).

D. Signature of multiple scattering
In order to analyze the influence of multiple scattering on the spectrograms by Monte-Carlo simulations, we need to specify the parameters characterizing a realistic cloud of particles.We consider the simple scenario used by Shi et al [20], and formulate the following assumptions: (1) the ejection is along the z direction, with particles propagating towards z > 0, (2) the cloud is translationally invariant in the x and y direction, (3) the particle positions z are directly given by vt, and (4) the particle radii and velocities are statistically independent, and do not depend on position and time, meaning that the probability density g takes the form g(r,t, a, v) = h(a) j(v).
Regarding the particle sizes, we assume a lognormal distribution [41] where µ = ln(µ 2 a / µ 2 a + σ 2 a ) and σ 2 = ln(1 + σ 2 a /µ 2 a ), µ a and σ a being, respectively, the mean and standard deviation of the particle radius a.For the velocity statistics, we assume a probability density where the parameter β gives the slope of the distribution, and v s is the velocity of the free surface [42].From these two distributions, we deduce that the number density of particles and the optical thickness take the form where V = (4π/3) a 3 h(a)da and σs = σ s (a)h(a)da are the average volume and average scattering cross section of the particles, respectively, ρ Sn is the volumetric mass den- sity and M s is the total ejected mass per unit area.We note that the expression (64) of the optical thickness is similar to that appearing in earlier works [13,[15][16][17]19], except that it is defined using the scattering mean-free path instead of the transport mean free path.We use the parameters [20] µ a = 0.75 µm, σ = 0.5, M s = 20 mg/cm 2 and β = 10, and we impose minimum and maximum values for the radii and the velocities a min = 0.1 µm, a max = 2.0 µm, v min = v s = 2250 m/s and v max = 4500 m/s.With these parameters the optical thickness is b = 42, showing that we consider the deep multiple scattering regime.In Fig. 2, we present numerical simulations of spectrograms for a given time in the dynamics of the ejecta.Figure 2 (a) compares a spectrum resulting from the full Monte-Carlo simulation (red solid line), and from a simulation restricted to single scattering events (black dashed line).The lineshape of the spectra are consistent with previously reported results [13,17].For high frequencies, the single scattering spectrum matches the spectrum computed with all multiple scattering sequences.This confirms that the front side of the ejecta scatterers chiefly in the single scattering regime, since the front side contains the fastest particles producing large Doppler shifts.For lower frequencies, as suggested by Franzkowiak et al. [14] and Andriyash et al. [13], the signal is dominated by high-order scattering.More specifically, for frequencies below the free-surface frequency ω s , the signal is only due to high-order scattering.This means that light that penetrates deep into the ejecta is inevitably multiply scattered.We also observe that while multiple scattering contributions could have been expected to generate random frequencies, they actually contribute in the same range as the single scattering spectrum.On closer inspection, we find that the multiple scattering sequences consist mainly of a single backscattering event surrounded by forward scattering events.Since forward scattering produces low Doppler shifts, according to Eq. ( 9), the cumulated Doppler shifts of forward scattering events remains small compared to the shift produced by backscattering. A a result, most multiple scattering sequences still encode the velocity of the particle they backscattered on, with an uncertainty produced by a broadening generated by the series of forward scattering events.The broadening remains small a the front side of the ejecta (where single scattering dominates) and becomes larger as light penetrates deeper into the particle cloud.
Figure 2 (b) compares spectra calculated for different optical thicknesses b of the particle cloud.To change the optical thickness, we take the surface mass M s as a tunable parameter, keeping all other parameters constant.We see that a peak corresponding to the velocity of the free surface emerges for b ≤ 4, and is not visible at larger optical thicknesses.This observation is consistent with results reported in previous works [15,17,20].

E. Inhomogeneities in the particle size distributions
The numerical simulations presented above included the full set of assumptions listed at the beginning of Sec.V D, including homogeneity in the particle size-velocity distribution.Molecular dynamics simulations [7][8][9] and holography measurements [41,43,44] suggest that fragmentation tends to produce smaller particles at the beginning of the ejection process and larger particles at the end.This results in a spatially inhomogeneous distribution of particle sizes.While the model based on Eq. ( 40) and the associated Monte-Carlo simulations handle inhomogeneities in the particle density, it might be relevant to include inhomogeneities in the statistical distribution of particle size.The only comparable works are from Andriyash et al. [16] and Kondratev et al. [17].Their model ap-plies to an ejacta expanding in air, in which the slowing down of particles depends on their size.This also results in an ejecta with an inhomogeneous size distribution, but of different nature.
In this work, we focus on the configuration studied by Sorenson et al. [43].An ejecta is produced in vacuum in a shock experiment, and a holographic setup allows them to measure the size distributions for five velocities in the ejecta.The experimental results confirm that the particles are larger at slower velocities.To include this degree of freedom in our model, we reformulate the basic assumptions as follows: the ejection is along the z direction, with particles propagating towards z > 0, (2) the cloud is translationally invariant in the x and y direction, (3) the particle positions z are directly given by vt, (4) the particle radii and velocities are statistically dependent, and do not depend on position and time, meaning that the probability density g takes the form g(r,t, a, v) = j(v)h(a, v).
Regarding the statistics of particle radii, we keep a lognormal distribution [41] with parameter µ(v) and σ (v) that now depend on the particle velocity.We fit these parameters on the five experimental distributions given in Ref. 43.Then we extrapolate this set of five couples (µ i , σ i ) by a fit of µ(v) and σ (v) on all velocities.This leads to µ(v) = ln 3.6185 From there, to retrieve j(v), we use the measured massvelocity distribution M(v) defined as where V (v) = (4/3)πa 3 ρ Sn h(a, v)da is the average volume of particles at velocity v, and V = V (v) j(v)dv is the average volume of particles in the entire ejecta.Experimentally M(v) follows an exponential law [42] such that M(v) = M s exp [−β (v/v s − 1)] with M s = 35 mg/cm 2 , β = 4.789 and v s = 2250 m/s.This gives This expression is similar to that given in Sec.V D, except that here the factor V / V (v) accounts for the inhomogeneity in the size distribution.As before, we can deduce the number density of particles and the optical thickness, that take the form where σs (v) = σ s (a)h(a, v)da is the average scattering cross section of the particles at velocity v.
Finally, to define a reference case for comparisons, we fit a lognormal size distribution over all particles in the ejecta.This describes an ejecta similar to that considered in Sec.V D, where inhomogeneities in the particle size distribution are not accounted for.In this case we get exp(µ) = 2.5750 × 10 −6 m and σ = 0.904, while keeping the same values for M s , β and v s .
Numerical simulations using this extended model are presented in Fig. 3.In Fig. 3 (a), we plot the size distributions at the front of the ejecta, at the back, and for a spatial average over the ejecta.We see that the inhomogeneous distributions are quite narrow compared to the distribution averaged over velocities, specifically at the front of the ejecta, and that they do not overlap much.This supports the idea that accounting for spatial inhomogeneity might be important.Indeed, the spectra shown in Fig. 3 (b) exhibit substantial differences.First, the real optical thickness in the inhomogeneous case is almost twice that obtained assuming a homogeneous distribution.This results from a rearrangement of matter in the ejecta that creates a dense and highly scattering region of small particles at the front, which screens the back part of the ejecta.While this could have been only the effect of a higher optical thickness, the homogeneous upscaled configuration in Fig. 3 (b) shows the opposite.Even at equal optical thicknesses, the visibility of the fastest particle at the front remains higher in the inhomogeneous case.This can be understood with the examination of Fig. 3 (c), where we plot the scattering mean-free path as a function of velocity.While shifting down the mean free paths using M s as a tuning parameter in order to reach the same optical thicknesses, the difference in slope of the ℓ s curve keeps a smaller mean free path at the front of the ejecta in the inhomogeneous case.In other words, the dense layer of small particles at the front of the ejecta generates a higher optical thickness.This effect is overlooked in models assuming homogeneous particle size distributions, introducing a potentially substantial bias in the interpretation of the spectrograms.

VI. CONCLUSION
In summary, we have developed a model for PDV measurements in shock experiments, emphasizing the influence of multiple scattering of the probe light.Starting from first principles of wave scattering, we have established rigorously the relationship between the specific intensity and the measured signal in PDV experiments, as well as the RTE that describes the evolution of the specific intensity upon scattering and absorption in the ejecta.The resulting RTE accounts for inelastic scattering (Doppler shifts) and inhomogeneities in the particle density, size and velocity in the quasi-homogeneous regime.The derivation provides a rigorous basis to the RTE model for PDV experiments, going beyond the usual phenomenological description, thus clarifying the underlying assumptions and limits of validity.We have also justified the use of a Monte-Carlo approach for numerical simulations, connecting the random-walk picture to the solution to the integral form of the RTE.Finally, we have used numerical simulations of realistic shock ejecta experiments to highlight the role of multiple scattering, and of inhomogeneities in the particle density and size-velocity distribution induced by the fragmentation process.We have shown that linking inhomogeneous source terms to PDV spectra is possible, and should be applied to other data sets than holographic measurements.In addition, since real ejecta have large optical thicknesses (b > 10), a diffusion approximation of the RTE could be derived in order to simplify the analyses in practice.These are potential lines to be followed in further investigations.

FIG. 1 .
FIG. 1. (a)Illustration of the micro-jet mechanism in a typical shock ejecta experiment.Upon reaching the machined free surface, the shock wave first comes into contact with the inwardly directed grooves.Under right angle conditions, the shock wave is reflected and the inward grooves become outward micro-jets.Due to the velocity gap between the jet-heads and the free surface, the micro-jets are stretched until surface tension is no longer sufficient to hold matter together and fragmentation begins.This results in the creation of an ejecta.(b) Schematic representation of a typical shock-loaded experiment with a PDV setup.The probe illuminates the ejecta and the free surface with a highly collimated laser beam (numerical aperture of 4.2 mrad and pupil size φ p = 1.3 µm).The backscattered field is collected by the probe as the measuring arm and interferes with the reference arm at the detector.The beating signal is registered with a high bandwidth oscilloscope before being analyzed.(c) Spectrogram of a tin micro-jetting experiment under pyrotechnic shock at P = 25 GPa.The free surface was engraved with 25 µm × 8 µm grooves.Independent Asay window measures gave an estimated surface mass M s = 5 mg/cm 2 .

FIG. 2 .
FIG. 2. (a) Comparison between single scattering (black dashed line) and multiple scattering spectra (red solid line).The simulations are carried out with the ejecta parameters defined in Sec.V D, and for a time t = 10 µs.(b) Spectra for different optical thicknesses.Same parameters as in (a).The variations in the optical thickness is obtained by changing M s , the total ejected mass per unit area.From left to right M s = 0.47 mg/cm 2 , M s = 1.41 mg/cm 2 , M s = 4.71 mg/cm 2 and M s = 20 mg/cm 2 resulting in b = 1, b = 4, b = 10 and b = 42.

1 FIG. 3 .
FIG. 3. (a)Comparison between the particle size distributions at the front of the ejecta (dotted line), at the back (dashed-dotted line), and for a spatial average over the ejecta (solid line).The parameters are given in Sec.V E. (b) Comparison between ejecta layouts.Same parameters as in Fig.3 (a), except for the upscaled homogeneous spectrum obtained by increasing M s in order to get an optical thickness b = 19.6.(c) Comparison between the scattering mean free paths for the three ejecta corresponding to the spectra in Fig.3 (b).