## 2 Quantum-Field Theory for Photonic Transport
### Bethe-Salpeter Equation for Photonic Transport in Samples of Disordered Active Mie Resonators  
We consider the propagating photonic intensity, see Figure 1c, as the field field correlator \(\Phi\) which is derived selfconsistently by the _Bethe-Salpeter equation_. The propagating electromagnetic wave in disordered random media is described by the single-particle Green's function Equation (3) that solves the non-linear wave equation, Equation (1). \[G_{k}^{\omega}=\frac{1}{\epsilon_{b}(\omega/c)^{2}-|\vec{k}|^{2}-\Sigma_{\vec {k}}^{\omega}} \tag{3}\]  
The scatterers are embedded in a background matrix, with the dielectric function \(\epsilon_{b}\), e.g., in air. The granularity of the random medium is implemented in the form of a spatially dependent potential of the permittivity \(\epsilon_{s}(\vec{r})\). As a result of our theoretical framework we derive the laser gain self-consistently and space-resolved. The order of non-linearity is a matter of the numerical self-consistency as the independent complex semiconductor scatterers are well described by the T-matrix [40] for spherical Mie resonators of matter with a self-consistent complex refractive index. We use a model of independent Mie scatterers, Figure 1b, here, in order to derive the self-energy contribution \(\Sigma_{k}^{\omega}\) for extended photo-active particles in Equation (3). The Mie scattering coefficients of n-th order are known from the literature [41] to have the following form  
\[a_{n} = \frac{m\Psi_{n}(my)\Psi_{n}^{\prime}(y)-\Psi_{n}(y)\Psi_{n}^{ \prime}(my)}{m\Psi_{n}(my)\xi_{n}^{\prime\prime}(y)-\xi_{n}(y)\Psi_{n}^{ \prime}(my)} \tag{4}\] \[b_{n} = \frac{\Psi_{n}(my)\Psi_{n}^{\prime}(y)-m\Psi_{n}(y)\Psi_{n}^{ \prime}(my)}{\Psi_{n}(my)\xi_{n}^{\prime\prime}(y)-m\xi_{n}(y)\Psi_{n}^{ \prime}(my)}.\]  
In this notation \(m=\frac{n_{s}}{n_{b}}\) is the relative refractive index between the scatterer and the background matrix, \(y=\frac{2\pi m_{s}\Gamma_{s}}{\lambda}\) is the size parameter depending on the scatterers' radius \(r_{S}\) as well as on the transport wavelength \(\lambda\). \(\Psi_{n}\), \(\xi_{n}\) are Riccati-Bessel functions. The refractive index \(n_{s}(\Phi)\) obtains a higher-order non-linearity as consequence of the self-consistency with respect to \(\Phi\), \(\text{Im}\epsilon_{s}(\Phi)\). The independent scatterer approach works very well in dense random media of filling of 35-55%. Positional dependent enhancement effects due to the high filling fraction are effectively mapped on the dynamics of \(n_{s}(\Phi)\) of the single scatterer in non-linear response, and we use the exact Mie scattering solution. The Green's function, Equation (3), for \(\Psi_{\omega}\), solves the wave equation Equation (1) and builds up the intensity field in the sample in the form of the field correlator \(\Phi_{\vec{k}\,\vec{k}^{\prime}}^{\omega}(\vec{Q},\Omega)\). The four point correlator  
\[\Phi_{\vec{k}\,\vec{k}^{\prime}}^{\omega}(\vec{Q},\Omega)=\langle\hat{G}_{\vec {k}_{+}\,\vec{k}_{+}^{\prime}}^{R}(\omega)\hat{G}_{\vec{k}_{-}\,\vec{k}_{-}^{ \prime}}^{A}(\omega)\rangle \tag{5}\]  
is thus expressed in terms of retarded and advanced Green's function \(\hat{G}^{R/A}\) denoting the field and its time-reversal. We utilize the transformation of coordinates in center-of-motion \((\vec{Q},\Omega)\) and relative \((\vec{k},\omega)\) momenta and frequencies [34]\(\vec{k}_{\pm}=\vec{k}\pm\vec{Q}/2\) and \(\omega_{\pm}=\omega\pm\Omega/2\). Considering the slab geometry to be extended within the \((x,y)\)-plane and finite in \(z\)-direction, the full Fourier transform as indicated for infinite samples [31] is replaced by a partial Fourier transform following the argument of the separation of the scales for the field \(\Psi\) and the intensity \(\Phi\).