Self-Consistent Quantum-Field Theory for the Characterization of Complex Random Media by Short Laser Pulses

We present a quantum-field theoretical method for the characterization of disordered complex media with short laser pulses. We solve this scheme of coherent transport in space and time with weighted essentially non-oscillatory methods (WENO), that are suitable for the determination of highly nonlinear and discontinuous procedures including interference effects and Anderson localization of light. The theory determines spatio-temporal characteristics of the scattering mean free path and the transmission cross section that are directly measurable in time of flight and pump-probe experiments. The results are a measure of the coherence of multiple scattering photons in passive as well as in optically soft random media. They are instructive in spectral regions where material characteristics such as the scattering mean free path and the diffusion coefficient are methodologically almost insensitive to gain or absorption and to higher order non-linear effects. Our method is applicable to optical coherence tomography (OCT) and advanced spectroscopy setups including samples of strongly scattering monoand polydisperse complex nanoand micro-resonators.


INTRODUCTION
The characterization of disordered media has been fascinating the community ever since, and groundbreaking analysis methods like coherent backscattering (CBS) [1], dynamic light scattering (DLS) [2][3][4][5][6][7], diffusing wave spectroscopy (DWS) [8] and optical coherence tomography (OCT) [9][10][11][12][13][14] have been developed on the basis of transport of classical electromagnetic waves in random media and photonic crystals [1,15,16]. Classical methods have been generalized with great success for polarizationsensitive optical coherence tomography [17][18][19][20], for applications in metrology [21][22][23][24] as well as for opto-medical imaging [25][26][27]. Sub-and hyperdiffusive random media [28] have attracted great interest. Quantum-optical coherence tomography (QOCT) using entangled-photonsources in a Hong-Ou-Mandel interferometer [29,30] has been demonstrated [31][32][33]. It yielded an improvement of the resolution of a factor of two compared to OCT. Probing the submicron scale characteristics of transport of light is a crucial aspect in the understanding of dynamic properties of disordered random media [34][35][36]. Many of these approaches however do not account for multiple scattering at all, thus self-interference effects which can yield up to a factor of two enhancement of coherently scattered light with respect to the incoming intensity in passive scatterers ensembles [1], the CBS peak, are neglected. All these methods have in common that a systematic incorporation of multiple scatter-ing processes of light by optically soft scatterers, so the decoherence of light due to light-matter interaction, absorption, all orders of nonlinear processes [37] as well as scattering losses, are not incorporated in a systematic way [7,[38][39][40]. The range from weakly to strongly scattering non-conserving media so far is not covered by a systematic methodology [41][42][43][44] while quantitative fluorescence spectroscopy (QFS) in the time and the frequency regime is broadly investigated in turbid media and soft matter [45]. Technological applications in solid state physics and soft matter such as novel light sources, random lasers and solar cells based on multiple scattering in disordered arrangements of active resonators may profit from such novel techniques [46][47][48][49]. Non-invasive and non-destructive methods of optical analysis and medical imaging that can detect reflected signals as small as 10 −10 of the incident intensity and beyond [9,[11][12][13][14] might be improved again by orders of magnitude. Their application range could be systematically improved in the fields of dynamic and non-conserving media, nematic liquid crystals, semiconductors for telecom applications, glasses and tissue [34,[50][51][52][53]. In this article we develop a quantum field theory for photonic transport in dense multiple scattering complex random media, see Fig. 1(a). The Bethe-Salpeter equation, which governs the propagation of the intensity [54,55], is solved for a propagating short laser pulse in random media with the help of a weighted essentially non-oscillatory solver (WENO) [56][57][58][59][60] in the space and time dependent framework [61], Fig. 1(c). We are going beyond the diffusion approximation by including interferences and repeated self-interferences in the sense of all orders of maximally crossed diagrams, the Cooperon [62,63]. This is generally associated with the Anderson transition of electrons as well as of light, sound and matter waves in multiple scattering random media [64][65][66][67][68][69][70][71][72][73][74][75], and thus with quantum effects [62,63]. The random medium is assumed to be optically complex and we are including nonlinearities, absorption and gain, which in consequence request the implementation of suitable conservation laws by means of the Ward-Takahashi identity for non-conserving media and resonators. The random medium can be comprised of ensembles of arbitrarily shaped particles as well as correlated disorder and glassy systems in principle [76]. We focus in this work on independent non-conserving Mie scatterers [77,78] in strongly scattering ensembles of a high filling fraction, for instance mono-and polydisperse complex TiO 2 powders [79,80]. Such ensembles are well known for showing pronounced Mie signature in their transport characteristics such as the scattering mean free path as it has been determined experimentally also by coherent backscattering for optically passive systems [1,16,41], Fig. 1(d). We derive in what follows self-consistent results for conserving and for non-conserving random media and we show that absorption and gain or non-linearities can be characterized in a time of flight experiment by the fraction of transmitted photons which experience a delay due to coherent multiple scattering. It will be shown that the deviation in the long time limit from the pure diffusive case provides fundamental knowledge about the subtle nature of the scattering ensemble and it's complexity in the sense of the resonator properties of the single scatterer.

Nonlinear Response
The electrodynamics for transport of light in random media is described basically by the wave equation where the polarizability in bulk matter P ( E) can be decomposed in linear and non-linear part = P L + P N L .
As opposed to the non-linear Schrödinger equation for matter waves [81] the wave equation allows in the presented quantum-field theoretical formalism for the straight forward incorporation of the Mie resonance [82] as a classical geometrical effect in the sense of a whispering gallery resonance of the light wave at the inner surface of the complex scatterer [83]. In general the polarizability is defined as P ( E) = 0 ( 1 − 1) E + P N L . 0 is the dielectric constant in free space, is defined in the literature as the material specific dielectric coefficient = 1 + χ, c is the speed of light. Higher order processes, for instance Kerr media [37] with the dielectric susceptibility χ (2) [79,80], are in electrodynamics classified by the dependency of P to the electrical field E [84] without loss of generality. It is well known that both the conductivity as well as the susceptibility contribute to the permittivity, so in general Im = Im(χ) + Re(σ/ω) is given. Absorption and optical gain are represented by a finite positive or negative imaginary part of the dielectric function, so in general Im = 0 is assumed. We take into account the Mie scatterer Fig.(1) for the determination of the single particle self-energy contribution Σ ω k of the quantum fieldtheoretical approach in what follows. The Mie scattering coefficients of n-th order are written as [82,85] where m here denotes the complex refractive index, y = 2π r scat /λ is the size parameter, λ is the wavelength of light and r scat is the spheres radius. Prime denotes the derivative with respect to the argument of the function. Ψ n and ξ n are Riccati-Bessel functions defined in terms of the spherical Bessel function and in terms of the Hankel function [82,85]. The characteristics of the active scatterer and the active embedding matrix are described by Im scat and Im b in the following.

Selfconsistent Bethe-Salpeter Equation for Inelastic Multiple Scattering of Photons
For the theoretical description we may use any distribution and shape of particles which may be described in the form of a scattering matrix. Here we consider in our results monodisperse spheres as well as a Gaussian distribution of spherical scatterers located at random positions [72][73][74][86][87][88][89][90][91][92][93], see Fig. (1). The scatterers as well as the background medium are described by the dielectric constants scat and b . In this work we use unpolarized light and therefore we consider the scalar wave equation which has been Fourier transformed from time t to frequency ω and reads where c denotes the vacuum speed of light and j ω ( r ) the current. The current j ω ( r ) may be expanded in orders of ω [94,95]. We do not take into account here a coupling to a microscopic model for dynamical feedback of the optically driven crystal in the non-equilibrium [96,97] or chaotic feedback [98]. The dielectric constant is spatially dependent, ( r ) = b + ∆ V ( r ), and the dielectric contrast is defined as ∆ = scat − b . The dielectric contrast describes in principal the arrangement of scatterers through the function V ( r ) = R S R ( r − R ), with S R ( r ) a localized shape function at random locations R.
The intensity is related to the field-field-correlation function Ψ( r, t )Ψ * ( r , t ) , where angular brackets . . . denote the ensemble or disorder average. To calculate the field-field-correlation the Green's function formalism is used, the (single-particle) Green's function is related to the (scalar) electrical field by The Fourier transform of the retarded, disorder averaged single-particle Green's function of Eq. (4) reads, where the retarded self-energy Σ ω q arises from scattering of the random potential −(ω/c) 2 ( scat − b )V ( r ). Using Green's functions the mode density, the local density of photonic states (LDOS), N (ω) may be expressed as N (ω) = −(ω/π)ImG ω 0 , where we use the abbreviation We study the transport of the already introduced field-field-correlation by considering the 4-point correlator, defined in terms of the non-averaged Green's functionsĜ,Ĝ * in momentum and frequency space as Φ ω q q ( Q, Ω) = Ĝ ω+ Here we introduced [83] the center-of-mass and the relative frequencies (Ω, ω) and momenta ( Q, q) with ω 1,2 = ω ± Ω/2 and q 1,2 = q ± Q/2. The variables Ω, Q are associated with the time and the position dependence of the averaged energy density, withQ = Q/| Q|, while ω ± = ω ± Ω/2 and q ± = q ± Q/2 etc. are the frequencies and momenta of in-and out-going waves, respectively. The intensity correlation, the disorder averaged particle-hole Green's function, Φ ω q q ( Q, Ω) is described by the Bethe-Salpeter equation By utilizing the averaged single particle Green's function, c.f. Eq. (6), on the left-hand side of Eq. (7) the Bethe-Salpeter equation may be rewritten as the kinetic equation, see, e.g., ref. [83], When we analyze the correlation function's long-time (Ω → 0) and long-distance (| Q| → 0) behavior, terms of O(Ω 2 , Q 3 , ΩQ) are neglected. Eq. (8) contains the total quadratic momentum relaxation rate 1/τ 2 = c 2 Im( b ω 2 /c 2 − Σ ω ) due to absorption/gain and due to impurity scattering in the background medium as well as the irreducible two-particle vertex function γ ω q q ( Q, Ω). In order to solve this equation we expand it into moments.
The energy conservation is implemented in the solution of the Bethe-Salpeter equation by a Ward identity (WI) for the photonic case, see Ref. [83]. The Ward identity is derived in the generalized form for the scattering of photons in non-conserving media. Non-linear effects, absorption and gain yield an additional contribution, and a form of the Ward-Takahashi identity for photons in complex matter [83,99,100] is derived. The additional contribution is not negligible and thus effectively present in all results of the transport characteristics of the selfconsistent framework [54,55,83]. For scalar waves the Ward identity assumes the following exact form The right-hand side of Eq. (9) represents reactive effects (real parts), originating from the explicit ω 2 -dependence of the photonic random potential. In conserving media (Im b = Im scat = 0) these terms renormalize the energy transport velocity v E relative to the average phase velocity c p without emphasizing the diffusive long-time behavior. [83,101] In presence of loss or gain, additional effects are enhanced by the prefactor which now does not vanish in the limit of Ω → 0.

Expansion of the Two-Particle Green's Function into Moments
The solution of the Bethe-Salpeter equation is derived by rewriting it in the form of a kinetic equation and by deriving a continuity equation. For this aim we expand the intensity correlator into its moments and we derive a diffusion pole structure out of the Bethe-Salpeter equation Eq. (7). The q integrated correlator is decoupled from the momentum dependent prefactors by some auxiliary approximation scheme. This approximation must obey the results of the ladder approximation and it must incorporate the set of physical relevant variables for the observed phenomena. In a first step we use the bare first two moments of the correlation function Φ q defined as These bare moments are related to physical quantities, the energy density correlation P ω E ( Q, Ω) and the currentdensity-correlation J ω E ( Q, Ω), by dimensional prefactors: J The projection of the correlator Φ q , Eq. (11), onto the bare moments Φ ρρ ( q, Ω) as in Eq. (12), and Φ jρ ( q, Ω), as in Eq. (13), is therefore then by The projection coefficients A( q ) and B( q ) are to be determined in the following. In this expansion the bare moments may be substituted by their physical counterparts, the energy density P ω E in Eq. (14) and the current density J ω E ( Q, Ω) from Eq. (15). The expansion coefficients A( q ) and B( q ) in Eq. (16) behave uncritically when the system localizes. Thus they can be determined by using the simple ladder approximation, where all expressions are known exactly. The ladder approximation of the two-particle vertex function is schematically explained in [54,55,83].
In the following we use the approximation to obtain the expansion coefficients from it. In the ladder approximation the zeroth bare moment is given by the superscript L refers to the ladder approximation.
The product has been expanded up to linear order in q. The renormalized vertex γ 0 is given bỹ where γ 0 is the bare vertex. f ω (Ω) is arising from the Ward identity and has been defined in Eq. (9). Within the simple ladder approximation the bare moment Φ L jρ ( Q, Ω) defined in Eq. (13) is thus given by We follow this strategy again and by expanding the prod- under the second integral up to first order in q we obtain the expression By employing the same expansion to the remaining product of the Green's function one eventually finds where the abbreviation ∆G ≡ G − G * has been introduced.
In the next step of determining the expansion coefficients A( q ) and B( q ), Eq. (16), we go back to the fieldfield correlation function Φ q q . In the uncritical ladder approximation the two particle Green's function is given by Employing the momentum expansion again, the equation, Eq. (22) can be simplified to yield By using the above given momentum expansion, Eq. (23), in combination with the expressions given in Eq. (21) and in Eq. (17) in connection with the described projection, or expansion into moments, Eq. (16), the following relation is eventually derived By comparison of coefficients in the relation, Eq. (24), the demanded coefficients A( q ) and B( q ) of the expansion into moments, Eq. (16), can now be determined as follows Employing those expressions for the expansion coefficients, we can eventually express the two-particle correlator Φ q q as follows The expression, Eq. (26), represents the complete expansion of the intensity correlator into its moments. This form is used now to decouple and therefore solve the Bethe-Salpeter equation.

General Solution of the Bethe-Salpeter Equation
We repeat the most important steps so far. The disorder averaged intensity correlation, the two-particle Green's function, obeys the Bethe-Salpeter equation, see Eq. (7) Φ The Bethe-Salpeter equation may be rewritten into the kinetic equation given in Eq. (8) To find the solution of Eq. (28), we first sum in Eq. (28) over momenta q, and we incorporate the generalized Ward identity as given in Eq. (9) and we expand the obtained result for small internal momenta Q and internal frequencies Ω. It is essential to employ the form of the two-particle correlator shown in Eq. (26). Eventually the generalized continuity equation for the energy density can be derived as The generalized continuity equation represents energy conservation in the presence of optical gain and/or absorption.
As the standard solution procedure the next step is to obtain a linearly independent equation which relates the energy density P ω E and the current density J ω E . This is realized in a similar way to above, one first multiplies the kinetic equation, Eq.(28), by the projector q ·Q and then follows the already outlined steps to eventually obtain the wanted second relation. This is the current relaxation equation relating the energy density P ω E and the energy density current J ω E . The memory function M (Ω) is introduced according to where γ ω p p ≡ γ ω p p ( Q, Ω) is the total irreducible twoparticle vertex, which will be discussed in detail in what follows.
So far, two independent equations, Eq. (29) and Eq. (30), have been obtained. Either of them is relating the current density J ω E with the energy density P ω E . Now one can eliminate one of the two variables in this linear system of equations. The two equation are combined to find an expression for the energy density that exhibits the expected diffusion pole structure for non-conserving media. Precisely in the denominator of Eq. (32) there appears an additional term as compared to the case of conserving media. This is the term ξ −2 a D, the mass term, accounting for loss (or gain) to the intensity not being due to diffusive relaxation. In Eq. (32) also the generalized, Ω-dependent diffusion coefficient D(Ω) has been introduced by the relation It shall be noted that Eq. (32) introduces the absorption or gain induced growth or absorption scale ξ a of the diffusive modes, which is to be well distinguished from the single-particle or amplitude absorption or amplification length. The diffusion constant without memory effects in Eq. (33), consists of the bare diffusion constant [101], and renormalizations from absorption or gain in the background medium (D b ) and in the scatterers (D scat ), The upper panel shows a diagrammatic expansion of the irreducible two-particle vertex γ. The lower panel displays the disentangled Cooperon with changed momentum arguments.
whereD 0 is the same as in Eq. (35), with (ImG ω q ) 2 replaced by Re(G ω 2 q ). In the above Eqs. (34)-(36) the following short-hand notations have been introduced, Vertex Function and Self-Consistency According to Eq. (31) and Eq. (33) the energy density, the two-particle function, given in Eq. (32) still depends on the full two-particle vertex γ ω q q . Before discussing the vertex function, we briefly recall our arguments with regards to dissipation. Dissipation breaks the time reversal symmetry [102,103] on the one hand side but on the other hand the dissipation rate itself is invariant under time reversal. As a general picture of this physics the damped harmonic oscillator can be mentioned, where the time reversed solution is still damped with the very same damping constant. Having this in mind we analyze the irreducible vertex γ ω q q for the self-consistent calculation of M (Ω), exploiting time reversal symmetry of propagation in the active medium. In the long-time limit (Ω → 0) the dominant contributions to γ ω q q are maximally crossed diagrams (Cooperons), which are valid as well as for conserving media, and they may be also disentangled.
In Fig. (2) the disentangling of the Cooperon into the ladder diagram is shown. The internal momentum argument of the disentangled irreducible vertex function in the second line of Fig. (2) is replaced by the new momentum Q = k + k . Thus γ ω q q acquires now the absorption or gain-induced decay or growth rate ξ −2 a D.
Finally the memory kernel M (Ω) reads Eqs. (33)-(37) constitute the self-consistency equations for the diffusion coefficient D(Ω) including the growth/decay length scale ξ a in presence of dissipation or gain.

Length and Time Scales
Within disordered systems a multitude of length an d time scales are defined, that are related to the single or the two-particle quantities respectively. An important length scale which can be directly measured in the experiment is the scattering mean free path l s defined in the single particle Green's function where the imaginary part of the self-energy introduces the decay length l s The decay length may equivalently be interpreted as the life time of the corresponding k-mode. In the case where the dielectric constant is purely real, so in the case of passive matter, the scattering mean free path l s describes the scale for determining the loss sole to scattering out of a given k-mode. In the other case, for gain and dissipation, the k-mode experiences amplification or absorption. In the case of gain this transport theory is valid for ImΣ(ω) < 0 while the flip of ImΣ(ω), ImΣ(ω) = 0, defines the point of the phase transition, i.e. the laser threshold, for the pumped single scatterer [54,55,[104][105][106][107].
We discuss in what follows the transport of the intensity and the scales related to it. The two-particle Green's function as given in Eq. (32 ) contains two obvious scales originating solely from finite values of the gain/absorption coefficient. These length scales may be defined by where a represents the amplification or absorption length of the intensity and osc marks the length over which the intensity oscillates, where ξ 2 a has already been defined in Eq. (34 ). The corresponding time scales may then be defined as Including the gain induced growth rate τ a as defined in Eq. (43), the intensity Green's function Eq. (32 ) may now be rewritten as where the coefficient α may symbolically contain all the factors explicitly shown and discussed in Eq. (32). Our aim in this article is to calculate the electrical fieldfield-correlator at different positions and frequencies Eq. (7) eventually leading to the evaluation the two particle Green's function given in Eq. (32). The momentum Q appearing in Eq. (32) represents in Fourier space a relative position within the sample. In three dimensions the momentum Q actually defines a volume unit within the sample. This volume is carefully to be distinguished from all other length scales e.g. the sample volume etc. It is merely the scale which determines the presence of correlation effects in photonic transport.
In analogy to the flip of the resonance, ImΣ = 0, in the single particle Green's function the equivalent threshold condition for the energy density is found as follows It leads to the critical length scale This length describes the volume where photonic transport, i.e. the energy density or intensity, may compensate diffusive losses by amplification due to the presence of some finite optical gain.

Weighted Essentially Non-Oscillatory Solver (WENO)
For the time and space dependent solution of the diffusion equation, Eq. (45), including a coherent laser pulse, see Fig.(1), we use a weighted essentially non-oscillatory method (WENO) in time in combination with a fourth order Runge-Kutta method in space.
Like the discontinuous Galerkin method [61] for hydrodynamic systems, the WENO method [56] has been specifically developed for discontinuous and rogue processes, shocks and steep gradients. Such procedures are well known to cause numerical problems or oscillations in the calculation of the first derivative. An efficient method is thus needed to refine the discretization of the problem locally in space and time. WENO is as such an upgrade of the essentially non-oscillatory method (ENO) [57,58] which has been developed for the calculation of hyperbolic conservation laws. ENO replaces the calculation of higher order difference quotients by the calculation of a bunch of lower order difference quotients which are of equal order. Whereas ENO incorporates only the difference quotient with the smallest approximation and thus always the influence of a part of the supporting points of a number of cells of the so-called stencil is neglected in the search for the stencil with the smoothest result for the interpolation, the WENO method is more sophisticated. It uses a convex combination of all candidates of lower order difference quotients of the stencil with an attributive weight g i [59] where D i are the smoothness indicators of the stencil. The variable h > 0 is defined as the machine accuracy which prohibits a division by 0. All weights are normalized to unity.g If the stencil contains a discontinuity the smoothness indicator should be essentially 0. The convergency of the stable solution is guaranteed by the Lax equivalence theorem [60].

Scattering Mean Free Path ls and Diffusion Constant D for Mono-and Polydisperse Passive and Active Scatterers
First we discuss here as results the scattering mean free path l s , see Eq. (40), and the diffusion constant D, see Eq. (33), as the self-consistent material characteristics of the disordered sample of complex Mie scatterers [108][109][110]. As the materials initial parameters we refer in what follows to the literature value of the passive refractive index for titania TiO 2 , n = 2.7.
In Fig. (3) we show the scattering mean free path l s for a Gaussion distribution of Mie scatterers, see Fig. (1), which is centered at the radius of r scat = 122.5 nm, the full width half maximum is 17.0 nm, and we compare it to results for monodisperse Mie scatterers of size r scat = 122.5 nm. The description takes advantage of the fact that scattering matrices of independent scatterers are additive in general. We find that the principal Mie characteristics of the central particle size r scat = 122.5 nm is qualitatively conserved but quantitatively reduced. This is intuitively clear due to the additivity of scattering matrices in the independent scatterers approach since no additional structural effects, e.g. in the sense of a varying concentration of surface defects or the occurrence of correlated clusters and glass transitions, are considered so far to induce any additional dependencies. The exact Mie resonance positions, the minima of l s , remain spectrally fixed for the Gaussian distribution of polydisperse scatterers ensembles compared to the monodisperse ensembles. What can be definitely deduced is that the scattering mean free path l s overall is prolonged for the Gaussian distribution which means that the scattering strength of the disordered sample is effectively reduced for polydisperse media. perse scatterers distributions compared to monodisperse ensembles. In Fig. (4) we show the results for the scattering mean free path l s for the Gaussian distribution of scatterers of the filling fractions 48.83 % and 55.10 % for passive Mie resonators as well as for absorbing scatterers. For the absorption we consider the literature value for TiO 2 bulk of Im(n) = 0.005932 as well as experimentally relevant absorption values for disordered granular arrangements. We find that the material characteristics of the scattering mean free path l s for an increase of the filling fraction from 48.83 % to 55.10 % is quantitatively reduced, whereas is is qualitatively confirmed. No crossover between both results is found all over the spectrum.  Fig.4. We find that the increase of the filling fraction rather moderately effects the quantitative behavior of ReD in the spectral region below λ < 490.0 nm, so λ < 4 · rscat, for λ > 490.0 nm we find, that the influence of the filling increases. We find that the effect of the bulk value for absorption of Im(n) = 0.005932 on the diffusion coefficient for moderate excitation is peaked at λ = 350.0 nm (green line), while it is almost undetectable in other spectral positions. The increase of absorption to τa = 0.9 ns brings the same crossing point as found for ls at the inflection points. Its physical meaning is rather important. Absorption or gain will rather hardly be measured at the inflection points. In the inset we show the crossing at λ = 583.0 nm which has been confirmed in [72,74] experimentally. It is found that a value of absorption, as it is commonly detected in experiments with TiO2 powders [72,74] can eliminate almost completely the size dependent influences of the scatterers geometry in the diffusion characteristics of random ensembles.
istics at λ = 345.0 nm to λ = 360.0 nm, λ = 430.0 nm to λ = 445.0 nm show a decrease of l s for moderate absorption. By increasing the absorption up to the value of Im(n) = 0.3 which is corresponding to τ a = 0.9 ns we find that the Mie resonances are further reduced. For the filling fraction of 55.10 % there is a crossover of all results of l s found. It can be deduced that for one specific distribution of Mie scatterers sizes and varying loss parameters and identical parameters otherwise spectral points or at least narrow spectral regions exist where the loss is barely detectable by the coherent backscattering or the coherent forward scattering experiment. This lossinsensitive point is located approximately at the turning points of the results for l s . In the dips and the peaks of l s loss and gain is most efficiently detected. Speaking in terms of the scattering strength of the disordered sample, it is interesting, that while absorption reduces the resonators influence in general and thus leads to the reduction of the scattering strength in the Mie resonance . The result is normalized to the peak transition time. We display the passive complex system (black), the same arrangement including the literature value for absorption of TiO2 bulk, Im (n) = 0.005932 (green), an experimentally relevant value for absorption of disordered random TiO2 media Im (n) = +0.3 (red), and the same scatterers arrangement with refractive index n = 2.7 as it is pumped and gain assumes a value of Im (n) = −0.3 (yellow). For all parameters a plateauing in the temporal evolution of σ 2 (t)/L 2 is found as a characteristics of the complex medium in the long time limit. Whereas for gain the spreading of the cross-section overall is increased, absorption generally seems to lead to an earlier inset of the plateauing effect. Results for the passive system (black), Im (n) = 0.005932 (green) and Im (n) = +0.3 (red) are to be compared to results of the scattering mean free path ls Fig. (3), and to the diffusion coefficient D,  For an absorbing polydisperse sample, τ a = 0.9 ns and Im(n) = 0.3, we find that almost any Mie characteristics is washed out due to the absorptive character of the single scatterer. A destructive interplay between absorption characteristics and the resonator properties is developed all over the electromagnetic excitation spectrum. Extremely pronounced is the difference in the magnitude of the diffusion coefficient D for wavelengths larger than λ = 700.0 nm. When the filling is enhanced from 48.83 % to 55.10 % the diffusion coefficient D at λ = 800.0 nm is reduced from 85 m 2 /s to 25 m 2 /s so approximately to 30 %. Is shall be pointed out here that there is no transition to an ordered sample considered and all results are derived for homogeneous disordered but polydisperse samples of the mentioned Gaussian distribution of scatterer radii. We display one of the crossover points in the inset of Fig. (5) the diffusion coefficient D which is in the vicinity of the wavelength λ = 583.0 nm.
When we discuss these results for l s and D, in terms of the Ioffe-Regel criterium with the benchmark of strong localization of light, we find kl < 1 could be experimentally derived for the monodisperse case, see Fig. (3), e.g. in the Mie resonance at λ = 490.0 nm, whereas for the polydisperse ensemble kl = 1.67 the condition is not strictly fulfilled and interference effects will play a subtle role. It can thus be concluded that the probability to find Anderson localized photons will be enhanced in the monodisperse ensemble where a factor of 2 enhancement with respect to the incoming intensity could be detected. The increase of the volume filling fraction for polydisperse samples will only have a moderate impact in the search for Anderson localized photons, whereas absorption, which is enhanced in disordered granular media as compared to the bulk case, is a crucial and limiting factor, see Fig. (4). Thus it is important to find a systematic theoretical method to distinguish micro-and macroscopic structural effects in the signal which is a mix of both.

Temporal Evolution of the Transmission Cross
Section and the Transmitted Intensity In section Quantum-Field Theory for Multiple Scattering of Photons we presented the theory to study the propagation and the localization of a laser pulse through a disordered ensemble of complex scatterers. The solution of this framework in the sense of the simulation of detectable scattered intensity is not restricted to a dimension or a direction. Transmission optical coherence tomography based measurements of optical material properties is one experimental platform for our methodology [111][112][113]. Here we present results for the temporal evolution of the transmission cross section σ, and the mean square width σ 2 respectively, Fig. (6), and the temporal dependency of the transmitted photonic intensity, Fig. (7). Our results are derived for homogeneous disordered samples with a Gaussian distribution of scatterer diameters centered at d = 245.0 nm and 55.10 % volume filling. This is the strongly scattering and strongly disordered regime. We display results for the excitation laser frequency of λ = 590.0 nm, where we know kl s = 3.62 from previous results, see Fig. (4). In terms of the Ioffe-Regel criterium this case should be far off any case of the Anderson localization regime, thus the aim of our considerations is to determine the amount of coherently interfering photons as the deviation from the purely diffusive case. The mean square width [72,74,88] is defined as the square of the up to the full-widthhalf-maximum limit integrated area σ of the transmitted transverse intensity distribution at the ensemble surface, see Fig.(1). We display this characteristics in Fig. (6) normalized to the square of the samples length L 2 .
We find that the transmission cross section, the mean square width, Fig. (6), is definitely depending on absorption and gain of the complex random medium, and it is a very sensitive measure. In Fig. (6) we display results for the temporal evolution of σ 2 /L 2 for the wavelength of the incident pulse of λ = 590.0 nm. Th case of λ = 590.0 nm is extremely close to the crossing point or turning point of the scattering mean free path l s , Fig.  (4), and the corresponding diffusion coefficient D, Fig.  (5). While the characteristics l s and D are almost insensitive to absorption at λ = 590.0 nm, a clear deviation of σ 2 /L 2 for the case of Im(n) = +0.3 from the passive case of about 12 % is found in the temporal limit of t/τ max = 4. In Fig. (7) we show the results for the transmitted intensity I/I incident in the long time limit. The purely diffusive case is marked in the log-plot as the dashed green line. It can be concluded for the strictly passive case (black line) that the increase of the transmission in the long-time limit in comparison to the diffusive case with the exponential decay is the amount of localized photons, that have been multiply and coherently scattering and interfering in the disordered medium in the sense of the maximally crossed procedures, represented diagrammatically by the Cooperon, see Fig. (1). It is further derived that absorption, incorporated by the literature value for bulk TiO 2 , Im(n) = +0.005932 (blue line), reduces the amount of localized photons, and the long time behavior of P E retrogrades towards the diffusive limit. This result as such has been expected. By comparison of the relative magnitude of the results for σ 2 /L 2 with the transmitted intensity P E , Fig. (7), it is however interesting to note, that the plateauing effect for the case of absorption is enhanced, precisely it shows an earlier onset. Thus the plateau as such seems not to be a signature of localization, however the magnitude of the deviation of the transmitted intensity P E from the diffusive limit can be interpreted as a sign of enhanced coherent multiple scattering and thus as an enhancement of interference effects in principle. The influence of the single Mie scatterer is noteworthy when we discuss the influence of gain. Gain as the negative imaginary part of the complex refractive index and as a negative part of the complex permittivity is a microscopic material characteristics equivalent to absorption. Whereas absorption however is a microscopic interaction where the life-times of light-matter bound states is irrelevant, this is different for the case of gain. Gain as such is achieved by an enhanced life-time of light matter bound states leading to an increase of the photon number in the incident wavelength. Gain and absorption are as such in our theory properties of the single complex Mie scatterer, and the microscopic materials characteristics is interacting with the resonator. In the case of strong external laser pulses the local density of states and thus the refractive index of TiO 2 bulk can be shifted and this effect will lead to a shift of the overall characteristics of the disordered granular medium l s , Fig. (4), and D, Fig. (5). When nonlinear effects, e.g. higher order harmonics of the incident wavelength ω, play a role, this processes might contribute to an enhancement of multiple scattering and of interference effects. These enhancements can result in time of flight experiments is a variety of observations, such as offcentered peak in the integrated spectrum. In this article we discuss gain in the central wavelength of the incident pulse, which can originate from light-matter bound states of higher harmonics, so non-linear effects play a role even though the system is far below any laser threshold. The gain is visible on the one hand side as a delayed onset of the plateauing of σ 2 /L 2 on the one hand, see Fig. (6) (yellow line), on the other hand we find an increase in the transmitted coherent photon intensity in the long-time limit, Fig. (7). The variation of P E due to enhanced multiple scattering and interference effects due to an increased filling of due to enhanced resonator properties of the geometrical scatterer can thus be well determined and distinguished by comparing the characteristics of l S , D, σ 2 /L 2 and P E from effects incorporating gain and absorption.

CONCLUSIONS
We have presented in this article a novel method for characterizing disordered complex random media based on a quantum-field theoretical approach, the Vollhardt-Wölfle theory, for photonic transport including interference effects in multiple scattering processes. Our theory incorporates the Ward-Takahashi identity for photonic transport and thus enables us to determine self-consistent results for the material characteristics of disordered granular complex media. We solve the theory in three dimensions space-and time-dependent with a weighted essentially non-oscillatory solver method (WENO). The solution with the WENO solver enables us to determine results for ultrashort analyzing pulses and pump-probe experiments, since we can deal with highly non-linear processes and discontinuities. We have presented here a systematic study of the scattering mean free path l s and the diffusion coefficient D. These characteristics can be compared directly to experimental results derived in a coherent backscattering experiment. They show that the resonator characteristics, e.g. the Mie resonance, plays a crucial role which is even more important, when non-linear effects, gain and absorption, are minimized by the choice of the scattering matter. It has also been shown that a polydispersity of the scatterers reduces the probability to reach an Anderson transition with transport of light. The incorporation of gain and absorption reveals, that all materials characteristics are very sensitive to such properties of complex matter. It has been shown that so called absorption-free measurements are bound to spectrally narrow areas where the resonator characteristics of the single scatterer leads to a minimized sensitivity with respect to a change of the complex refractive index and the complex permittivity. Our results for ab initio simulations of time of flight experiments yield the characteristics of the normalized transmission or reflection cross section and the absolute as well as the normalized number of coherently scattered photons. We presented results random mono-and polydisperse ensembles of TiO 2 Mie scatterers in a transmission optical coherence tomography setup. It has been demonstrated that these characteristics offer an increased sensitivity to any microscopic and the macroscopic structural modification compared to the coherent backscattering experiment. The underlying theory paths a way towards the detection of subtle interference effects due to multiple scattering events in OCT setups that may lead to an increase of the sensitivity of OCT of orders of magnitude. This effect can be enhanced by occurring scatterers resonances. We conclude that our combinatory analysis of underlying transport theory and its results, the scattering mean free path the diffusion constant, and the derived characteristics in the temporal evolution is suitable to distinguish between perfectly coherent multiple scattering and interferences and on the other hand between influences of the complex random medium in the full spectrum of the analysis. We provide a consistent method, which is able to characterize disordered media in the weakly and the strongly scattering regime, the approach is suitable to incorporate ultrashort and intense light pulses and the resulting subtle local and non-local light-matter interactions on a broad temporal range. The method is ab initio not limited to light, it can be performed with the full spectrum of electromagnetic excitations and it can be transferred to any other type of wave propagation as matter waves and sound. It will be subject of subsequent work to investigate further influences of multiple scattering and higher order nonlinear effects in ensembles of clusters and composite scatterers such as shells, as well as of macromolecules that can be random or quasi-ordered and may form so called meta glasses.
Authors contributions: All authors were equally involved in the preparation of the manuscript. All the authors have read and approved the final manuscript.