An energy conservative hp-scheme for light propagation using Liouville’s equation for geometrical optics

In this contribution an alternative method to standard forward ray-tracing is briefly outlined. The method is based on a phase-space description of light propagating through an optical system. The propagation of light rays are governed by Hamilton’s equations. Conservation of energy and étendue for a beam of light, allow us to derive a Liouville’s equation for the energy propagation through an optical system. Liouville’s equation is solved numerically using an hp-adaptive scheme, which for a smooth refractive index field is energy conservative. A proper treatment of optical interfaces ensures that the scheme is energy conservative over the full domain.


Introduction
Standard forward ray-tracing methods for geometrical optics, based on Monte Carlo methods with bin counting, are in general not energy conservative and exhibit a slow convergence. A recently new approach is based on a phasespace description of light, where each light ray propagating through an optical system evolves according to an optical Hamiltonian [1,2]. The 'time' parameter for Hamilton's equations can be the arc length of the ray, or in our case it will be the distance along a chosen optical axis, denoted by z. A light ray in phase-space is described by position q ∈ R d and momentum p ∈ R d coordinates, where the momentum is related to the direction of a light ray. For two-dimensional optics d = 1, while for three-dimensional optics d = 2. Furthermore, the momentum is restricted to Descartes' disk, by |p| ≤ n(z, q) where n is the refractive index field.
The Hamiltonian description of light allows us to apply Liouville's theorem, where we now consider a differential element of phase-space. In the context of optics such elements in phase-space are called étendue. It describes a beam of light instead of individual light rays. Moreover, in a lossless optical system, energy must be conserved, implying that the luminous flux of a beam must be conserved when moving through an optical system. Furthermore étendue is conserved. Étendue dU and luminous flux dΦ are related by another invariant, known as basicluminance ρ which is defined by The basic-luminance is measured in lumen per projected area per steradian, and is a function of the distance along * e-mail: r.a.m.v.gestel@tue.nl the optical axis z, the position coordinates q and momentum coordinates p, i.e., ρ = ρ(z, q, p). The invariance of ρ implies This allows us to derive a Liouville's equation for the basic-luminance, which reads where h = h(z, q, p) = − n(z, q) 2 − | p| 2 is the optical Hamiltonian. Liouville's equation holds where h is sufficiently smooth. Near an optical interface h will be discontinuous and we have to apply Snell's law or the law of specular reflection [2]. A proper treatment of the optical interfaces is needed, since Snell's law and the law of specular reflection cause for non-local boundary conditions in phase-space, at which the momentum p is discontinuous.

Numerical method and results
As a proof of concept we restrict ourselves to twodimensional optics (d = 1). In general, Liouville's equation cannot be solved for complex systems, hence a numerical method is required. For spatial (phase-space) discretisation we apply the discontinuous Galerkin-spectral element method (DG-SEM) by Kopriva [3], while for z-integration we use a low-storage fourth order explicit Runge-Kutta [4] method. The spatial discretisation is hpadaptive, meaning we can increase the accuracy either by using more elements via h-refinement, or by increasing the order of accuracy via p-refinement. For light moving towards an optical interface, there is no condition on ρ, whilst for light moving away from an optical interface  Figure 1. Initial basic-luminance distribution ρ(z, q, p) at z = 0.
we have a Dirichlet boundary condition where the value is to be determined from the values on the other side of the interface, according to Snell's law or the law of specular reflection. However, simply taking the values from the other side does not suffice. A special approach related to the treatment of non-conforming geometries [5] is required, with an additional constraint for luminous flux conservation.
Test case: As a proof of concept we will apply the strategy outlined above to a test-case dubbed the 'bucket of water' from [6]. First, Liouville's equation is made dimensionless by some appropriate scaling, such that q ∈ [−1, +1] and that the maximum value of the initial condition is equal to 1. The test-case consists of a simple jump in the refractive index, i.e, n(q) = 1.4 if q ≤ 0 and n(q) = 1 if q > 0. Note that the optical axis is parallel to the optical interface. As an initial condition, we take the phase-space distribution shown in Figure 1, and after solving Liouville's equation till z = 0.7 we obtain the distribution displayed in Figure 2. Note that the distribution is perfectly discontinuous along the optical interface. The distribution involves three different segments, the first one is the source part in {q ≤ 0, p ≥ 0}, which has only propagated through the medium with refractive index n 1 . The second one is in the so-called reflected region {q ≤ 0, p < 0}, where a part of the distribution ends up due to reflection. The third segment is in the refracted region {q > 0, p ≥ 0}, where a part of the distribution ends after being refracted at the optical interface. Moreover, to prove energy conservation, the relative error in the total luminous flux as a function of z is plotted in Figure 3. We observe energy conservation up to machine precision.

Conclusion
The described variable order method for light propagation through an optical system based on Liouville's equation is promising. A special treatment of the optical interfaces ensures that the scheme is fully energy conservative. Future research will be conducted to extend the method to