## 3 Gauss-Maxwell Beams
### Linearly Polarized Gauss-Maxwell Beams
Once we have the correct paraxial vector field solutions to [PERSON]'s equations, we can compute the averaged paths that will be followed by photons, according to (46). In this regard, let us first start by the lowest level of approximation, namely the first order in \(\Delta\mathbf{r}_{\perp}\), and compute the time-averaged Poynting vector (45). Thus, neglecting the terms containing \((\Delta x)^{2}\), \((\Delta y)^{2}\) and \(\Delta x\Delta y\) in (60) and (61), we obtain  
\[\mathbf{S}=\frac{1}{2}\,E_{0}H_{0}|\psi|^{2}\left( \begin{array}{c}\frac{ \Delta x}{2k\sigma_{0}\sigma_{z}}\ \sin\varphi_{z}\\ \\ \frac{\Delta y}{2k\sigma_{0}\sigma_{z}}\ \sin\varphi_{z}\\ \\ 1\end{array} \right), \tag{66}\]  
where  
\[\sin\varphi_{z}=\frac{z}{2k\sigma_{0}\sigma_{z}}. \tag{67}\]  
The ray equation that follows from (66) is  
\[\frac{d\Delta\mathbf{r}_{\perp}}{dz}=\frac{z}{4k^{2}\sigma_{0}^{2}}\frac{ \Delta\mathbf{r}_{\perp}}{\sigma_{z}^{2}}, \tag{68}\]  
where the left-hand side term has been recast in terms of \(\Delta\mathbf{r}_{\perp}\) instead of \(\mathbf{r}_{\perp}\), without loss of generality (there is no effect on the \(z\)-derivative). As it can be noted, this equation corresponds to Equations (21) and (22) for the usual Gaussian beam. That is, we find that the energy spreads spatially in the same way as a standard Gaussian beam at the lowest level of approximation, although the polarization state changes along both the transverse and the longitudinal directions, even if it was specifically defined at the input plane. In order to determine the deviations from the usual Gaussian beam approach, let us consider the full expression for the fields (60) and (61). The computation of the time-averaged Poynting vector renders  
\[\mathbf{S}=\frac{1}{2}\,E_{0}H_{0}|\psi|^{2}\left( \begin{array}{c}\frac{ \Delta x}{2k\sigma_{0}\sigma_{z}}\ \sin\varphi_{z}\left[1+\frac{1}{2}\left(\frac{\|\Delta \mathbf{r}_{\perp}\|}{2k\sigma_{0}\sigma_{z}}\right)^{2}\right]\\ \\ \frac{\Delta y}{2k\sigma_{0}\sigma_{z}}\ \sin\varphi_{z}\left[1+\frac{1}{2} \left(\frac{\|\Delta\mathbf{r}_{\perp}\|}{2k\sigma_{0}\sigma_{z}}\right)^{2} \right]\\ \\ 1-\frac{1}{4}\left(\frac{\|\Delta\mathbf{r}_{\perp}\|}{2k\sigma_{0}\sigma_{z}} \right)^{4}\end{array} \right), \tag{69}\]  
where it can be seen that the \(z\)-component is diminished, precisely, in an amount equivalent to the energy flux going into the transverse directions. Keeping in the ray equation terms depending on the transverse displacement, \(\Delta\mathbf{r}_{\perp}\), up to third order, we obtain  
\[\frac{d\Delta\mathbf{r}_{\perp}}{dz}=\frac{z}{4k^{2}\sigma_{0}^{2}}\frac{ \Delta\mathbf{r}_{\perp}}{\sigma_{z}^{2}}\left[1+\frac{1}{2}\left(\frac{\| \Delta\mathbf{r}_{\perp}\|}{2k\sigma_{0}\sigma_{z}}\right)^{2}\right]. \tag{70}\]In order to get an idea of the importance of the second term within the square bracket, we can evaluate it taking into account the orders of magnitude of the different quantities involved in it. To this end, we can use values from the experiment reported in Reference [17]. Accordingly, in meters, \(\lambda\sim 10^{-6}\), while \(\sigma_{0}\), \(\sigma_{z}\) and \(\|\Delta\mathbf{r}_{\perp}\|\) are the three of the order of \(\sim 10^{-3}\). With these values, the correction term is of the order of \(10^{-6}\), which fully justifies that, in a good approximation, a scalar Gaussian beam description can be used instead of a more exact vector approach.