1 Introduction↩︎

The weighted least action principle of geometrical optics for scalar wave equations was introduced in [1] as a generalization of the Fermat principle of least time. It was further shown [2] that the geometrical optics (semiclassical) limit of many dispersive wave equations is related to the optimal transport paradigm [3] proposed by Monge in the late 18th century. The general framework of optima (mass) transport problem is as follows: Consider two distributions, or intensities, \(I_1\) over a domain \(\Omega_1\) on a first plane \(z=0\), and \(I_2\) over a domain \(\Omega_2\) on a second plane \(z=h\), such that \(\int_{\Omega_1} I_1 = \int_{\Omega_2} I_2\). Let \({\boldsymbol{T}}({\boldsymbol{x}})\) be a mapping from \(\Omega_1\) to \(\Omega_2\). We say that \({\boldsymbol{T}}\) transforms \(I_1\) to \(I_2\) if it satisfies \[\label{1} \int_{\Omega_1} I_1({\boldsymbol{x}})\phi({\boldsymbol{T}}({\boldsymbol{x}})) = \int_{\Omega_2} I_2({\boldsymbol{x}}) \phi({\boldsymbol{x}})\tag{1}\] for any function \(\phi\) continuous on \(\Omega_2\). The class of such mappings \({\boldsymbol{T}}\) is denoted \({\cal C}\). Note that if \({\boldsymbol{T}}\in {\cal C}\) is differentiable then (1 ) implies \[I_1({\boldsymbol{x}})= I_2({\boldsymbol{T}}({\boldsymbol{x}})) |J({\boldsymbol{T}})|, \label{p2}\tag{2}\] where \(J\) is the Jacobian of the mapping \({\boldsymbol{T}}\). Let \(C({\boldsymbol{x}},{\boldsymbol{y}})\) a real valued function on \(\Omega_1\times\Omega_2\) be a cost function. The Monge problem is to find a mapping \(\bar{{\boldsymbol{T}}} \in {\cal C}\) that minimizes \[\int_{\Omega_1} I_1({\boldsymbol{x}}) C({\boldsymbol{x}},{\boldsymbol{T}}({\boldsymbol{x}})) \;d{\boldsymbol{x}}\label{p4}\tag{3}\] within the class \({\boldsymbol{T}}\in{\cal C}\). In the optics context, the mapping \(\bar{\boldsymbol{T}}\) is induced by the rays connecting points on \(\Omega_1\) to points on \(\Omega_2\). In [1] it is shown that for a proper selection of the cost function \(C\), the minimization problem above provides the full ray mapping between two distributions of light. The specific cost function \(C\) is determined by the Hamiltonian structure of the wave problem. For example, the optimal mapping \(\bar{\boldsymbol{T}}\) of the Monge problem for the cost function \(C({\boldsymbol{x}},{\boldsymbol{T}}({\boldsymbol{x}}))=n\sqrt{1+|{\boldsymbol{T}}-{\boldsymbol{x}}|^2}\) provides the ray mapping of the geometrical optics limit of the Helmholtz equation \(\Delta u+k^2n^2 u=0\), where \(n\) is the refraction index and \(k\) is the wave number. Another canonical example is \(C({\boldsymbol{x}},{\boldsymbol{T}}({\boldsymbol{x}}))=|{\boldsymbol{T}}-{\boldsymbol{x}}|^2\) where the optimal \(\bar{\boldsymbol{T}}\) provides the ray mapping for the paraxial Fresnel equation.

In addition to the theoretical significance of the weighted least action principle, it also serves to solve the phase from intensity problem. The difficult problem of estimating the phase of a wave is solved from two relatively easy intensity measurements, and a single optimization problem [4].

In all the examples listed in [2] the rays serve two roles. First, they are directly related to the gradient of the phase \(s(x)\). Second, the rays are the transport vectors for the propagation of the intensity. For this reason the optimal mapping \(\bar{T}\) provides a solution to the phase from intensity problem.

Our goal in this paper is to extend the theory to more general wave problems where there are two distinct families of rays. One family consists of the normal to the wavefronts of the wave, and a different family of rays that are tangent to the direction of the energy propagation. The chief example of such wave problems is Maxwell equations, and we will concentrate on this case. However, we point out that our analysis and conclusions apply to other anisotropic wave problems, such as those arising in elasticity [5]. The geometrical optics limit of Maxwell equations in an isotropic medium is the same as the limit of the scalar wave equation, where the two families of rays coincide. Thus, we are interested here in the anisotropic case.

In the next section we will recall some essential properties of the geometrical optics limit of Maxwell equations. In section 3 we will derive the equivalence of this limit and an associated mass transport problem. This will enable us to solve the phase from intensity problem for Maxwell equations. Examples of media that are common in practice, and where the calculations can be made explicit are presented in section 4.

A comment on our notation: We use \({\boldsymbol{x}}.{\boldsymbol{p}},{\boldsymbol{q}},{\boldsymbol{T}}\), etc. to denote two dimensional vectors, and \(x,p,q\) etc. to denote three dimensional vectors.

2 The geometrical optics limit of Maxwell equations↩︎

This section follows the exposition of [6] and [7]. We consider time harmonic solutions of the Maxwell equations. The geometrical optics limit of these equations includes two families of curves: Fresnel wave normals, denoting here by \(p\) and Fresnel rays denoted here \(q\). The curves \(p\), as their name indicates, are the normals to the wavefronts. Denoting the phase by \(s(x)\) we have \(p=\nabla s\). The rays \(q\) are associated with the Poynting vectors, and are the carriers of energy. Denoting the electromagnetic energy \(G\), the transport equation is \[\nabla \cdot (qG)=0. \label{p6}\tag{4}\] The geometrical optics limit is equipped with a Hamiltonian structure that has an important geometric characterization. Denoting the Hamiltonian \(H\) (not to be confused with the magnetic field!), the Hamilton equation is \[H(\nabla s,x)=H(p,x)=0. \label{p8}\tag{5}\] The surface \(H(p,x)=0\) for fixed \(x\) is called Fresnel surface of wave normals. In an isotropic medium the surface is simply the sphere \(|p|^2- n^2=0\), where \(n\) is the refraction index of the medium. In general anisotropic media the Fresnel surface of normals consists of two nested shells, whose shape is determined by the dielectric \(\epsilon\) and induction \(\mu\) matrices. In practice most materials are magnetically isotropic, and to simplify the presentation we assume \(\mu=1\). Similarly, the \(q\) curves are associated with another surface, Fresnel surface of rays \(H^*(q,x)=0\), that is determined by the inverse matrices \(\epsilon^{-1},\, \mu^{-1}\). In section 4 below we write explicit expressions for \(H\) and \(H^*\) in a reference frame where the matrix \(\varepsilon\) is diagonal.

Maxwell equations imply an important geometric connection between \(p\) and \(q\): \[q=\frac{\nabla_p H}{p \cdot \nabla_p H}. \label{p18}\tag{6}\] This implies \(p\cdot q=1\). Except for rare nongeneric instances where \(p\) lies at the intersection of the two Fresnel sheets, we can assume that the phase satisfies one of the two shells. The dual relation to (6 ) \[p=\frac{\nabla_q H^*}{q \cdot \nabla_q H^*}\label{p18b}\tag{7}\] holds as well.

Since the characteristic equations for the differential equation \(H(\nabla s,x)=0\) are \[x_{\tau} = \nabla_p H,\;\; p_{\tau}=-\nabla_x H, \label{p14d}\tag{8}\] we obtain that when \(H=H(p)\) then \(p(\tau)=p(0)\). Since \(p\) is constant along the characteristic, the first equation of (8 ) implies that \(x(\tau)\) (and the rays \(q\) too following (6 )) are straight lines.

3 The Weighted Least Action principle↩︎

Consider two planes \((x_1,x_2,x_3=z=0),\; (x_1,x_2,x_3=z=h)\). Since \(q\) is the direction of propagation, the impact, or wave action, of the electromagnetic wave on each plane is given by \(a = G q_3\). We thus define the vector \[{\boldsymbol{v}}=(\frac{q_1}{q_3},\frac{q_2}{q_3}), \label{p19}\tag{9}\] and rewrite the transport equation (4 ) in the form \[\frac{\partial a}{\partial z}+ \nabla_{{\boldsymbol{x}}} \cdot ({\boldsymbol{v}}a)=0\label{p20}\tag{10}\] subject to \[\label{I1I2} a(0,{\boldsymbol{x}}) =I_1({\boldsymbol{x}}) \;\;, \;\;a(h,{\boldsymbol{x}})=I_2({\boldsymbol{x}}) \;.\tag{11}\]

It is convenient to use (5 ) to express \(p_3\) (=\(s_z\)) in terms of \({\boldsymbol{p}}:=(p_1,p_2)\): \[p_3=D({\boldsymbol{p}},{\boldsymbol{x}}). \label{p22}\tag{12}\] We can therefore write \({\boldsymbol{v}}\) as \[{\boldsymbol{v}}=-\nabla_{{\boldsymbol{p}}} D({\boldsymbol{p}},{\boldsymbol{x}}). \label{p24}\tag{13}\]

Let \({\boldsymbol{T}}^*:\Omega_1\rightarrow \Omega_2\) be the ray mapping obtained by integrating the orbit \[{\boldsymbol{x}}^*_z={{\boldsymbol{v}}}(z,{\boldsymbol{x}}^*),\;\;\;{\boldsymbol{x}}^*(0)={\boldsymbol{x}}, \label{veq}\tag{14}\] where \({\boldsymbol{v}}\) is determined by (13 ), and setting \[{\boldsymbol{T}}^*({\boldsymbol{x}}) = {\boldsymbol{x}}^*(h). \label{p24b}\tag{15}\] Recalling equation (10 ), we obtain \({\boldsymbol{T}}^*\in {\cal C}\). Equation (13 ) implies that the Legendre transform of \(D({\boldsymbol{p}},{\boldsymbol{x}})\) is \[D^*({\boldsymbol{v}},{\boldsymbol{x}})=max_{\boldsymbol{p}}D({\boldsymbol{p}},{\boldsymbol{x}})+{\boldsymbol{p}}\cdot {\boldsymbol{v}}\label{p26}\tag{16}\] where the equality holds for \({\boldsymbol{p}}=\nabla_{{\boldsymbol{x}}} s\), i.e. \[D^*({\boldsymbol{v}},{\boldsymbol{x}})=D(\nabla_{\boldsymbol{x}}s,{\boldsymbol{x}})+\nabla_{\boldsymbol{x}}s\cdot {\boldsymbol{v}}\label{p66}\tag{17}\] Examples of explicit forms of \(D^*\) will be given at the next section.

We are now at a position to write down the Weighted Least Action principle for the geometrical optics limit of Maxwell Equations. Consider two points \({\boldsymbol{x}}_1 \in \Omega_1,\;\; {\boldsymbol{x}}_2 \in \Omega_2\). The action of an orbit \({\boldsymbol{x}}(z)\) connecting them is defined as \[Q({\boldsymbol{x}}_1,{\boldsymbol{x}}_2)=\min \int_{0}^{h} D^*({\boldsymbol{x}}_z,{\boldsymbol{x}}) dz, \label{p30}\tag{18}\] where the minimization is over all orbits \({\boldsymbol{x}}(z)\) connecting the two endpoints. Let the wave action at \(\Omega_1\) and \(\Omega_2\) be \(I_1({\boldsymbol{x}})\) and \(I_2({\boldsymbol{x}})\) respectively. Let \({\cal C}\) be the family of mapping \({\boldsymbol{T}}\) that transform \(I_1\) to \(I_2\) in the sense of equation (2 ). Consider now the problem of minimizing the Monge, or weighted action, functional \[M(T)=\int_{\Omega_1} I_1({\boldsymbol{x}}) Q({\boldsymbol{x}},{\boldsymbol{T}}({\boldsymbol{x}}))d{\boldsymbol{x}}, \label{p32}\tag{19}\] over all maps \({\boldsymbol{T}}\in {\cal C}\).

The ray mapping \({\boldsymbol{T}}^*\) defined in equation (15 ) is a minimizer \(\bar{\boldsymbol{T}}\) of (19 ).

Consider the complete \(z\) derivative of the phase function \(s\) along any orbit \({\boldsymbol{x}}(z)\): \[\frac{d s({\boldsymbol{x}}(z),z)}{dz} = p_3 + \nabla_{\boldsymbol{x}} s \cdot {\boldsymbol{x}}_z =D(\nabla_{\boldsymbol{x}} s,{\boldsymbol{x}})+ \nabla_{\boldsymbol{x}}s \cdot {\boldsymbol{x}}_z\leq D^*({\boldsymbol{x}}_z,{\boldsymbol{x}}). \label{p34}\tag{20}\] where we used \(s_z=p_3\) in the first equality, equation (12 ) in the second equality and equation (16 ) for the inequality. If we substitute the orbit \({\boldsymbol{x}}^*\) defined in equation (14 ), and recall equation (17 ), we get an equality in (20 ): \[\frac{d s({\boldsymbol{x}}^*(z),z)}{dz} = D^*({\boldsymbol{x}}^*_z,{\boldsymbol{x}}).\]

Integrating this equation from \(z=0\) to \(z=h\) we obtain:

\[s({\boldsymbol{T}}^*({\boldsymbol{x}}),h)-s({\boldsymbol{x}},0)= \int_0^h D^*({\boldsymbol{x}}^*_z,{\boldsymbol{x}})dz \geq Q({\boldsymbol{x}}, T^*({\boldsymbol{x}}))\label{p36}\tag{21}\] where last inequality follows from the definition of the action \(Q\). Multiplying both sides by \(I_1({\boldsymbol{x}})\) and integrating we obtain, using the mass transport property (1 , 15 ) of \(T^*\) \[\int s({\boldsymbol{x}},h)I_2({\boldsymbol{x}}) - s({\boldsymbol{x}},0)I_1({\boldsymbol{x}})=\int \left(\int_0^h D^*({\boldsymbol{x}}^*_z,{\boldsymbol{x}}) dz\right) I_1({\boldsymbol{x}})d{\boldsymbol{x}}\geq \int Q({\boldsymbol{x}},{\boldsymbol{T}}^*({\boldsymbol{x}}))I_1({\boldsymbol{x}})d{\boldsymbol{x}}, \label{p38}\tag{22}\]

Now, let \(\bar{{\boldsymbol{T}}}({\boldsymbol{x}})\) be the optimal mapping in the theorem, and let \(\bar{{\boldsymbol{x}}}(z)\) be the orbit connecting \({\boldsymbol{x}}\) and \(\bar{{\boldsymbol{T}}}(x)\) and achieves the optimal \(Q({\boldsymbol{x}},\bar{{\boldsymbol{T}}}({\boldsymbol{x}}))\). If we integrate the flow (20 ) along \(\bar{\boldsymbol{x}}(z)\) we obtain \[\int s({\boldsymbol{x}},h)I_2({\boldsymbol{x}}) - s({\boldsymbol{x}},0)I_1({\boldsymbol{x}}) \leq \int \left(\int_0^h D^*(\bar{\boldsymbol{x}}_z,{\boldsymbol{x}}) dz\right) = \int Q({\boldsymbol{x}},\bar{{\boldsymbol{T}}}({\boldsymbol{x}}))I_1({\boldsymbol{x}})d{\boldsymbol{x}}. \label{p40}\tag{23}\] Equations (22 ) and (23 ) together imply \(\int Q({\boldsymbol{x}},\bar{{\boldsymbol{T}}}({\boldsymbol{x}}))I_1({\boldsymbol{x}})d{\boldsymbol{x}}\geq \int Q({\boldsymbol{x}},{\boldsymbol{T}}^*({\boldsymbol{x}}))I_1({\boldsymbol{x}})d{\boldsymbol{x}}\), hence \({\boldsymbol{T}}^*=\bar {\boldsymbol{T}}\) is the minimal map.

We proceed to apply the weighted least action principle to solve the problem of determining the phase of the wave from two measurements of the energy at two planes along a preferred propagation axis. In the scalar wave equation this problem can be solved directly from the weighted least action principle, since the ray directions obtained from the optimal mapping \(\bar{{\boldsymbol{T}}}\) are directly related to the normals to the wavefronts. In the present anisotropic medium this is no longer the case. We demonstrate the new solution method for the homogenous case where \(H=H(p)\), or \(D=D({\boldsymbol{p}})\). As stated in the preceding section, when \(D\) depends only on \({\boldsymbol{p}}\), the rays are straight lines. Therefore we can write \[{\boldsymbol{v}}=(\bar{{\boldsymbol{T}}}({\boldsymbol{x}})-{\boldsymbol{x}})/h. \label{p44}\tag{24}\] The optimization problem (19 ) becomes to find minimum of \[M({\boldsymbol{T}})=\int_{\Omega_1} I_1({\boldsymbol{x}}) D^*({\boldsymbol{T}}({\boldsymbol{x}})-{\boldsymbol{x}})d{\boldsymbol{x}}, \label{p46}\tag{25}\] over all maps \({\boldsymbol{T}}\in {\cal C}\). Since \(D({\boldsymbol{p}})\) is a concave function, \(D^*\) is convex and the minimization problem is well-defined and the minimizer is unique (see, e.g [8] sec. 1.3 ). The optimal mapping \(\bar{{\boldsymbol{T}}}\) can be found by a variety of known numerical schemes (e.g. [8], Ch 6). The optimal mapping \(\bar{{\boldsymbol{T}}}\) implies, through equation (24 ), two components for the 3d vector \(q\). A third equation is the Fresnel ray surface \(H^*(q)=0\). This determines completely the vectors \(q({\boldsymbol{x}})\) in the domain \(\Omega_1\). Equation (7 ) determines then the vectors \(p({\boldsymbol{x}})\) on \(\Omega_1\) and thus the initial phase.

It is important to note that while the solution of the Monge problem (25 ) is unique, the mapping \({\boldsymbol{T}}\) realizing the illumination problem may not be so. Indeed, as observed in [9], each critical point of the functional \(M\) is a legitimate solution. Moreover, each sheet of the Fresnel surfaces implies its own solution.

4 Explicit expression for the action \(Q\)↩︎

We present a few examples where the action \(Q\), and thus the functional \(M\), can be computed explicitly. Assume that the dielectric matrix \(\varepsilon\) is diagonalized with the third direction \(z=x_3\) being the along the normal to the planes \(z=0\) and \(z=h\). Using the eigenvectors to determine the reference frame, the matrix \(\varepsilon\) is diagonal with elements \((\varepsilon_1,\varepsilon_2,\varepsilon_3)\). Introducing the two auxiliary functions [6], [7] \[\begin{align} \Phi= \frac{1}{2}\left( (\frac{1}{\varepsilon_2}+\frac{1}{\varepsilon_3}) p_1^2+(\frac{1}{\varepsilon_1}+\frac{1}{\varepsilon_3}) p_2^2+ (\frac{1}{\varepsilon_1}+\frac{1}{\varepsilon_2}) p_3^2 \right), \tag{26} \\ \Psi=(p_1^2+p_2^2+p_3^2)\left(\frac{p_1^2}{\varepsilon_2 \varepsilon_3}+ \frac{p_2^2}{\varepsilon_1 \varepsilon_3}+\frac{p_3^2}{\varepsilon_1 \varepsilon_2}\right), \tag{27} \end{align}\] the two sheets of Fresnel surface are given by \[H=\Phi-1 \pm \sqrt{\Phi^2-\Psi}=0. \label{p10}\tag{28}\] Observe that the definitions of \(\Phi\) nd \(\Psi\) imply that both sheets are homogenous of degree two. Therefore \(p \cdot \nabla_p H =2\), and thus the relation (7 ) becomes \(q=\nabla_p H/2\).

The direct computation of \(D^*\) as the Legendre transform of \(D({\boldsymbol{p}},{\boldsymbol{x}})\) is hard. Instead we compute it directly via the Fresnel ray surface equation. For this purpose we notice that the function \(D^*\) presented in equation (16 ) can be written as \[D^* = p_3+{\boldsymbol{p}}\cdot {\boldsymbol{v}}= (p_1,p_2,p_3) \cdot (q_1/q_3, q_2/q_3 , 1) = p \cdot q /q_3. \label{p51}\tag{29}\] Recalling the identity \(p \cdot q=1\) the action \(Q\) of a ray connecting \({\boldsymbol{x}}_1 \in \Omega_1,\;\; {\boldsymbol{x}}_2 \in \Omega_2\) can be written as \[Q({\boldsymbol{x}}_1,{\boldsymbol{x}}_2)=\min \int_{0}^{h} \frac{dz}{q_3({\boldsymbol{x}}_z,{\boldsymbol{x}})}, \label{p53}\tag{30}\]

To find \(q_3({\boldsymbol{v}},{\boldsymbol{x}})\) we write the two sheets of the Fresnel ray surface \(H^*(q,x)=0\) in the form \[q_3^2(\alpha \pm \sqrt{\alpha^2-\beta})=1, \label{p55}\tag{31}\] where (suppressing the \(x\)-dependence of the dielectric coefficients) \[\begin{align} \alpha = \frac{1}{2}\left((\varepsilon_2+\varepsilon_3){\boldsymbol{v}}_1^2+(\varepsilon_1+\varepsilon_3){\boldsymbol{v}}_2^2+(\varepsilon_1+\varepsilon_2)\right), \tag{32} \\ \beta=({\boldsymbol{v}}_1^2+{\boldsymbol{v}}_2^2+1)\left(\varepsilon_2 \varepsilon_3 {\boldsymbol{v}}_1^2+\varepsilon_1\varepsilon_3 {\boldsymbol{v}}_2^2 +\varepsilon_1\varepsilon_2 \right). \tag{33} \end{align}\] For example, in the isotropic case \(\varepsilon_i\equiv \varepsilon=n^2\), we ottain \[Q({\boldsymbol{x}}_1,{\boldsymbol{x}}_2)=\min \int_{0}^{h} n \sqrt{1+{\boldsymbol{x}}_z^2}dz. \label{xwqupkzs}\tag{34}\] Another special case is uniaxial material where \(\varepsilon_1=\varepsilon_2 \neq \varepsilon_3\). Here one sheet of the Fresnel ray surface is \(q^2=\varepsilon^{-2}\) and \(Q\) is as in equation (33 ). The action for the other sheet is easily found from equation (31 ): \[Q({\boldsymbol{x}}_1,{\boldsymbol{x}}_2)=\min \int_{0}^{h} n \sqrt{1+\lambda{\boldsymbol{x}}_z^2}dz. \label{p61}\tag{35}\] where \(\varepsilon_1=n^2\) and \(\lambda=\varepsilon_3/\varepsilon_1\).

5 Conclusions↩︎

We derived a weighted least action principle for the anisotropic Maxwell equations, as well as for other anisotropic dispersive vector-valued wave problems such as the Cauchy equations of linear elasticity. The theory here differs from the earlier derivation of this principle for scalar (or complex) dispersive waves since in the anisotropic case there are two systems of rays. We used the Hamiltonian structure of Maxwell equations to define a Monge functional related to Fresnel rays. Then, the duality between the Fresnel wave normal surface and the Fresnel ray surface was applied, together with the new least action principle, to solve the problem of phase from intensity. We explicitly computed the cost function for the new variational principle for several common situations.

Finally, we comment that Monge formulated the problem for civil engineering applications and in his original paper he used the cost function \(C({\boldsymbol{x}},{\boldsymbol{T}}({\boldsymbol{x}}))=|{\boldsymbol{T}}-{\boldsymbol{x}}|\). It is interesting to notice [10] that shortly after Monge published his paper, Dupin proposed to him that the mapping associated with this cost function might be related to optics. However, while Dupin’s intuition is remarkable, his suggestion is not valid since as pointed out above, geometrical optics is actually related to different cost functions.

References↩︎