Theory of high-energy correlated multiphoton x-ray diffraction for synchrotron radiation sources

We present a theoretical formulation for the multiphoton diffraction phenomenology in the nonrelativistic limit, suitable for interpreting high-energy x-ray diffraction measurements using synchrotron radiation sources. A hierarchy of approximations and the systematic analysis of limiting cases are presented. A convolutional representation of the diffraction signal allows classification of the physical resources contributing to the correlation signatures. The formulation is intended for developing a theoretical description capable of describing plausible absence or presence of correlation signatures in elastic and inelastic diffractive scattering. Interpreting these correlation signatures in terms of the incoming field modulated many-body electronic density correlations provides a novel perspective for structural imaging studies. More essentially, it offers a framework necessary for theoretical developments of associated reconstruction algorithms.


I. INTRODUCTION
In the x-ray photon-induced diffractive scattering measurements where signals are monitored in the far-field zone, one of the central goals lies in extracting information about the hitherto unknown structure of the materials. In these measurements, high-resolution structural information in the real space are routinely obtained by monitoring the photon scattering patterns and applying suitable inversion schemes [1][2][3][4]. Obtaining atomically resolved structural information, which has been one of the central endeavors of the x-ray crystallographic measurements, has seen enormous progress in recent years. Such improvements in precision and resolution were facilitated by the deployment of advanced light sources [5][6][7][8][9][10][11][12][13]. For structural studies, hard x-ray sources whose wavelength corresponds to the atomic resolution are preferred as a tool. The fact that the probe wavelength is comparable to the typical length scales related to relative distances between the scattering centers (e.g., 1Å corresponds to 12.4 keV) allows imaging at the atomic resolution. The operational regime of the hard x-ray sources corresponds to the photon energies ranging from several keV up to hundreds of keV [14]. The higher brilliance and the coherence of the present-day synchrotron radiation sources available at storage rings or free-electron lasers (FEL) allow fewer-shot measurements for non-reproducible objects. Storage rings, owing to their peak brightness being several orders of magnitude lower than the FEL sources, typically operate at lower flux conditions. Therefore, even though a typical measurement scheme may involve multiple pulses, the mean number of photons deployed in each shot is comparatively lower.
As a consequence, the propensity of scattering events leading to ionization during individual pulse interaction, per atom remains lower compared to FEL sources. The higher brilliance of the latter, which aids the signal observation, also induces scattering processes that involve multiple, correlated photon-matter interactions. In this paper, we focus on the multiphoton nature of the diffractive scattering at the storage-ring based synchrotron radiation sources [15,16]. In relevant situations, the full statistical distribution of the scattered photons involving all the higher-order moments is required for extracting information about the full set of diffractive scattering events. The detection schemes based on single-photon monitoring, which effectively observes the reduced many-particle density operator of the scattered photon field, may not be an adequate quantity. It calls for a theoretical formulation that can systematically account for a multiphoton scattering and detection scenario and, more importantly, examines the physical situations which may enable the validity of a truncated subspace monitoring of the signal. Investigating the cases where a particular reduced space detection may hold also offers rationale in terms of the underlying structural features of specific materials. In typical measurements, pixel-resolved diffracted photon intensities are sampled in the reciprocal space of the spatial variables. Following several acquisition steps, these outcomes are composed together to generate diffraction patterns which, in turn, correspond to the parametric variations of the electronic properties of scattering centers in real space. The signal is routinely interpreted in terms of the time-dependent, many-body electronic density of the scattering centers [17][18][19]. Consequently, while inverting the diffraction pattern, the commonly deployed inversion schemes assume that the diffraction patterns have originated solely from the time-dependent electronic density snapshots. Previous studies have demonstrated a generalization of this phenomenology by establishing that in the one-photon diffractive scattering measurements, the signal is in fact proportional to the many-body electronic density-density correlation function [9,[20][21][22]. In this communication, we establish a further generalization by extending it to the case of n-photon diffractive scattering scenario.
In Sec. II, we introduce the Hamiltonian and develop the integral equation-based formulation for multiphoton diffractive scattering. We present a detailed analysis for the interpretation of such signals involving the higher-order correlation functions of the incoming sources, electronic density, and detection modes. Subsequently, in Sec. III, we analyze various plausible physical features responsible for correlation signatures in the signal. We systematically develop a hierarchy of approximations, discuss their validity, and recover, in the asymptotic limits, the commonly used expressions. A set of physical conditions that may justify the usage of such a limit are described. It will be shown that the factorization of correlation functions, often applicable, may not be a generic prescription for interpreting the observed signal. The paper concludes with a perspective and outlook on the theory, a summary of the assumptions, presented in Sec. IV. detection, which is considered in the paper. The scattered modes contain a wealth of information regarding correlations among the incoming field-modulated scattering centers.
These can be revealed by focusing on the multiphoton diffraction characteristics.

A. Hamiltonian
Multiphoton diffraction phenomenology involving interaction between electronic matter and high-energy x-ray photons in the nonrelativistic limit is described by the minimal- In the above, the first two terms in the H matter correspond to the electronic matter modes.
The first term is the one-body operator consisting of electron kinetic energy and the electronnuclear potential, whereas the second term is a two-body operator describing the electronelectron interactions. We introduced the electronic field operators ψ(σr) (associated spin index is denoted by σ) which follow fermionic anti-commutation relations {ψ(σr), int containing n(r)A 2 (r) [where n(r) = σ ψ † (σr)ψ(σr) is the electronic density operator] dominantly contributes to the scattering processes in the high-energy regime away from resonance, often taking place beyond the ionization thresholds. In the appropriate parameter regime, the role of these terms in the dynamics becomes relevant with varying degrees of importance, and consequently, the relative contributions towards the scattering cross-sections vary. In general, the scattering events described by the Hamiltonian component H (2) int effectively take place from the electronic degrees of freedom driven by the H (1) int . Often, the ionization processes triggered by the latter lead to the onset of the radiation damage [26]. Restriction on the fluence of the incoming sources and plausible sparsity of the distribution of the scattering centers allows for a parameter window within which the diffraction remains largely unaffected by ionization events. We focus on the diffractive scattering processes particularly suitable for structural imaging and allow this analysis to be restricted to the scatterings induced by H (2) int . In the following, we separate the term H (2) int into initially occupied pump modes and initially unoccupied signal modes. Subsequently, as an approximation, we omit the scattering processes originating from pure pump or signal modes yielding, n(r)A 2 (r) = (A s (r)+A p (r)) 2 n(r) = (A 2 s (r)+A 2 p (r)+2A s (r)A p (r))n(r) ≈ 2A s (r)A p (r)n(r). Diffractive scattering via pure pump or signal modes may allow, as proposed in [22], selfheterodyning of the signal. Further, we adopt the interaction representation defined by is the propagator. Using the mode-expansion of the vector potential in terms of plane-waves, and assumed rotating wave approximation, we can express the relevant Hamiltonian component as H (2) int (t) = α 2 dr ks,kp,µs,µpc √ ω ks ω kp |ε * ks,µs · ε kp,µp | a † ks,µs a kp,µp e +iq·r n(r, t) e −iω k t + h.c..
In the above, we defined the polarization vector of the mode (indexed by k, µ) as ε k,µ and usedc = 2π/V α 2 , where V , and ω p (ω s ) are the mode volume and frequency of the incoming (scattered) mode, respectively. The mode creation (annihilation) operators a † k,µ (a k,µ ) follow the bosonic commutation relations [a k,µ , a † k ′ ,µ ′ ] = δ kk ′ δ µµ ′ . The difference frequency between the incoming and scattered photon modes is denotedω k = ω kp − ω ks . In Eq. (4), we note the appearance of the electronic density operator and the position-dependent phase term in the exponent, which incorporates the difference between the incoming and scattered photon momenta,q = (k p − k s ).
The elastic components of the diffraction signal, where the difference between the moduli of the incoming and scattered photon momenta is zero, provide structural information. In comparison, the inelastic components where the incoming and scattered photon momenta differ in magnitude contain structural and dynamical information regarding the electronic excitations. The applicability of elastic scattering extends to techniques such as small-angle x-ray scattering (SAXS) or wide-angle x-ray scattering (WAXS) which involves averaging over the structural disorder. They become crucial for interfaces, amorphous materials with short-range order and molecular aggregates, and proteins in the solution phase.

B. Multiphoton diffraction signal and classification
In the high-intensity photon-matter interaction at the relevant energy regime, the typical one-pixel detection scenario monitors fewer number of modes among the full set of diffracted modes. As envisaged in Fig. 1, the scattered photonic modes would impinge on the pixelated array detector where photon intensities are registered. The pixel-resolved photon intensity distribution, obtained by averaging over repeated measurements that monitor individual pixels separately, constitutes the diffraction pattern. We present a theoretical framework that is capable of describing the generalized diffraction and detection scenario appearing in such cases. Towards that goal, we start by defining the diffractive scattering signal as the time-space integrated intensity of the electric field at the pixelated detector given by, S (n) out = dt n dr n · · · dt 1 dr 1 ⟨E (−) (r n , t n ) · · · E (−) (r 1 , t 1 )E (+) (r 1 , t 1 ) · · · E (+) (r n , t n )⟩.
Here, the integral expression describes the observable corresponding to the n-photon diffraction signal while the operator expectation is taken over the final state of the combined system consisting of the electronic matter, driving field, and detector. The signal, defined in this manner, can describe both the standard detections where the diffraction pattern is generated via averaging over several observed pixel distributions and the nonstandard situations where multi-pixel coincident detection is possible. The latter remains a forward-looking scenario due to current technological limitations. However, both cases can be described by employing different averaging procedures. In the Keldysh-Schwinger formulation, the operator expectation value can be evaluated by accounting for the forward-backward evolution of the combined system [27][28][29][30][31]. In other words, the final state can be systematically expanded in the powers of the interaction Hamiltonian to yield different scattering configurations. Each of these configurations corresponds to a particular order of the photon-matter scattering events. Since, the operator defined above is related to the n-photon-matter inter-action events, an expansion to the n-th order in H (2) int (t) for the bra and ket generates the desired signal. A compact expression for the signal can be obtained by collecting individual path-ordered operators and averaging them over the initially correlated combined state ϕ in .
Furthermore, if a factorizibility assumption on the initially correlated state is made, it implies ϕ in ≡ ϕ matter,in ⊗ ϕ field,in where the term ϕ matter,in represents the state of the correlated electronic matter prior to the scattering events, and ϕ field,in represents the initial state of the incoming photonic sources. For the latter, it has been assumed that the scattered modes are in a vacuum while the rest of the modes are in an arbitrary field state. The prior choice of neglecting H (1) int during the expansion signifies the absence of externally induced electronic current. The lack of nontrivial field-induced correlation among the electronic modes justifies the factorized form of the initial state. In this limit, the operator expectation value defined in Eq. (5) can be expressed as In the above, the arguments of the operators have been distinguished to indicate the time and path ordering. In the expression, we introduced a set of multipoint, space-time dependent correlation functions for the sake of brevity. In particular, we defined the correlation function of the scattered photon and detector modes. It can be expressed as the sum of terms consisting of the products of elementary two-point correlation functions given by Here the subscript d ∈ {1 · · · n} and j ′ , j ′′ ∈ {1 ′ · · · n ′ ; 1 ′′ · · · n ′′ } are associated with the detection and scattering events, respectively. The elementary two-point correlation functions are expressed as correlation function, which can be expressed as The spatio-temporal coherence and statistical properties of the incoming photonic field encoded in this term dictate the nature of the correlation between the scattering events [32][33][34][35].
Hence, a prior characterization of its functional nature is essential for carrying out controlled diffractive scattering as well as a reliable posterior data inversion. Furthermore, we defined the electronic density correlation function composed of the space-time dependent density operators n(r, t) expressed as Since at this stage, the electronic density operators do not assume any ad hoc separation in terms of a noninteracting and a correlated part, in principle, this function retains the full set of information regarding the correlation and fluctuation properties of the electronic degrees of freedom [36].  The scattering events that are spatially close (i.e., at length scales shorter than the coherence length) have been denoted by the orange-colored vertex, and the ones that are spatially far apart (i.e., at length scales larger than the coherence length) are denoted by the blue-colored vertex. These lead to different subclasses of diagrams.

B. Role of factorizable incoming photonic field correlation function
In certain cases, the spatio-temporal coherence properties of the incoming photonic sources become the dominant governing factor that determines the nature of correlated scattering processes. Typically, the incoming photonic modes are scattered from the distribution of scattering centers. Information regarding the time-dependent relative spatial distribution of the scattering centers is contained in the correlation and fluctuation characteristics of the electronic densities. The latter alters the statistics of the scattered photonic modes and gives rise to diagnostic interference patterns at the detector. If the temporal and transverse coherence lengths of the incoming photonic modes are shorter than the spatiotemporal coherence lengths associated with the internal electronic processes, the multipoint incoming field correlation function can be represented in a factorized form without loss of generality. In its simplest form, the factorized form can be expressed in terms of a sum containing products of two-point correlation functions as given below Here out,1 = dt n dr n · · · dt 1 dr 1 dr n ′′ dt n ′′ · · · dr 1 ′′ dt 1 ′′ dr n ′ dt n ′ · · · dr 1 ′ dt 1 ′ × D (n) s (r n , t n , · · · r 1 , t 1 ; r 1 ′ , t 1 ′ · · · r n ′ , t n ′ ; r n ′′ , t n ′′ · · · r 1 ′′ , t 1 ′′ ) p,w (r j ′ , t j ′ ; r j ′′ , t j ′′ ) . we focus on the previously introduced n-point density correlation function, in Eq. (10). This expression can not be reexpressed in terms of the lower-order factorized counterparts due to the many-particle nature of the function. However, exceptions may arise due to the phys-ically motivated assumptions leading to the separation of space-time variables. One such condition may originate while considering the joint spatio-temporal characteristics of the incoming field and electronic correlation properties. To describe such a scenario, we introduce the notion of material-specific spatio-temporal coherence length and define the corresponding parametric variable λ(r 1 · · · r n , t 1 · · · t n ), which is a function of the distribution of scattering centers in the real space. Using the definition, without any loss of generality, we may recast the expression for the density correlation function as where K 0,g (r j ′ , t j ′ , r j ′′ , t j ′′ ) is the reference density correlation function defined in the limit where the multipoint density correlation function factorizes. The associated exponent de- properties of the electronic wave functions, the interatomic potentials, and the nature of the disorders. The latter two features are particularly important for liquids and amorphous materials, respectively. In the absence of significant inter-particle interactions, the coherence lengths become comparatively short-ranged. As a consequence, the cumulant function may decay faster than the reference density correlation function. In such cases, we would replace the cumulant function with a c-number and arrive at a partially factorized limit. In general, any physical assumption that lead to the systematic separation of multipoint space-time variables would result in the factorization of the electronic correlation functions in terms of their lower-order counterparts. A diagrammatic interpretation that may guide such intuitive length scale separation has been described in Fig. 4. Therein, a classification scheme based on whether the separation between the scattering centers is larger (shorter) than the coherence length scale parameter is shown to generate several subclasses of diagrams for one particular configuration of the detection and incoming field. Furthermore, in cases where the properties of the incoming field correlation allows separation of the integral variables into non-interacting space-time clusters, the multipoint electronic correlation functions undergo factorization that is dominantly governed by the same. It leads to the field and electronic density correlation functions being associated with the same set of combinatorial indices, w.
Furthermore, due to the weak presence of correlation in the incoming field, the cumulant functional in the signal expression can be trivially set to unity, which yields the completely factorizable limit. In that limit, the reference density correlation function can be presented as a sum of the products of renormalized two-point density correlation functions. The terms that form the relevant combinatorial set are explicitly given in Eqs. (A5) and (A7). These terms and the ones given in Eqs. (A6) and (A8) put together would constitute the larger subset {g} indicated in Eq. (13). Incorporating the assumptions stated above, we recover the signal given in Eq. (6) in the factorized limit as S (n) out,2 = dt n dr n · · · dt 1 dr 1 dr n ′′ dt n ′′ · · · dr 1 ′′ dt 1 ′′ dr n ′ dt n ′ · · · dr 1 ′ dt 1 ′ × D (n) s (r n , t n , · · · r 1 , t 1 ; r 1 ′ , t 1 ′ · · · r n ′ , t n ′ ; r n ′′ , t n ′′ · · · r 1 ′′ , t 1 ′′ ) Even though the signal expression is evaluated in the factorized limit, we note that the correlation contributions are combinatorial in nature. The information remains encoded into the various products of the 2-point correlation functions involving electronic densities, incoming and scattered field operators. In order to make the role of the involved photonic field modes explicit, we insert the mode-decomposition of the incoming, scattered, and detected field modes in Eq. (14) and recast the signal expression as S (n) out,2 = (α 4 ) n dt n dr n · · · dt 1 dr 1 dr n ′′ dt n ′′ · · · dr 1 ′′ dt 1 ′′ dr n ′ dt n ′ · · · dr 1 ′ dt 1 ′ where we used the notation {µ s j , µ p j } and {k s j , k p j } to indicate the set of polarization values µ s 1 · · · µ sn , µ p 1 · · · µ pn and momentum values k s 1 · · · k sn , k p 1 · · · k pn lying in either of the branches over which the summation would be carried out. We note that the factorized Moreover, we consider a situation where the bandwidth of the incoming photonic source encompasses the energy range of the electronic dynamics of interest. Further, assuming that the envelope of the incoming pulses has a narrow bandwidth and a smaller angular spread, it leads to a case where the largest spectral weights of the field modes are around a mean wavevector k p j and, polarization µ p j . It allows the replacements, ω kp j ′ ω kp j ′′ ≈ ω kp j and ε kp j ′ ,µp j ′ = ε kp j ′′ ,µp j ′′ ≈ ε kp j ,µp j . If the incoming pulses have negligible spatial dependence compared to the material specific coherence length scales, a replacement of the position dependence of the field profile by its value at r 0 can be made. This leads to the factorization of the field intensities, keeping only the terms which are commensurate with the approximation on electronic density operators. We also assume, ω ks j ′ = ω ks j ′′ ≈ ω ks j and ε ks j ′ ,µs j ′ = ε ks j ′′ ,µs j ′′ ≈ ε ks j ,µs j . Alongside, it is convenient to move to a new set of time variables defined ast j = (t j ′′ + t j ′ )/2 andτ j = (t j ′′ − t j ′ ). Finally, we arrive at the expression S (n) out,3 = (α 4 ) n dt n dr n · · · dt 1 dr 1 dr n ′′ · · · dr 1 ′′ dr n ′ · · · dr 1 ′ × dt n dτ n · · · dt 1 dτ 1 (1) 0 (r n ′ , t n ′ ; r n ′′ , t n ′′ ) exp(−iq 1 ′ · r 1 ′ · · · − iq n ′ · r n ′ + iq n ′′ · r n ′′ · · · + iq 1 ′′ · r 1 ′ I (r 0 ,τ 1 ,t 1 ) · · · I (r 0 ,τ n ,t n ) exp(−iω k 1τ 1 · · · − iω knτn ).
Here, F 1,j (r 0 ,t j ) = |E 0 (r 0 ,t j )| 2 denotes the intensity of the incoming field and F 2,j (r 0 ,τ j ) = exp (−Γ 0,jτ 2 j ), where Γ 0,j = 2 ln 2/τ 2 0,j , with τ 0,j being the temporal width of the pulse. A more constrained approximation can be invoked at this point by factorizing the two-point density correlation functions in terms of products of means of electronic density operators.
Additionally, we assume that the scattered modes are assigned to a particular set of predesignated detection modes. It amounts to neglecting the combinatorial nature of the corresponding correlation function. Moreover, if the mean of electronic density operators does not depend significantly on the temporal delay variablesτ j , we may drop the argument in the expression and proceed to carry out the corresponding Fourier integral. If the detected frequencies are assumed to lie within the vicinity of the incoming frequencies, i.e., ω ks ≈ ω kp , the Fourier integral of F 2,j (r 0 ,τ j ) contains a frequency domain Gaussian function that is finite over a negligibly narrow range of frequencies. The exponent in the latter can be replaced by unity, i.e., π/Γ 0,j exp(−ω 2 k j /Γ 0,j ) ≈ π/Γ 0,j . Using these approximations, the expression in Eq. (18) can be symmetrized. For an approximation that drops time-dependence of the means of density operators entirely, the two temporal integrals decouple and can be carried out separately. These hierarchical sets of approximations are presented in the form Eq. (20b), for the case of n = 1, yields the diffraction signal which frequently appears in the literature. That expression may be thought of as a limit recoverable following a set of approximations in which dressing of electronic densities by classical-driving has been assumed. Traditionally, this limiting case, for the lowest order, was derived by considering the first-order perturbation in the photon-matter interaction under the Born approximation.
Such considerations are solely based on the arguments regarding the field-matter interaction strength. The Born-approximation for the multiphoton case would imply the mutual independence of the scattering events. The presumptive application of such phenomenology leads to the factorizability of the scattering cross section in terms of mutually independent contributions. Within our treatment, this limit is recovered via a more rigorous and systematic procedure. In passing, we note that a decomposition of the average electron density in terms of contributions from atomic centers located at R p can be used to re-express the density scattering factors in terms of a linear combination of contributions at positions p arising from atomic density scattering factors. Such exercises have been taken up in several articles previously, notably in [13,25]. Following algebraic simplifications, it was concluded that the contributions to the diffractive imaging can be decomposed in terms of one-and twoparticle contributions. In the general framework adopted in this article, the scattering target is treated as a single object, even it consists of several constituents e.g., many molecules of the same species. Therefore, we have refrained from adopting such a decomposition in terms of one-and two-particle contributions.

IV. CONCLUSIONS AND OUTLOOK
In this article, we formulated a theory of multiphoton x-ray diffraction, with a focus on the role of higher-order photonic and electronic correlations. The correlations and fluctuations associated with the outgoing photonic modes were shown to contain useful and sometimes crucial information for the electronic matter correlations. We also pointed out the physical origins of the correlation resources that may be exploited to design controlled diffractive imaging measurements. For specific physical situations, these correlation-generating resources may act in unison to give rise to novel signal components that contain a wealth of information regarding the underlying material structure. These measurements may become crucial for designing advanced structural probes for complex correlated matter [37][38][39][40], e.g., quantum materials with novel topological properties [41], strange metals [42][43][44] and quantum materials and molecules in engineered electromagnetic environments [45,46]. In However, we note that the incorporation of dynamical pathways related to the abovementioned features may restrict the application of the factorization hierarchy. In practice, distinct signatures that point towards the existence of factorizable limits may need the full suite of correlation analyses. Simulation of the multipoint electronic correlation functions is a challenging task [47][48][49][50]. Rate equations involving time-dependent electronic configurations have shown early promise [51][52][53], which may provide a multiscale simulation toolbox towards the goal in the specified parameter regime. An extension of the rate-equation framework incorporating the inter-configuration coherence is underway. Such an algorithm, still numerically expensive, may provide a realistic way to investigate the factorization hierarchy.
This, in turn, would require unique data-inversion algorithms.
In recent years, theoretical proposals concerning the usage of novel multimode radiation sources [22] as investigative tools in structural biology have been made. The presented integral equation based formulation can easily accommodate them via the basic modification of the incoming field correlation functions. Moreover, optimally tailored incoming photonic field properties that selectively amplify specific parts of the signal can be found by employing control algorithms [54][55][56][57][58].
The inaccessibility of the phase in the intensity-detected scattering measurements has been a persistent problem that prevents the unambiguous inversion of the signal. It has been addressed in the case of crystalline and non-crystalline finite-sized samples via oversampling and iterative algorithms [59,60]. In a set of forward-looking proposals, which may become challenging to implement at the current state-of-the-art, Mukamel et al. and others have proposed to resolve the problem with the help of phase-sensitive coincidence-detected diffractive scattering measurements [22,61]. This proposal holds significant promise if it can be ingeniously combined with advanced detector technologies. The theoretical development in this paper requires minimal modification to incorporate and extend the proposed approach.
One of the consequences of monitoring fewer detection modes was that the available reconstruction algorithms assume the instantaneous electronic densities to be the central quantities for the posterior data inversion. It was shown that the multipoint density correlation functions are the generalized quantities carrying information about the distribution of the scattering centers. The presented paper lays the groundwork for the identification of specific terms in the matter correlation functions related to the relevant features in the diffraction signals. The analysis, possibly, sets the ground for more accurate interpretations of the snapshot diffraction measurements, and reconstruction algorithms beyond density-based heuristics [19,[62][63][64][65][66][67] and provides a theoretical framework for further developments pertinent to the single-particle imaging using FEL radiation sources [26,68]. One of the promising applications of such an elaborate description lies in time-domain structural biology. Work in this direction is already underway by the authors. The correlation function that indicates the combinatorial detection pathways, for n = 2, can be expressed as In the above, we note that the scattered modes can be paired with the detection modes in four (i.e., m = 4) distinct ways. Similarly, for the case of n = 3 we have m = 36 terms which can be expressed as, s (r 2 , t 2 ; r 3 ′ , t 3 ′ )D (1) s (r 3 , t 3 ; r 2 ′ , t 2 ′ ) + D (1) s (r 1 , t 1 ; r 2 ′ , t 2 ′ )D (1) s (r 2 , t 2 ; r 3 ′ , t 3 ′ )D (1) s (r 3 , t 3 ; r 1 ′ , t 1 ′ ) + D (1) s (r 1 , t 1 ; r 3 ′ , t 3 ′ )D (1) s (r 2 , t 2 ; r 1 ′ , t 1 ′ )D Following an analogous exercise, the factorization of the incoming field correlation functions can be performed as well. In doing so, only the 2-tuples arising from the pairings across branches were kept. This particular approximation neglects the effect of extended spatiotemporal field coherence on densities. For the two-point incoming field correlation function it yields two terms (i.e., w = 2) given bỹ = K (1) 0 (r 1 ′ , t 1 ′ ; r 1 ′′ , t 1 ′′ )K (1) 0 (r 2 ′ , t 2 ′ ; r 2 ′′ , t 2 ′′ ) + K (1) 0 (r 1 ′ , t 1 ′ ; r 2 ′′ , t 2 ′′ )K (1) 0 (r 2 ′ , t 2 ′ ; r 1 ′′ , t 1 ′′ ).