## 3 Gauss-Maxwell Beams  
### General Considerations on the Propagation Procedure  
Gaussian beams are exact solutions to the paraxial Equation (5), as seen above. They are typically interpreted as the amplitude associated with electric and magnetic field vectors, transverse to the beam axis at any \(z\) value and everywhere polarized in the same direction. However, they are not exact paraxial solutions to [PERSON]'s equations. Solutions that satisfy this requirement were found by [PERSON] et al. [28], who called them Gauss-Maxwell beams (for other alternative but equivalent procedures, see References [35; 36; 37; 38]). These beams can be determined according to the prescription that we shall describe now, and that will be considered from now on. Consider the general form of a scalar solution to the paraxial Helmholtz equation, Equation (5),  
\[\psi(\mathbf{r}_{\perp};z)=\hat{\mathcal{U}}(\hat{\mathbf{r}}_{\perp},\hat{ \mathbf{p}}_{\perp})\psi(\mathbf{r}_{\perp};0), \tag{28}\]  
where  
\[\hat{\mathcal{U}}(\hat{\mathbf{r}}_{\perp},\hat{\mathbf{p}}_{\perp})=e^{iz \hat{\mathbf{p}}_{\perp}^{2}/2k} \tag{29}\]  
is the free-space propagator, with \(\hat{\mathbf{r}}_{\perp}=\mathbf{r}_{\perp}\) and \(\hat{\mathbf{p}}_{\perp}=-i\
abla_{\perp}\) denoting the transverse position and momentum operators, respectively. According to Equation (28), the spatial distribution of electromagnetic energy can be easily determined by considering an input beam \(\psi(\mathbf{r}_{\perp};0)\) at \(z=0\) and then evolving it along \(z\) by means of the propagator (29). According to the discussion in the previous section, this propagator provides us with the transverse distribution of electromagnetic energy at each \(z\)-plane. Although this procedure is for scalar fields, a similar procedure can also be followed for the vector fields governed by [PERSON]'s equations by recasting the electromagnetic field as a six-component vector field (to some extent, analogous to the so-called Riemann-Silberstein electromagnetic vector [8]),  
\[\mathbf{F}(\mathbf{r}_{\perp};z)=\frac{1}{\sqrt{2}}\left( \begin{array}{c} \sqrt{\epsilon_{0}}E_{x}(\mathbf{r}_{\perp};z)\\ \sqrt{\epsilon_{0}}E_{y}(\mathbf{r}_{\perp};z)\\ \sqrt{\epsilon_{0}}E_{z}(\mathbf{r}_{\perp};z)\\ \sqrt{\mu_{0}}H_{x}(\mathbf{r}_{\perp};z)\\ \sqrt{\mu_{0}}H_{y}(\mathbf{r}_{\perp};z)\\ \sqrt{\mu_{0}}H_{z}(\mathbf{r}_{\perp};z)\\ \end{array} \right), \tag{30}\]  
where \(\mathbf{H}=\mathbf{B}/\mu_{0}\) is used instead of \(\mathbf{B}\) for simplify, as it will be seen below. The vector field \(\mathbf{F}(\mathbf{r}_{\perp};z)\) at the output plane \(z\) arises after propagating a distance \(z\) an input vector field \(\mathbf{F}(\mathbf{r}_{\perp};0)\), with the evolution being described by a certain operator that has to be determined. To this end, notice that if the evolution of an input scalar field is described by the propagator (29), then the evolution of \(\mathbf{F}\) should be described by a 6 \(\times\) 6-matrix operator, henceforth denoted by \(\underline{\mathbb{U}}\). This operator is obtained by replacing the transverse position vector \(\mathbf{r}_{\perp}=(x,y)\) in the scalar operator (29) by a more general position vector operator [28], with matrix elements given by  
\[\mathbb{R}_{\perp}=\mathbf{r}_{\perp}\mathbb{I}+k^{-1}\mathbb{G}_{\perp}=(x \mathbb{I}+k^{-1}\mathbb{G}_{x},\ y\mathbb{I}+k^{-1}\mathbb{G}_{y}). \tag{31}\]  
In other words, we have a transformation  
\[\hat{\mathcal{U}}(\hat{\mathbf{r}}_{\perp},\hat{\mathbf{p}}_{\perp})\psi( \mathbf{r}_{\perp};0)\longrightarrow\underline{\mathbb{U}}(\hat{\mathbf{r}}_{ \perp}\mathbb{I}+k^{-1}\mathbb{G}_{\perp},\hat{\mathbf{p}}_{\perp}\mathbb{I}) \mathbf{F}(\mathbf{r}_{\perp};0), \tag{32}\]  
where \(\mathbb{I}\) is a \(6\times 6\) identity matrix, and the \(\mathbb{G}\)-matrices are defined as  
\[\mathbb{G}_{x}=\frac{1}{2}\left( \begin{array}{cc}-\mathbb{S}_{2}&\mathbb{S} _{1}\\ -\mathbb{S}_{1}&-\mathbb{S}_{2}\\ \end{array} \right),\qquad\mathbb{G}_{y}=\frac{1}{2}\left(