Bohmian-Based Approach to Gauss-Maxwell Beams

: Usual Gaussian beams are particular scalar solutions to the paraxial Helmholtz equation, which neglect the vector nature of light. In order to overcome this inconvenience, Simon et al. ( J. Opt. Soc. Am. A 1986 , 3 , 536–540) found a paraxial solution to Maxwell’s equation in vacuum, which includes polarization in a natural way, though still preserving the spatial Gaussianity of the beams. In this regard, it seems that these solutions, known as Gauss-Maxwell beams, are particularly appropriate and a natural tool in optical problems dealing with Gaussian beams acted or manipulated by polarizers. In this work, inspired in the Bohmian picture of quantum mechanics, a hydrodynamic-type extension of such a formulation is provided and discussed, complementing the notion of electromagnetic ﬁeld with that of (electromagnetic) ﬂow or streamline. In this regard, the method proposed has the advantage that the rays obtained from it render a bona ﬁde description of the spatial distribution of electromagnetic energy, since they are in compliance with the local space changes undergone by the time-averaged Poynting vector. This feature confers the approach a potential interest in the analysis and description of single-photon experiments, because of the direct connection between these rays and the average ﬂow exhibited by swarms of identical photons (regardless of the particular motion, if any, that these entities might have), at least in the case of Gaussian input beams. In order to illustrate the approach, here it is applied to two common scenarios, namely the diffraction undergone by a single Gauss-Maxwell beam and the interference produced by a coherent superposition of two of such beams.


Introduction
One of the most appealing features of geometrical optics, and also a remarkable and convenient advantage, is perhaps the fact that this model relies on the well-defined and very intuitive concept of ray. This concept, which simplifies the description and analysis of optical processes (e.g., imaging) and phenomena (e.g., reflection, refraction, miracles), has a direct physical meaning, as the path along which the electromagnetic energy in the form of light flows. Accordingly, it has helped us to understand how light goes from one place to another in a simple fashion, so it is not strange that descriptions and methodologies based on the notion of ray have become a valuable tool, for instance, to infer properties of the medium traversed by light (e.g., refractive index, curvature, thickness, etc.) or, conversely, based on such properties, to determine how geometrical arrangements of optical devices (e.g., lenses, mirrors, prisms, etc.) can be devised to act on or control light. However, when we move to the realm of electromagnetic or wave optics, this concept blurs up. That is, rays are still considered, but typically in a virtual sense, as an auxiliary tool to evaluate relatively complex integrals, as it is associated with such a beam propagate along the optical axis and spread across the transverse directions? Or, how does the polarization state influence the interference between two of such beams while they evolve along the optical axis? With the purpose to provide an answer to these questions, here we develop a hydrodynamic description for Gauss-Maxwell electromagnetic beams, developed in 1989 by Simon et al. [28], in terms of electromagnetic energy flow lines or rays. The motivation behind this analysis is to provide a ray-based description for localized electromagnetic fields (Gaussian-type beams), where both position and polarization degrees of freedom are present, and hence it can readily be applied to the analysis and interpretation of diffraction and interference experiments of the kind mentioned one. As is well known, usual Gaussian beams are exact solutions to the Helmholtz equation in paraxial form, but not to Maxwell's equation under paraxial conditions. The approach proposed by Simon et al., though, shows the specific functional form that the beam has to satisfy in order to be an exact solution to Maxwell's equations in paraxial form (in vacuum). In this regard, we would also like to mention that there are other approaches in the literature worth exploring with the same methodology (i.e., where the concept of Bohmian trajectory could be exported) in order to determine the corresponding ray equations, such as nondiffracting Helmholtz-Gauss beams [29], in the case of scalar paraxial fields, or vector Helmholtz-Gauss beams and related families [30,31], in the case of vector paraxial fields, solutions to Maxwell's equations.
In this work, in particular, we have focused on and investigated the formal aspects of the rays that describe the spatial development of the electromagnetic energy density in the case of single Gauss-Maxwell beams and interference in the coherent superposition of two Gauss-Maxwell beams. Furthermore, we also consider the limit where the behavior of these vector beams and their superpositions can be described, in a good approximation, by the usual Gaussian beams, hence neglecting the associated polarization state. Particular interest is thus paid to the transverse momentum, which is eventually the observable quantity in an experiment and, therefore, the quantity of interest to take the methodology here proposed to the analysis of real laboratory experiments. Accordingly, this work has been organized as follows. To be self-contained and, at the same time, to offer a wider contextualization of the work, in Section 2 we introduce a general overview of the treatment for standard scalar Gaussian beams, starting from a revision of some general aspects connected to monochromatic scalar fields in vacuum, and how a ray ("photon" trajectory) equation can be properly specified for these fields. Then, the case of Gauss-Maxwell beams is analyzed and discussed in Section 3, first for a linearly polarized beam and then extended to any polarization state. Section 4 is devoted to the extension of this methodology to the superposition of two coherent Gauss-Maxwell beams with the same polarization state and also with mutually arbitrary polarization states-a case that may occur if a beam is diffracted by a two slit and, immediately afterwards, each diffracted beam acquires a different polarization state. To conclude, a series of final remarks are summarized in Section 5.

General Aspects for Monochromatic Scalar Fields
Consider a general electromagnetic monochromatic scalar field in vacuum [2], with wavelength λ = c/ω. Within a generalized framework, the amplitude of this field, Ψ(r), can be specified as a complex-valued static (time-independent) field satisfying Helmholtz's equation, plus the corresponding boundary conditions. From this equation, the components of the wave vector k = (k x , k y , k z ) provide us with valuable information about how the field Ψ(r) changes spatially along each space direction, while the monochromaticity condition requires that the wave number, k = k = k 2 x + k 2 y + k 2 z = 2π/λ, remains constant. Because we are dealing with light, let us assume that there is a preferential direction for its space propagation. This direction is going to be referred to as the longitudinal or parallel ( ) direction, which eventually defines the system optical axis. Any other orthogonal direction is going to be denoted as the perpendicular or transverse (⊥) direction. Thus, in Cartesian coordinates, if we choose the longitudinal direction along the z-axis, the x and y axes specify the transverse (mutually orthogonal) directions. Consequently, the position vector can be recast as r = (x, y, z) = (r ⊥ , z) and the wave vector as k = (k ⊥ , k z ), with the wave number being k = k ⊥ 2 + k 2 z . If paraxial conditions are assumed, that is, space variations of the field are slower along the longitudinal direction than along the transverse ones, then k z k ⊥ and hence k ≈ k z . Accordingly, the scalar field Ψ(r) can be recast as a plane wave propagating along the z-direction, modulated by a spatially-dependent amplitude, that is, The substitution of the ansatz (3) into the Helmholtz Equation (2) renders where ∇ 2 ⊥ = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 is the transverse Laplacian. Physically, paraxiality implies that the amplitude of the scalar field has to smoothly change along the longitudinal (z) direction. This allows us to neglect the second derivative in z in Equation (4), thus allowing us to simplify this equation, which reads as As it can readily be noticed, this equation is formally equivalent to the time-dependent Schrődinger equation for a free particle of mass m in two dimensions [32]. Actually, in this latter case, if it is assumed that propagation along z is classical, we find a relationship between this coordinate and time: where p =hk. For simplicity and convenience, we have chosen z 0 = 0 in (6), although this is not necessary. It is worth highlighting that, in spite of its apparent simplicity, Equation (6) brings in two important consequences at a conceptual level: 1. This relation enables a direct switch from propagation in time of an extended wave, namely the scalar field Ψ(r), to space diffusion along the longitudinal (axial) direction of a particular "slice" of such a wave. More specifically, if we consider a transverse section or plane of the full wave [consider, for instance, that such a wave describes a pulse with amplitude ψ(r)] within the XY plane for a given value z 0 (in other words, the input plane z = z 0 , analogous to considering t 0 in a time-propagation), we shall obtain its spatial redistribution or accommodation to the corresponding boundaries at subsequent planes, with z increasing. Consequently, if the pulse has an extension along the z direction, considering different "slices" of the pulse (i.e., different z 0 planes), we may easily determine the shape of the pulse at a further distance by just combining all the resulting "slices" z f . Notice that this fact also allows to establish the validity limit for the approximation, which is going to remain correct provided dispersion along the z direction can be assumed to be negligible (i.e., as long as all "slices" travel with nearly the same speed). 2. On the other hand, such a relation enables a simple, direct link between wave optics and matter-wave optics, which arises from the formal relation established between Schrődinger's equation and the Helmholtz equation in paraxial form, both being parabolic differential equations describing the transport of the quantity Ψ (regardless of the nature of such a quantity, that is, whether it describes an electromagnetic amplitude or a probability amplitude) with an imaginary diffusion constant (this complex valuedness is precisely the fundamental trait that allows interference in the solution of these equations, but not in the heat equation, although it is also a differential equation of the parabolic class). This is also a rather convenient issue both analytically and computationally, because it explicitly shows that optical and matter-wave problems ruled by the same equation form have the same solution and, eventually, the same interpretations [33].
Taking into account point 2 above, Equation (5) can readily be recast in the form of a Schrődinger-type equation by replacing the coordinate z with the value indicated by the relation (6)-notice, however, that this gives rise to a two-dimensional Schrődinger equation, since the only space coordinates, after substitution, are x and y. This direct analogy can be taken a step further to introduce a guidance equation in the Bohmian form, namely with describing the space phase variations of the complex field amplitude ψ. Since dr ⊥ /dz ≈ k ⊥ /k, this equation allows us to describe the transverse distribution of electromagnetic energy in terms of streamlines by providing the corresponding initial conditions and integrating it along z [32]. Further, observe that also in this case, the connection between Equation (7) and the Bohmian equation for matter waves is straightforward, since the latter can be directly obtained by considering the relation (6) and the change S ψ → S/h, from the phase of the amplitude ψ to the phase of the matter wave (withh emphasizing the fact that S has units of action).

Gaussian Beam Propagation
Consider that at the input plane z 0 = 0, the electromagnetic field is described by a beam with Gaussian amplitude (on the z 0 plane), which may represent a Gaussian mode released from an optical fiber or the light coming out from a laser pointer. The beam (9) is centered at r ⊥,c = (x c , y c ) and its width is related to its waist, w 0 , by the simple relation σ 0 = w 0 /2. In this latter regard, we take here the conventional definition for the waist of a Gaussian beam, as its size at the point of its focus (here located at the input plane, z 0 = 0), which corresponds to the radius of the 1/e 2 irradiance contour at the plane (z 0 ) where the wavefront is flat. In order to provide some typical values, we can take them from the experiment reported in Reference [17]. In this experiment, a coherent superposition of two nearly Gaussian beams (in a good approximation [12]) is generated, with waists w 0 = 0.608 mm, wavelength λ = 943 nm, and their centers separated a distance d = 4.69 mm. The propagation of the amplitude (9) along the z-axis is obtained by acting on it with the free-space propagatorÛ , that is, which is equivalent to considering the integral (for a derivation of this expression, see Appendix A). Accordingly, the substitution of the ansatz (9) into the integral (11) renders where is a complex-valued norm factor, ϕ z is the well-known Gouy phase in optics (typical of Gaussian beams), andσ is an also complex-valued spread factor. The dispersion of the beam at a distance z from the slits is given by the expression related, by means of the simple relation σ z = w z /2, to radius of the 1/e 2 contour when the wave has propagated a distance z, As it can be noticed from (16), the expansion of the beam along the optical axis can be characterized in terms of the Rayleigh length, z R = 2kσ 2 0 .
which constitutes a critical or characteristic length that sets a separation between two dynamically different regimes regarding the beam expansion, namely the Fresnel or hyperbolic expansion regime and the Fraunhofer or linear expansion regime [34]. The intensity or irradiance obtained from (12) is The streamlines that describe the spatial propagation of this intensity along z can be obtained, according to Equation (7), from the equation of motion which, particularized to each transverse coordinate, leads to The integration of Equations (21) and (22) is straightforward, leading to where x 0 and y 0 denote the position of the optical streamline at the input plane z = 0, while its position on the output plane z is given by x z and y z . Notice that the solutions rendered by Equations (23) and (24) display radial symmetry with respect to the center of the Gaussian beam, (x c , y c ). Accordingly, if their initial positions distribute around this point in a circle of radius ρ 0 = r ⊥,0 − r ⊥,c , they will distribute on a plane at a distance z around a circle of radius ρ z = r ⊥,z − r ⊥,c that has increased its size in a factor σ z /σ 0 with respect to ρ 0 . This is better seen if, from the definitions of ρ 0 and ρ z , the two equations are rewritten in a more compact form as: Relatively close to the input plane, this increase is negligible and then, as z starts becoming important compared to the critical value z R , there is a quadratic dependence on z, as expected in a typical Fresnel regime. In contrast, asymptotically, far from the input plane, the increase of the radius is linear with z, in agreement with a propagation in a Fraunhofer regime. This can also be noticed from the expression describing the expansion rate of the beam, which can be obtained from Equations (21) and (22) [or also directly from (25)] and describes the beam expansion in relation to the actual beam size. This expression displays a maximum for z = z R . Thus, before reaching the position of this critical plane, the expansion rate is first negligible (for z ≪ z R ) and then starts increasing nearly linear with z (for z z R ). This implies that the radius will undergo a sort of initial fast boost that will accelerate its increase and then the expansion of the beam. However, after surpassing the critical plane, the rate falls with z, which will keep the radius increasing at a constant rate, and hence the divergence displayed by any ray with respect to the center of the input beam, consistent with the fact that far from the input plane (z → ∞), the radial speed, becomes a constant. This is in compliance with the linear increase of the cross section that we observe in Gaussian beams far from the source (e.g., a non-collimated laser beam). Furthermore, it is also worth noting that the expression for the expansion rate (26) is also the same as the wavefront curvature of the Gaussian beam, with its inverse providing the radius of curvature of the latter. In this regard, we have that wavefronts are flat (infinite radius of curvature) at the input plane z = 0 (i.e., at the beam waist) and also as z becomes large compared to z R (i.e., in the Fraunhofer regime). On the contrary, they undergo maximal curvature when z = z R (with the radius of curvature being 2z R ), just at the plane where the beam expansion rate (26) gets its maximum value, 1/2z R .

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 Maxwell's equations. Solutions that satisfy this requirement were found by Simon 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), is the free-space propagator, withr ⊥ = r ⊥ andp ⊥ = −i∇ ⊥ 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 ψ(r ⊥ ; 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 Maxwell'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]), where H = B/µ 0 is used instead of B for simplify, as it will be seen below. The vector field F(r ⊥ ; z) at the output plane z arises after propagating a distance z an input vector field F(r ⊥ ; 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 F should be described by a 6 × 6-matrix operator, henceforth denoted byÛ. This operator is obtained by replacing the transverse position vector r ⊥ = (x, y) in the scalar operator (29) by a more general position vector operator [28], with matrix elements given by In other words, we have a transformation where I is a 6 × 6 identity matrix, and the G-matrices are defined as The G-matrices satisfy the properties: for a, b, c = x, y, which can be proven taking into account the following matrix relations: with n ≥ 1 and where [·, ·] + denotes the anticommutator. Note that the action ofÛ on the vector field F can be accomplished in two steps: 1. First, there is a space translation, from r ⊥ I to R ⊥ , by an effective amount k −1 G ⊥ , proportional to λ. 2. Then, the beam undergoes a boost, accounted for by the action of the momentum operatorp ⊥ I, while it is freely propagating along the z-direction.
The above two-step prescription allows us to determine in a relatively simple fashion the evolution (along z) of Gauss-Maxwell beams (or any linear combination of them) once the input amplitude, F(r ⊥ ; 0), is known. Monitoring this evolution with rays, which are in compliance with the paraxial form of Maxwell's equations, can be done now with the aid of the time-averaged Poynting vector, since where E (r, t) and H(r, t) denote, respectively, the electric and magnetic monochromatic vector fields, solutions to the Maxwell equations, and T is the average over the period of the radiation. Notice that this averaging allows us to also recast the time-averaged Poynting vector just in terms of the time-independent amplitudes E(r) and H(r), which arise from the generalization of E (r, t) and H(r, t), respectively, to the complex domain [2], but that here are directly determined from the output beam F(r ⊥ ; z), at a distance z from the input (transverse) plane. In this case, the guidance equation is defined as where S ⊥ and S z are, respectively, the transverse and longitudinal components of the time-averaged Poynting vector. As before, Equation (46) corresponds to a phase velocity, with its validity being determined by the fact that we are dealing with vacuum, where the phase velocity and the components of the time-averaged Poynting vector are proportional (in general, for nondispersive media, phase velocity and group velocity point in the same direction [2]).

Linearly Polarized Gauss-Maxwell Beams
As seen above, the ray r ⊥,z = (x z , y z ) described by Equations (23) and (24) allows to understand the distribution of the electromagnetic energy along the z axis when the latter is specified by only a scalar field. In this regard, such rays have analogous properties to those that we find for quantum wave packets [34,39]. However, they do not contain any vector-type information, because they have not been obtained from the paraxial form of Maxwell's equations, but from the paraxial approximation applied to Helmholtz's equation. To this end, first we need to determine the electric and magnetic field components associated with the evolution of the input Gaussian amplitude (9), which is done with the aid of the transformation relation (32), and then the rays are determined from Equation (46). In this latter regard, it is interesting to compare the outcome from this equation for a Gauss-Maxwell beam with the rays described by Equations (23) and (24) for a bare Gaussian beam.
For simplicity in the analysis, here we are going to consider the case of horizontal polarization (along the x axis; vertical polarization is taken along the y-axis), leaving the case of arbitrary polarization for the next section. Accordingly, the input electromagnetic vector field, F(r ⊥ ; 0), reads as where we have made use of the relation between amplitudes H 0 = 0 /µ 0 E 0 . Following the transformation relation (32), and taking into account the functional form displayed by the output Gaussian amplitude (12), the action of the operatorÛ on the scala field ψ, generates the matrix operator where the argument of (12) has been replaced by the first argument ofÛ (the second argument ofÛ is the generator of the transformation from z = 0 to finite z), with ∆r ⊥ = r ⊥ − r ⊥,0 . For convenience, we have defined with From now on, because it is more compact, we shall consider the last expression for q z in the derivations. Before further proceeding, it is convenient to simplify the argument of the exponential in (48), just to determine how the matrices involved will act on the column vector (47). To this end, making use of the properties satisfied by the G-matrices, we find with ∆r ⊥ = ∆r ⊥ , and where the property (36) has been used. Now, we consider a Taylor series expansion of the exponential matrix in (48), where, again for convenience, the term ∆r ⊥ · G ⊥ is going to be explicitly separated as ∆xG x + ∆yG y . Thus, we have According to the property (37), the only surviving terms in the second sum of (52) are those with m ≤ 2, hence such an expression can be further simplified to which can be recast in a more compact form as where the prefactor corresponds to the Gaussian function describing the beam (12), while the matrix terms provide us with the correct expression for the electric and magnetic fields, according to the rule (32). After evaluating the combined action of the G-matrices and their products on the column vector (47), we obtain the expression for the output vector field F at a distance z, or, in a more conventional fashion, in terms of the separate electric and magnetic components, as which include up to second-order corrections in the transverse coordinates, withx,ŷ andẑ denoting the unit vectors along the three Cartesian directions. Next, we are going to analyze the physical implications carried by these contributions.
Let us firstly start by neglecting the second-order contributions. Accordingly, the field vector (56) reads as At this level of approximation, we already find the first-order corrections to the paraxial approximation that lead to usual Gaussian beams, with the electric and magnetic fields (57) and (58) being respectively (notice that the π/2 phase factor is not an extra factor, but it arises from the expression for q z introduced above, in Equation (50)). As it can be noticed, the two fields do not remain in-plane (i.e., contained within the same constant z plane), as one would expect from a usual scalar field guess, but they contain a small out-of-plane component, along the longitudinal (z) direction. From these expressions, it is now possible to obtain some information about the polarization state of the electromagnetic field in the near and far fields, that is, in terms of the value of z compared to that of z R . Without loss of generality, in this regard, let us consider the electric component (same holds for the magnetic one). The two cases of interest worth stressing are: • If z z R , then σ z ≈ σ 0 and ϕ z ≈ z/z R π/2. So, in the near field regime, we have and the field is elliptically polarized, with the polarization plane being the XZ-plane (perpendicular to the y transverse coordinate), since As it can be noticed, only at the center of the beam there is horizontal polarization (parallel to the x transverse coordinate); anywhere else, the z-component of the polarization increases linearly with the distance with respect to the center of the beam. Since this a regime where the expansion of the beam is still negligible, this polarization effect should be particularly relevant for distances (from the center) of the order of the beam waist, w 0 , where the fast decrease of the intensity might complicate its detection (particularly, if we consider typical values used in diffraction experiments).

•
On the other hand, for z z R , in the far field regime, we have σ z ≈ zσ 0 /z R and ϕ z ≈ π/2. Accordingly, the electric field component reads as that is, the field is still contained within the XZ-plane, but the elliptical polarization state has degenerated into linear, since which becomes negligible very quickly as z increases, recovering the linear polarization along the x transverse coordinate at any distance from the center of the beam. Therefore, contrary to the case of a monochromatic plane wave, transversality is not ensured unless we consider the limit case σ 0 → ∞, which is precisely the case near the maximum of the Gaussian, where the wavefront could be approximated by a nearly plane wave, as can be seen from Equations (60) and (61), when the z-component contribution is neglected. This model was considered in a previous work [8] to analyze the problem of interference, which is revisited and extended in next section for the case of interference between two Gauss-Maxwell beams.
Once we have the correct paraxial vector field solutions to Maxwell'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 ∆r ⊥ , and compute the time-averaged Poynting vector (45). Thus, neglecting the terms containing (∆x) 2 , (∆y) 2 and ∆x∆y in (60) and (61), we obtain where The ray equation that follows from (66) is where the left-hand side term has been recast in terms of ∆r ⊥ instead of r ⊥ , 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 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, ∆r ⊥ , up to third order, we obtain (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, λ ∼ 10 −6 , while σ 0 , σ z and ∆r ⊥ are the three of the order of ∼ 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.

Arbitrarily Polarized Gauss-Maxwell Beams
In the previous section, for simplicity, we have considered a simple case of linear polarization at the input plane. As seen, the polarization of this beam changes with z. Now we are going to extend the approach to a general input polarization state with the polarization plane perpendicular to the longitudinal direction. Such a state can be described by the superposition with θ ∈ (0, π) and φ ∈ (0, 2π) defined on the Poincaré sphere [40], and where the vector states |H (|P for θ = 0) and |V (|P for θ = π) denote, respectively, horizontal and vertical polarization with respect to the x-axis. For simplicity, from now on, we shall consider the coefficients α = cos(θ/2)e −iδ/2 and β = sin(θ/2)e iφ+iδ/2 , which also include any possible relative phase shift imprinted on the state by some external action (e.g., as a result of a weak measurement). As it can be noticed, these coefficients satisfy the relation |α| 2 + |β| 2 = 1 (unit radius on the Poincaré sphere).
To start with, let us consider the input electric and magnetic field vectors which are expressed in terms of an arbitrary polarization vector field specified by the coefficients α and β introduced above. The corresponding six-dimensional electromagnetic vector field at the input plane z = 0 is which reduces to (47) when α = 1. Taking into account how the different G-matrices and their products operate over the column vector (73), the expression for the propagated polarized vector field becomes As before, from this vector we readily obtain the electric and magnetic vector fields along the z-axis: As seen above, dealing with these full expressions does not lead to any important difference with respect to retaining only up to the first order in the transverse displacement, ∆r ⊥ , and, in the present case, may lead to a rather complex expression for the time-averaged Poynting vector (which does not imply any significant advantage both at the conceptual level and at the methodological one). Therefore, let us thus keep terms at the lowest level in ∆r ⊥ , which reduces the six-component electromagnetic vector field to and the electric and magnetic fields to In order to determine the flux of associated rays, we first compute the time-averaged Poynting vector from the fields (78) and (79), which reads as From (80), we obtain the following set of ray equations which are coupled by virtue of the polarization factor sin θ sin(φ + δ). As it can readily be seen, if the polarization state of the input beam is linear, regardless of its vibration direction (φ = 0 or π), and has not received any extra kick (δ = 0), then the coupling disappears and we recover the same situation that we had in the previous section for the horizontally polarized input Gaussian beam. Otherwise, if the polarization state is not either horizontal or vertical, the second term in the polarization factor provides an extra spin component that makes the rays to get out the plane where they are initially contained, namely the XZ-plane or the YZ-plane. Accordingly, we have a clear description of the energy flux without appealing to a standard picture of rotating electric field vectors, where the electromagnetic energy streamlines already display such a rotation following the spin imprinted by the polarization vector. For instance, in the case of circularly right-handed polarized light (θ = π/2 and φ = −π/2, with δ = 0), described by the (non-normalized) polarization vector which is the case for which the polarization factor is maximal, we have which can be recast as with R ϕ z being the rotation matrix that is, the position vector ∆r ⊥ is acted by a rotation ϕ z − π/2.

Young-Type Interference with Gauss-Maxwell Beams
Now we are going to consider the case of a coherent superposition of two linearly polarized Gaussian beams at the input plane z = 0. This is the case, for instance, of the two diffracted beams produced by two slits (Young-type diffraction), which arise from the incidence of light of a conventional laser pointer linearly polarized. As mentioned above, the standard Gaussian beam description can be used to describe the eventual spatial distribution of radiation, but it does not provide us with any clue on polarization, in case we need it. That is, we would have a two-beam coherent superposition with denoting the diffracted Gaussian beams, with i = 1, 2. For simplicity, although without loss of generality, we are going to consider that both beams are identical (σ 0,i = σ 0 ), they are symmetrically placed with respect to x = 0 (x 0,1 = −x 0,2 = x c ), and aligned along the x-axis (y 0,1 = y 0,2 = 0). As for the prefactor N in (90), it is just a z-independent global norm factor that amounts to approximately 1/ √ 2 when the overlap of the two input beams is nearly zero, as it can readily be seen from its functional form, Without any need for considering the effect of polarization (which is the same to say that both beams are equally polarized, as we have already seen above), the intensity of the radiation described by the scalar superposition field (89), except for constant factors, is going to be proportional to the density profile along the x-direction, where Note that we have denoted (92) with the term density, because of its direct relation to the time-averaged electromagnetic energy density (except for a proportionality constant) [2]. Furthermore, it can readily be seen that, in the limit z z R , where large x values are involved (x x c ), the density (92) can be approximated to with The density (94) describes a Young-type fringe pattern (along the x-direction) modulated by a Gaussian envelop, arising from the Gauss shape of the input beams, and with maxima at positions where d = 2x c . Actually, if we assume that x n /z = sin θ n , with θ n being the angular position of each maximum with respect to the coordinate origin at the input plane, the above condition becomes that is, the usual two-slit interference condition [2]. Let us now consider the case of polarized input beams, particularly, the case of horizontal polarization (in next section, we shall introduce the case of arbitrarily polarized superimposed beams). Following the procedure described in Section 3.2, we find that the electric and magnetic fields (within the paraxial approximation) associated with each partial wave (90) are where, again, we have assumed identical beams (E 0,1 = E 0,2 = E 0 and H 0,1 = H 0,2 = H 0 ). The total electric and magnetic fields that result from the superposition of both beams read as where ∆x i = x − x c,i and ∆ψ = ψ 1 − ψ 2 . From these expressions, we now compute the time-averaged Poynting vector, which reads as As before, the ray equation is now obtained from the components of (103), since flow of the beam at each point is described by the direction of the time-averaged Poynting vector right on that point, which is also the direction along which the wave vector k points [2,41]. Accordingly, we obtain which involves the factor This factor leads to a reduction of the fringe visibility as the polarization states |P 1 and |P 2 become more different, with the maximum loss, with P 2 |P 1 , when these states are orthogonal (e.g., horizontal vs. vertical polarization, or left-handed vs. right-handed circular polarization), in agreement with the aforementioned Arago-Fresnel laws.
By inspecting (108), we readily notice that, again, the fringe distribution is along the x coordinate. In the case of the horizontally polarized beams discussed above, we could get an idea of why this was the behavior: leaving aside the scalar fields themselves, the total electric and magnetic fields explicitly depended on only one transverse coordinate, namely the electric field on the x-coordinate and the magnetic field on the y-coordinate, as it can be seen in Equations (100) and (101), respectively. This is in clear contrast with the fields described by Equations (A12) and (A13), which both depend on both coordinates (in their z-component). That is, not only coordinates and polarization coefficients are mixed up, but there is also a coupling between different field components, which is going to play an interesting role from the point of view of the topology displayed by the corresponding rays. Such a mixture of coordinates becomes more apparent from the time-averaged Poynting vector or the corresponding ray equations along the x and y directions in terms of the longitudinal coordinate z (see expressions in Appendix B). The latter equations reduce to Equations (103) and (104), respectively, in the particular case α 1 = α 2 = 1 and β 1 = β 2 = 0 (i.e., two input beams horizontally polarized at the input plane z = 0). However, unlike the case described by Equations (103) and (104), now the generality of the current equations makes difficult to get a clear clue on what is going on with the dynamics they describe. Nonetheless, still we notice that there is a correlation between transverse coordinates (mediated, in turn, by the coupling with the polarization state) that may induce out-of-plane dynamics (i.e., dependence of one transverse component of the rays on the other, and vice versa), which are not present in the case of either horizontal or vertical polarized input beams [8].

Final Remarks
In this work we have explored a reliable methodology to tackle the study and analysis of experiments performed with Gaussian beams from a ray-based perspective. We could have just remain at the level of standard Gaussian beams, which are acceptable paraxial solutions to Helmholt'z equation and, as it has been shown, they admit a simple treatment in terms of such rays. However, despite its convenience, this approach is limited by the fact that it does not include the description of the polarization state of light, which requires a vector treatment. With such a purpose, we have considered the theoretical framework provided by the approach developed by Mukunda et al. [28] of Gauss-Maxwell beams, which are paraxial vector solutions to Maxwell's equations. Accordingly, we have a formulation that nicely combines the field distribution in coordinates at the same time that accounts for the polarization state, optimal whenever polarization plays a major role, as it is the case, for instance, in the experiment on weak measurements reported by Kocsis et al. [17].
Here, in particular, because light undergoes fast oscillations with time, we have considered a time-averaged approach, which has allowed us to determine rays analogous to those provided by geometrical optics. However, contrary to geometrical rays, the advantage of the rays here introduced is that they follow the evolution of the electromagnetic field and, therefore, provide us with an accurate description of diffraction and interference phenomena. This is achieved by relating them with the components of the time-averaged Poynting vector, which describe how electromagnetic energy distribute spatially. Of course, because we have considered paraxial conditions, we have been able to relate this distribution with the longitudinal component of such a vector instead of assuming a full time-dependent picture-which would be more accurate as well as more complicated, without, however, adding any new physics. Nonetheless, it is worth mentioning that the ray-based picture here provided is able to account for results like those reported in Reference [17] without further appealing to quantum mechanical notions, since the basic ingredients are already contained in standard (classical) electromagnetism. In other words, although Maxwell's equations say nothing about probabilities, still they are useful to reproduce behaviors typically associated with quantum mechanics. Actually, without getting too deeper into the question of whether photons have or have not a wave function, we note that the above approach is somehow in compliance with Bohm's view (see p. 98 in Reference [45]) that, in absence of currents, like the quantum probability density, electromagnetic energy (light) "acts like a fluid, which flows continuously without loss or gain from one point to another". Accordingly, the rays here defined just play the role of the corresponding streamlines that allow us to understand how the energy spatially transfers from one place to another in presence of diffraction, interference and polarization, the three typical traits of electromagnetic optics. Funding: This research and the APC have been both funded by the Spanish Agencia Estatal de Investigación (AEI) and the European Regional Development Fund (ERDF) grant number FIS2016-76110-P.

Conflicts of Interest:
Authors declare no conflict of interest.

Appendix A. General Solution to the Paraxial Equation (5)
The integral (11) is the general solution to the paraxial Equation (5), as it was mentioned in Section 3. To prove this, let us consider without loss of generality the one-dimensional case, and then we shall proceed with the generalization to two dimensions. In such a case, Equation (5)  The solution to ψ(x; z) can be expressed as a linear combination of plane waves, as ψ(x; z) = 1 √ 2π e iκxψ (κ; z)dκ.
The integral over κ can be easily done to yield which accounts for the transverse propagation of the amplitude ψ as a function of the (longitudinal) z coordinate. The generalization of the integral (A8) to two dimensions (the case considered in Section 3) is straightforward. Given there is no correlation between the x and y coordinates, the full solution is just the direct product of (A8) evaluated for x and y, respectively: ψ(r ⊥ ; z) = 1 iλz e ik r ⊥ −r ⊥ 2 /2z ψ(r ⊥ ; 0)dr ⊥ .