{'fraction_non_alphanumeric': 0.059678813759136515, 'fraction_numerical': 0.014199261219250216, 'mean_word_length': 4.224031750376351, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Extending geometrical optics: A Lagrangian theory for vector waves

Even when neglecting diffraction effects, the well-known equations of geometrical optics (GO) are not entirely accurate. Traditional GO treats wave rays as classical particles, which are completely described by their coordinates and momenta, but vector-wave rays have another degree of freedom, namely, their polarization. The polarization degree of freedom manifests itself as an effective (classical)"wave spin"that can be assigned to rays and can affect the wave dynamics accordingly. A well-known manifestation of polarization dynamics is mode conversion, which is the linear exchange of quanta between different wave modes and can be interpreted as a rotation of the wave spin. Another, less-known polarization effect is the polarization-driven bending of ray trajectories. This work presents an extension and reformulation of GO as a first-principle Lagrangian theory, whose effective-gauge Hamiltonian governs the aforementioned polarization phenomena simultaneously. As an example, the theory is applied to describe the polarization-driven divergence of right-hand and left-hand circularly polarized electromagnetic waves in weakly magnetized plasma.


A. Motivation
Geometrical optics (GO) is a reduced model of wave dynamics [1,2] that is widely used in many contexts ranging from quantum dynamics to electromagnetic (EM), acoustic, and gravitational phenomena [3][4][5]. Mathematically, GO is an asymptotic theory with respect to a small parameter that is a ratio of the wave relevant characteristic period (temporal or spatial) to the inhomogeneity scale of the underlying medium. Practical applications of GO are traditionally restricted to the lowest-order theory, where each wave is basically approximated with a local eigenmode of the underlying medium at each given spacetime location. Then, the wave dynamics is entirely determined by a single branch of the local dispersion relation. However, this approximation is not entirely accurate, even when diffraction is neglected. If a dispersion relation has more than one branch, i.e., a vector wave with more than one polarization at a given location, then the interaction between these branches can give rise to important polarization effects that are missed in the traditional lowest-order GO.
One interesting manifestation of such polarization effects is the polarization-driven bending of ray trajectories. At the present moment, it is known primarily in two contexts. One is quantum mechanics, where polarization effects manifest as the Berry phase [6] and the associated Stern-Gerlach force experienced by vector particles, i.e., quantum particles with spin. Another one is optics, where a related effect has been known as the Hall effect of light; namely, even in an isotropic dielectric, rays propagate somewhat differently depending on polarization if the dielectric is inhomogeneous (see, c.f. Refs. [7][8][9][10][11]). But the same effect can also be anticipated for waves in plasmas, e.g., radiofrequency (RF) waves in tokamaks. In fact, since for RF waves in laboratory plasma is typ-ically larger than that for quantum and optical waves, the polarization-driven bending of ray trajectories in this case can be more important and perhaps should be taken into account in practical ray-tracing simulations. However, ad hoc theories of polarization effects available from optics are inapplicable to plasma waves, which have more complicated dispersion and thus require more fundamental approaches. Thus, a different theory is needed that would allow the calculation of the polarization-bending of the ray trajectories for plasma waves and, also more broadly, waves in general linear media.
Relevant work was done in Refs. [12,13], where a systematic procedure was proposed to asymptotically diagonalize the dispersion operator for linear vector waves. Polarization effects emerge as O( ) corrections to the GO dispersion relation. However, this approach excludes mode conversion, i.e., the linear exchange of quanta between different branches of the local dispersion relation. Since the group velocities of the different branches eventually separate, mode conversion is typically followed by ray splitting and, in this particular context, was studied extensively (see, e.g., Refs. [14][15][16][17][18][19][20]). However, these works considered wave modes that are resonant in small, localized regions of phase space. Hence, the nonadiabatic dynamics was formulated as an asymptotic scattering problem between two wave modes, so the polarizationdriven bending of ray trajectories was not included.
The main message of this work is that mode conversion and the polarization-driven bending of ray trajectories are two sides of the same coin and can be considered simultaneously within a unified theory. The first general theory that captures them both simultaneously was proposed in Ref. [21]. This theory was successfully benchmarked [22] against previous theories describing the Hall effect of light [7][8][9][10][11]. However, the formulation in Ref. [21] is still limited since it requires that the wave equation be brought to a certain (multisymplectic) form resembling the Dirac equation. Although any nondissipative vector wave allows for such representation in prin-ciple [21,23], casting the wave dynamics into the specific framework adopted in Ref. [21] can be complicated. Thus, practical applications require a more flexible formulation that do not rely on this specific framework.
Here we propose such a theory. In addition to generalizing the results of Ref. [21], we also introduce, in a unified context and an instructive manner, some of the related advances that were made recently in Refs. [24][25][26]. It is expected that the comprehensive analysis presented in this work will facilitate future practical implementations of the proposed theory, particularly in improving ray-tracing simulations.

B. Outline
We consider general linear nondissipative waves determined by some Hermitian dispersion operator. Using the Feynman reparameterization and the Weyl calculus, we obtain a reduced Lagrangian for such waves. In contrast with the traditional GO Lagrangian, which has an accuracy of O( 0 ) in the GO parameter , our Lagrangian is O( 1 )-accurate, so it captures polarization effects. As an example, we apply the formulation to study polarization effects on the propagation of EM waves propagating in weakly magnetized plasma. (The case of strongly magnetized plasma will be discussed in a separate paper.) The advantages of our theory are as follows. (i) The theory is derived in a variational form, so the resulting equations are manifestly conservative. (ii) Through the use of the Feymann reparameterization, we can obtain the dynamics of continuous waves and of their rays directly from a variational principle. (iii) The theory assumes no specific wave equation, so quantum spin effects and classical polarization effects can be studied on the same footing. (iv) Moreover, a related formalism [21] is applicable to develop new reduced theories for relativistic spinning particles [21,24]. This paper is organized as follows. In Sec. II, the basic notation is defined. In Sec. III, the variational formalism used to describe vector waves is presented. In Sec. IV, a general procedure is proposed to block-diagonalize the wave dispersion operator. In Sec. V, we reparametrize the wave action to facilitate asymptotic analysis. In Sec. VI, the leading order GO approximation is discussed. In Sec. VII, the more accurate model that includes polarization effects is discussed. In Sec. VIII, the theory is applied to describe polarization effects on the propagation of EM waves in weakly magnetized plasma. In Sec. IX, our main results are summarized. Finally, Appendix A presents a brief introduction to the Weyl symbol calculus.

II. NOTATION
The following notation is used throughout the paper. The symbol " . =" denotes definitions, "c. c." denotes "complex conjugate," and "h. c." denotes "Her-mitian conjugate." The identity N × N matrix is denoted by I N . The Minkowski metric is adopted with signature (+, −, −, −). Greek indices span from 0 to 3 and refer to spacetime coordinates x µ = (x 0 , x) with x 0 corresponding to the time variable t. Also, partial derivatives on spacetime will be denoted by ∂ x , where the individual components are ∂ µ . = ∂/∂x µ = (∂ t , ∇) and d 4 x . = dt d 3 x. Latin indices span from 1 to 3 and denote the spatial variables, i.e., x = (x 1 , x 2 , x 3 ) and ∂ i . = ∂/∂x i . Summation over repeated indexes is assumed. In particular, for arbitrary four-vectors a and b, In Euler-Lagrange equations (ELEs), the notation "δa :" means that the corresponding equation is obtained by extremizing the action integral with respect to a.

A. Wave action principle
The dynamics of any nondissipative linear wave can be described by the least action principle, δS = 0, where the real action S is bilinear in the wave field [23]. We represent a wave field, either classical or quantum, as a complex-valued vector Ψ(x). We allow this vector field to have an arbitrary number of componentsN . In the absence of parametric resonances [27], the action can be written in the form [28] S .
where D is aN ×N Hermitian matrix kernel [D(x, x ) = D † (x , x)] that describes the underlying medium. Varying S with respect to Ψ † leads to Similarly, varying with respect to Ψ gives the equation adjoint to Eq. (2), which we do not need to discuss. It is convenient to describe the wave Ψ(x) as an abstract vector |Ψ in the Hilbert space of wave states with inner product [23,29] Υ|Ψ In this representation, Ψ(x) = x|Ψ , where |x are the eigenstates of the coordinate operatorx such that We also introduce the momentum (wavevector) operatorp such that x|p µ |x = i∂[δ 4 (x − x )]/∂x µ in the x-representation [30]. Thus, the action (1) can be rewritten as whereD is the Hermitian dispersion operator such that D(x, x ) = x|D|x . Treating Ψ| and |Ψ as independent variables [23] and varying the action (4) gives δ Ψ| :D |Ψ = 0 (plus the conjugate equation), which is the generalized vector form of Eq. (2). Specifically, Eq. (2) is obtained by projecting Eq. (5) with x| and using the fact that the operator d 4 x |x x| =1 is an identity operator.

B. Extended wave function
As shown in Refs. [21,23], reduced models of wave propagation are convenient to develop when the action is of the symplectic form; namely, wherep 0 = i∂ t (in the x-representation) and "the wave Hamiltonian"Ĥ = H(x,p) is some Hermitian operator that is local in time, i.e., commutes witht. (For extended discussions, see Refs. [23,31].) In order to cast the general action (4) into the symplectic form (6), let us perform the so-called Feynman reparameterization [32,33] that lifts the wave dynamics governed by Eq. (4) from R 4 to R 5 . Specifically, we let the wave field depend on spacetime and on some parameter τ so that Ψ(τ, x) = x|Ψ(τ ) . Note that |Ψ(τ ) belongs to the same Hilbert space defined in Sec. III A. Thus, the inner product remains the same; i.e., Υ(τ )|Ψ(τ ) = d 4 x Υ † (τ , x)Ψ(τ, x). We consider the following "extended" action: where L .
and ∂ τ Ψ(τ, x) = x|∂ τ Ψ(τ ) . Note that the Lagrangian L is local in the parameter τ ; i.e., the abstract vectors are all evaluated at τ . From hereon, all fields will be evaluated at τ , and we will avoid mentioning the dependence of |Ψ on τ explicitly. The ELE corresponding to the action S X is given by Note that Eq. (9) can be interpreted as a vector Schrödinger equation in the extended variable space, whereD acts as the Hamiltonian operator. The dynamics of the original system described by Eq. (5) is a special case of the dynamics governed by Eq. (9), which corresponds to a steady state with respect to the parameter τ ; i.e., ∂ τ Ψ = 0. The advantage of the representation (7) is that the action has the manifestly symplectic form, so we can proceed as follows.

A. Variable transformation
We introduce a unitary τ -independent transformation Q that maps |Ψ to someN -dimensional abstract vector |ψ yet to be defined: Inserting Eq. (10) into Eqs. (8) leads to In what follows, we seek to con-structQ such that the operatorD eff is expressed in a block-diagonal form. The procedure used is identical to that given in Refs. [12,13]. However, in order to account for resonant-mode coupling,D eff will be made only blockdiagonal, instead of fully diagonal as in Refs. [12,13].

B. Weyl representation
Let us consider Eq. (11b) in the Weyl representation. (Readers who are not familiar with the Weyl calculus are encouraged to read Appendix A before continuing further.) In this representation, L D is written as [28] L D = Tr d 4 where 'Tr' represents the matrix trace. The Wigner tensor W (τ, x, p) corresponding to |ψ is defined as and D eff (x, p) is the Weyl symbol [Eq. (A1)] corresponding to the operatorD eff . It can be written explicitly as where ' ' is the Moyal product [Eq. (A6)] and D(x, p), Q(x, p), and [Q † ](x, p) are the Weyl symbols corresponding toD,Q, andQ † , respectively. Also, the Weyl representation of the unitary condition,Q †Q =ÎN , is which will be used below.

C. Eigenmode representation
Let us assume that the symbols D eff and Q can be expanded in powers of the GO parameter where ω and |k| are understood as the characteristic wave frequency and wave number, respectively. Also, T and are the characteristic time and length scales of the background medium, correspondingly. Hence, we write where (Λ, U, Q 0 , Q 1 ) areN ×N matrices of order unity.
To the lowest-order in , the Moyal products in Eqs. (14) and (15) reduce to ordinary products, so By properties of the Weyl transformation, the fact thatD is a Hermitian operator ensures that D(x, p) is a Hermitian matrix. Hence, D(x, p) hasN orthonormal eigenvectors e q (x, p), which correspond to some real eigenvalues λ (q) (x, p). Let us construct Q 0 (x, p) out of these eigenvectors so that where the individual e q form the columns of . Hence, the matrix Λ has the following diagonal form: To the next order in , Eq. (15) reads as follows: Here we assumed that term that involves the Poisson bracket {Q † 0 , Q 0 }, which arises from the expansion of the Moyal star product [Eq. (A9)], is of the first order in . Following Ref. [12], we let In order to determine G(x, p), we write Eq. (14) to the first order in . Introducing the bracket and noting that DQ 0 = Q 0 Λ, we obtain [12] U where δU (x, p) is aN ×N matrix given by , as done in Refs. [12,13]. However, this method is appli- , which is in violation of the assumed ordering in Eq. (17b). Hence, instead of diagonalizing U , we propose to only blockdiagonalize U as follows. When |λ (m) − λ (n) | O(1), we choose the off-components of G m n so that U m n = 0. (We call such modes nonresonant.) When |λ (m) −λ (n) | O( ), we let G m n = 0. (We call such modes resonant.) By following this prescription and permutating the ma-trix rows, we obtain U in the following form: where U j (x, p) are n j × n j Hermitian matrices and J is the total number of blocks, so J j=1 n j =N . Note that, in the particular case where only nonresonant modes are present, U (x, p) is diagonal, and one recovers the results obtained in Refs. [12,13].
Since the matrix D eff ≈ Λ+ U is made block-diagonal, the Lagrangian (12) is unaffected by the matrix elements W m n with indices (m, n) such that U m n = 0. Thus, without loss of generality, we can write Now that blocks of mutually nonresonant modes are decoupled, let us focus on the dynamics of modes within a single block of some size N . Hence, the block index will be dropped, and we adopt Here ψ is a complex-valued function with N components, and one obtains which is a N × N Hermitian matrix. Furthermore, it is convenient to split [[D eff ]] as follows: where we used the bracket introduced in Eq. (24) and the subscript 'A' denotes "anti-Hermitian part;" i.e., for any matrix M , then The expression in Eq. (34) can also be written more explicitly as Examples of physical systems, where these simplified formulas are applicable, include spin-1/2 particles [21,24] and EM waves propagating in isotropic dielectrics [22].

B. Parameterization of the action
In order to derive the corresponding ELEs, let us adopt the following parameterization: Here θ(τ, x) is a real variable that serves as the rapid phase common for all N modes (remember that all modes within the block of interest are approximately resonant to each other). Also, I(τ, x) is a real function, and z(τ, x) is a N -dimensional complex unit vector (z † z = 1), whose components describe the amount of quanta in the corresponding modes. (Since we parameterize the N -dimensional complex vector ψ by the N -dimensional complex vector z plus two independent real functions θ and I, not all components of z are truly independent. For an extended discussion, see Ref. [21].) After substituting the ansatz (36) into Eq. (29a), the Lagrangian L τ is given by (Here we formally introduce to denote that z is a slowlyvarying quantity; however, this ordering parameter will be removed later.) Now, we calculate the Wigner tensor (30). Substituting Eq. (36) into Eq. (30), we obtain where we introduced the four-wavevector k µ (τ, x) .
, which is considered a slow function. [Accordingly, the contravariant representation is k µ (x, τ ) = (ω, k).] Inserting Eq. (38) into Eq. (29b) and integrating over the momentum coordinate, we obtain where we integrated by parts and used z † z = 1. Here is the zeroth-order (in ) group velocity of the wave. We then introduce the convective derivative d dτ Summing Eqs. (37) and (39), we obtain the action S = dτ L + O( 2 ), where the Lagrangian is given by Equation ( In what follows, we discuss the consequences of this theory and provide an example, where we apply the theory to study polarization effects on EM waves in weakly magnetized plasmas.

A. Continuous wave model
To lowest order in , the Lagrangian (42) can be approximated simply with which one may interpret as a Hayes-type representation [34] of the GO wave Lagrangian in the extended (τ, x) space. This Lagrangian is parameterized by just two functions, the rapid phase θ and the total action density I. Thus, varying the action S GO = dτ L GO , we obtain the following ELEs: where v µ (τ, x) is the GO four-group-velocity (40). As mentioned in Sec. III, the dynamics of the physical wave propagating in spacetime is obtained by adopting ∂ τ Ψ = 0, which also corresponds to ∂ τ I = ∂ τ θ = 0. Hence, Eqs. (44) become λ(x, k) = 0.
Equation (45a) is the action conservation theorem, or the photon conservation theorem. Equation (45b) is the local dispersion relation. For an in-depth discussion of these equations, see, e.g., Refs. [1,4].

B. Point-particle model
The ray equations corresponding to the above field equations can be obtained as the point-particle limit. In this limit, I can be approximated with a delta function Here I 0 denotes the total action, which is conserved according to Eq. (45a). The value of I 0 is not essential below so we adopt I 0 = 1 for brevity. In this representation, the wave packet is located at the position X(τ ) in space-time, and the independent parameter is τ . [This means that at a given τ , the wave packet is located at the spatial point X(τ ) at time t(τ ).] When inserting Eq. (46) into Eq. (43), the first term in the action gives the following: where P µ (τ ) . = −∂ µ θ(τ, X(τ )). Similarly, Thus, the point-particle action is expressed as This is a covariant action, where X(τ ) and P (τ ) serve as canonical coordinates and canonical momenta, respectively. Treating X and P as independent variables leads to ELEs matching Hamilton's covariant equations These are the commonly known ray equations; for instance, see Ref. [1]. They can also be written as Note that the first term in the integrand in Eq. (49) represents the symplectic part of the canonical phasespace Lagrangian, and the second term represents the Hamiltonian part. Since the Hamiltonian part λ(X, P ) does not depend explicitly on τ , then dλ(X, P )/dτ = 0 along the ray trajectories. Thus, the ray dynamics lies on the dispersion manifold defined by λ(X, P ) = 0. (51) As a reminder, λ(x, p) is defined as the average eigenvalue of the resonant block, i.e., λ .
The GO action (49) is only accurate to lowest order in ; hence, one can approximate λ(x, p) λ (n) (x, p), where λ (n) is any particular resonant eigenvalue. This occurs because the resonant eigenvalues differ by O( ) and because the polarization coupling is also O( ).

VII. EXTENDED GEOMETRICAL OPTICS
In this section, we explore the polarization effects determined by the Lagrangian (42). For the sake of conciseness, we only discuss the point-particle ray dynamics. For an overview of the continuous-wave model, see Ref. [21].

A. Point-particle model
The ray equations with polarization effects included can be obtained as a point-particle limit of the Lagrangian (42). As in Sec. VI B, we approximate the wave packet to a single point in spacetime [Eq. (46)]. As shown in Refs. [21,25], the Lagrangian (42) can be replaced by a point-particle Lagrangian so the action is where Z(τ ) . = z(τ, X(τ )) is the point-particle polarization vector and we dropped the GO ordering parameter . In the complex representation, Z and Z † are canonical conjugate, and Even though the components of Z are not independent by definition (Sec. V B), it can be shown [21] that treating them as independent in this point-particle model leads to correct results provided that the initial conditions satisfy Eq. (53). Hence, the independent variables in S XGO are (X, P, Z, Z † ), and the corresponding ELEs are δZ : Together with Eqs. (31)-(33), Eqs. (54) form a complete set of equations. The first terms on the right-hand side of Eqs. (54a) and (54b) describe the ray dynamics in the GO limit. The second terms describe the coupling to the mode polarization. Equations (54c) and (54d) describe the wave-polarization dynamics. As in Sec. VI B, the Hamiltonian part of Eq. (52) is constant along the ray trajectories. As before, the ray dynamics lies on the dispersion manifold defined by setting the Hamiltonian part to zero; i.e., λ(X, P ) + Z † U(X, P )Z = 0.

B. Precession of the wave spin
Let us also describe the rotation of Z(τ ) as follows. Since U(X, P ) is a traceless Hermitian N × N matrix, it can be decomposed into a linear combination of N 2 − 1 generators T u of SU(N ), which are traceless Hermitian matrices, with some real coefficients −W u [35]: Then, we introduce the (N 2 − 1)-dimensional vector so that Z † UZ = −S·W. The components of S(τ ) satisfy the following equation: where f abc are structure constants. They are defined via [T a , T b ] = if abc T c so that the structure constants f abc are antisymmetric in all indices [35].
For example, consider the case when only two waves are resonant. Then, N 2 − 1 = 3, T v are the three Pauli matrices divided by two (so |S| 2 = 1/2), and f wuv is the Levi-Civita symbol, so f wvu S v W u = (S × W) w . For a Dirac electron, which is a special case, such S is recognized as the spin vector undergoing the well known precession equation, d τ S = S × W [21]. In optics, this is an equation for the Stokes vector that was derived earlier to characterize the polarization of transverse EM waves in certain simple media [7,36,37].
Hence, it is convenient to extend this quantum terminology also to N resonant waves. We will call the corresponding (N 2 − 1)-dimensional vector S a generalized "wave-spin" vector and express f wvu S v W u symbolically as (S * W) w , where ' * ' can be viewed as a generalized vector product. Notably, using the concept of spin vector S, one can rewrite Eqs. (54c) and (54d) as follows: which is understood as a generalized precession equation.
In the particular case when S is conserved (we call such waves "pure states"), then Eqs. (54a), (54b), and (59) form a closed set of equations, and λ − S · W serves as an effective scalar Hamiltonian. The dynamics of Z and Z † does not need to be resolved in this case, so one can rewrite S XGO as a functional of (X, P ) alone: An example of the dynamics described by such action will be discussed in Sec. VIII E.
A more general case is when S is close to some eigenvector w of W that corresponds to some nondegenerate eigenvalue Ω w . If Ω w is large enough, then S(τ ) will remain close to w(τ ) and will only experience smallamplitude oscillations. These oscillations can be understood as a generalized zitterbewegung effect [38], and they are transient, i.e., vanish whenẆ becomes zero. In this regime, no mode conversion occurs at τ → ∞. In contrast, if Ω w is not large enough, the change of S governed by Eq. (59) is not necessarily negligible. This corresponds to mode conversion and causes ray splitting at τ → ∞ (see, e.g., Refs. [14][15][16][17][18][19][20]). This is discussed below.
= −i|β| 2 /α, the equation above can be written as a Weber equation whose solutions are the parabolic cylinder functions D n (w). In Refs. [17,39], the matrix connecting the waves entering and exiting the resonance are obtained by analyzing asymptotics of D n (w). Specifically, where where Γ is the Gamma function and η . = β/ √ α. The transmission and conversion coefficients for the wave quanta are, correspondingly, (Also see Ref. [16] for a somewhat different approach leading to the same answer.) This calculation shows that mode conversion, in the way as commonly described in literature [14][15][16][17][18][19][20], is nothing but a manifestation of the wave-spin precession described by Eqs. (54c) and (59). Note that the present point-particle model cannot capture ray-splitting because it introduces only one ray for the whole field. However, this theory does predict the transfer of wave quanta, which is a prerequisite for ray-splitting. For a complete analysis on ray-splitting mode conversion, please refer to Refs. [1,18,19].

VIII. DISCUSSION: WAVES IN WEAKLY MAGNETIZED PLASMAS
A simplified form of the theory above was applied to describe spin-1/2 particles [21,24] and waves in isotropic dielectrics [22]. Here we present another example of its application, namely, EM waves in weakly magnetized cold plasmas. (The case of strongly magnetized plasmas will be discussed in a separate paper.) We assume that the plasma response is determined by particles of just one type, e.g., electrons. The generalization to multicomponent plasma is straightforward to do.

A. Dispersion operator
The linearized equations of motion are [40] Here q, m, n 0 (x), and v(t, x) are the particle charge, mass, unperturbed background density, and flow velocity, respectively. Also, E(t, x) denotes the perturbation electric field, B(t, x) is the perturbation magnetic field, B 0 (x) is the background magnetic field, and c is the speed of light. We introduce a re-scaled velocity field v(t, x) where ω p (x) . = [4πq 2 n 0 (x)/m] 1/2 is the plasma frequency and Ω (x) . = qB 0 (x)/(mc) is the gyrofrequency.
Let us write Eqs. (70) using the abstract Hilbert space notation. Let |v be a state vector representing the velocity field such that v(x) = x|v . Likewise, we introduce |E and |B as the state vectors of E(x) and B(x), respectively. Then, Eqs. (70) can be written as follows: whereω p .
= Ω(x). (As a reminder, p 0 = i∂ t andp = −i∇ are the components of the fourmomentum operator in the x-representation.) Also, α . = α 1 , α 2 , α 3 are 3 × 3 Hermitian matrices [41] These matrices serve as generators for the vector product. Namely, for any two column vectors A and B, one has where the superscript 'T' denotes the matrix transpose. The next step is to construct a dispersion operator for the electric field state |E . Starting from Eq. (71a), we solve for the velocity field in terms of the electric field. Hence, we formally obtain the following: Similarly, we obtain |B = −ic(α ·p)p −1 0 |E from Eq. (71c). Substituting these results into Eq. (71b), we obtain wherê serves as the dispersion operator for |E . (For convenience, we let c = 1.) Since ω p (x) and Ω(x) are independent of time, thenp 0 commutes withω p andΩ, sô D is manifestly Hermitian. The corresponding action (4) for the electric field is S = E|D|E , and the extended action (7) is Note that E is a three-dimensional vector field, soN = 3.

B. EM waves in weakly magnetized plasma
We now follow the procedure given in Secs. IV and V to block-diagonalize the dispersion operator. The Weyl symbol ofD is For the sake of simplicity, we consider the case of a wave propagating in a weakly magnetized plasma. (The general case will be described in a separate paper.) Thus, supposing that the typical wave frequency is much larger than the gyrofrequency (ω ∼ p 0 Ω), we expand the dispersion symbol (78) in powers of Ω/p 0 : where To simplify the following calculation, we assume that D 1 ∼ O(Ω/p 0 ) is comparable in magnitude to the GO parameter , but this is not essential. Hence, we will consider D 1 as a perturbation only. Following Sec. IV C, the next step is to identify the eigenvalues and eigenmodes of the dispersion symbol D 0 (x, p). The corresponding eigenvalues are where p · p = p 2 0 − p 2 . These eigenvalues correspond to the dispersion relations of two transverse EM waves and of longitudinal Langmuir oscillations, respectively. The matrix Q 0 defined in Eq. (20) is given by where e 1 (p) and e 2 (p) are any two orthonormal vectors in the plane normal to e p (p) . = p/|p|. A right-hand convention is adopted such that e 1 × e 2 = e p . One can easily verify that these vectors are indeed eigenvectors of D 0 (x, p). For example, where Eq. (73a) was used. Similar calculations follow for the other two eigenmodes e 2 and e p . We now analyze the dynamics of the transverse EM waves.
From Sec. V, the eigenvalue is λ(x, p) = −p · p + ω 2 p (x), and Ξ(p) = [e 1 (p), e 2 (p) ] is a 3 × 2 matrix. Since Ξ(p) only depends on the spatial momentum coordinate, then the polarization-coupling Hamiltonian U in Eq. (35) is given by where e q is the dual to e q , so e q e r = δ q r . (Specifically, e q is a row vector, whose elements are complex-conjugate of those of e q .) Since ∇λ = ∇ω 2 p , we can also write Eq. (84) in the form where σ y is the y-component of the Pauli matrices and F(p) is a vector with components given by For example, one may choose so that where p ⊥ .
= p x , p y , 0 T ; or, more explicitly, (The specific choice of e 1 and e 2 does not affect the resulting equations within the accuracy of the present theory. For more details, see Sec. VIII G.) Returning to the perturbation caused by the background magnetic field, the projection of the eigenmodes on the matrix D 1 (x, p 0 ) is given by where we used Eq. (73b).

C. Ray dynamics
Now, let us discuss the point-particle ray dynamics. Following Sec. VII, we substitute λ(x, p) = −p·p+ω 2 p (x), Eq. (85), and Eq. (91) into Eq. (52). We then obtain the point-particle action where the polarization-coupling matrix is given by and Z(τ ) is a complex-valued vector with two components that describe the degree of polarization along the vectors e 1 and e 2 . It is normalized such that Z † Z = 1.
In the action (92), the two polarization modes are coupled through the Pauli matrix σ y . However, these modes can be decoupled when using the basis of circularly polarized modes. We introduce the variable transformation where and Γ(τ ) is a new vector with components denoted as Inserting Eq. (94) into the action (92) leads to −P · P + ω 2 p (X) + Σ(X, P )Γ † σ z Γ , (97) where σ z is another Pauli matrix,  104)]. The blue and red lines represent the ray trajectories for the right-hand and left-hand polarized rays, respectively. For simplicity, nonmanetized plasma is considered, so the Faraday effect is absent. The plasma frequency is given by ω 2 p (x) = y 2 + z 2 . The initial location of the ray trajectories is X0 = (0, 1, 0), and the initial momentum is P0 = (5, 0, 1). (The units are arbitrary, since the figure is a general illustration only.) For this simulation, the GO parameter is roughly ∼ 1/ |P0| ∼ 0.2. Due to the radial gradient in the plasma frequency, the wave rays follow helical trajectories along the x axis.
Here Γ ± (τ ) represent the wave quanta belonging to the right-hand and left-hand circularly polarized modes, respectively (as defined from the point of view of the source). Also, Γ is normalized such that Γ † Γ = 1.
Treating X(τ ), P (τ ), Γ(τ ), and Γ † (τ ) as independent variables, we obtain the following ELEs: δΓ : Together with Eq. (93), Eqs. (99) form a complete set of equations. The first terms on the right-hand side of Eqs. (99a) and (99b) describe the ray dynamics in the GO limit. The second terms describe the coupling of the mode polarization and the ray curvature.

D. Restating the Faraday effect
To better understand the polarization equations, let us rewrite Eq. (99c) as an equation in the basis of linearly polarized modes: [This equation could also be obtained if the ray equations were derived directly from the action (92).] Since Σ is a scalar and σ y is constant, this can be readily integrated, where Θ(τ ) . = τ 0 dτ Σ(X(τ ), P (τ )) is the polarization precession angle and Z 0 . = Z(τ = 0). This result can be also be expressed explicitly as follows: It is seen that the polarization of the EM field rotates at the rate Σ(t) in the reference frame defined by the basis vectors (e 1 , e 2 ). The first term in Eq. (93) is identified as the rate of change of the wave Berry phase [6]. (In optics, the rotation of the polarization plane caused by the Berry phase is also known as the Rytov rotation [37,43,44].) The second term in Eq. (93) is identified as the rate of change due to Faraday rotation.

E. Dynamics of pure states
If a ray corresponds to a strictly circular polarization such that σ z Γ = ±Γ, the action (97) can be simplified to S XGO = dτ L ± , where the Lagrangian is given by Here the Lagrangian L ± governs the propagation of righthand and left-hand polarization modes, respectively. The corresponding ELEs are Hence, with the use of the noncanonical coordinates (x, p), the equations of motion no longer depend on the specific choice of F(p); i.e., they are invariant with respect to the choice (88) of vectors e 1 and e 2 . Note that the same equations could be obtained directly from the point-particle limit of Eq. (84), if one substitutes −∇λ =ṗ. For an extended discussion of pure states governed by noncanonical Lagrangians, see Ref. [12].

IX. CONCLUSIONS
Even when neglecting diffraction, the well-known equations of geometrical optics (GO) are not entirely accurate. Traditional GO treats wave rays as classical particles, which are completely described by their position and momentum coordinates. However, vector waves have another degree of freedom, namely, their polarization. Polarization dynamics are manifested in two forms: (i) mode conversion, which is the transfer of wave quanta between resonant eigenmodes and can be understood as the precession of the wave spin, and (ii) polarization-driven bending of ray trajectories, which refers to deviations of the GO ray trajectories arising from first-order corrections to the GO dispersion relation. They are easily understood by drawing parallels with quantum mechanics, where similar effects (yet involving ) are known as spin rotation and spin-orbital coupling.
In this work, we propose a first-principle variational formulation that captures both types of polarizationrelated effects simultaneously. We consider general linear nondissipative waves, whose dynamics are determined by some dispersion operatorD. Using the Feynman reparameterization and the Weyl calculus, we obtain a reduced Lagrangian model for such general waves. In contrast with the traditional GO Lagrangian, which is O( 0 )accurate in the GO parameter , our Lagrangian is O( )accurate. In our procedure, polarization effects are contained in the O( ) corrections to the GO Lagrangian. These corrections may be especially significant for modeling RF waves in laboratory plasmas because such waves can have not-too-small (as opposed, for instance, to quantum particles whose spin effects are typically weak). As an example, we apply the formulation to study the polarization-driven divergence of RF waves propagating in weakly magnetized plasma. Assessing the importance of polarization effects on waves propagating in strongly magnetized plasma will be discussed in a separate paper. Likewise, the method of including dissipation [26] in the above theory will also be described separately. This appendix summarizes our conventions for the Weyl transform. (For more information, see the excellent reviews in Refs. [1,[45][46][47].) The Weyl symbol A(x, p) of any given operatorÂ is defined as where p·s = p 0 s 0 −p·s and the integrals span over R 4 . We shall refer to this description of the operators as a phasespace representation since the symbols (A1) are functions of the eight-dimensional phase space. Conversely, the inverse Weyl transformation is given bŷ (A2) Hence, A(x, x ) = x|Â|x can be expressed as In the following, we outline a number of useful properties of the Weyl transform.
• For any operatorÂ, the trace Tr [Â] . = d 4 x x|Â|x can be expressed as • If A(x, p) is the Weyl symbol ofÂ, then A † (x, p) is the Weyl symbol ofÂ † . As a corollary, the Weyl symbol of a Hermitian operator is real.
Here ' ' refers to the Moyal product, which is given by andL is the Janus operator The arrows indicate the direction in which the derivatives act, and ALB = {A, B} is the canonical Poisson bracket in the eight-dimensional phase space, namely, Provided that ALB is small, one can use the following asymptotic expansion of the Moyal product: • The Moyal product is associative; i.e., • Now we tabulate some Weyl transforms of various operators. (We use a two-sided arrow to show the correspondence with the Weyl transform.) First of all, the Weyl transforms of the identity, position, and momentum operators are given bŷ For any two functions f and g, one has f (x) ⇔ f (x), g(p) ⇔ g(p).