as $$\label{WDF}
W_{\Gamma}(x,p) 
= \int {\rm d}s \, \e^{-\rmi k p s}\, \Gamma\left(x+\frac{s}{2}, x-\frac{s}{2}\right)\, ,$$ with back transformation given by $$\label{WDFback}
  \Gamma(x_1,x_2) =  \left ( \frac{k}{2 \pi} \right )^{d-1} \int {\rm
    d}p \, 
\e^{\rmi k (x_1-x_2)p} \; W_\Gamma \left(\frac{x_1+x_2}{2}, p \right).$$ Here, $k$ represents a wave number such as introduced in (\[helm\]), (\[BoT\]) and $d$ is the dimension of the full, interior problem. It is shown in \[app:A\], that for sufficiently smooth initial phase space distributions $W_{\Gamma_0} = \rho_0$, the averaged Wigner transform of $\hat{K}_n$ in (\[K\]) is given in terms of the classical flow equations, see also [@GCT15; @MWH01].

Classical phase space densities are driven by the phase space dynamics, in our case a boundary map or more generally the Poincaré map of an SOS. The initial density $\rho_0 = W_{\Gamma_0}$ can then be associated with a boundary density of rays arriving directly from a source distribution in the interior. Mapping the source ray density through subsequent reflections leads to the iterated densities $\rho_0\to\rho_1\to\cdots\rho_n\to\cdots$ which can be described in terms of the (linear) integral operator $\cal L$ defined in the lossless limit as [@book] $$\label{FP}
{\cal L} [\rho_n](X) = \rho_{n-1}(\varphi^{-1}(X)) =  
\int {\rm d} X' \delta(X-\varphi(X'))\rho_{n-1}(X')$$ (see below for 