A discontinuous Galerkin method to solve Liouville’s equation of geometrical optics

. We present an alternative method to ray tracing that is based on a phase space description of light propagation. Liouville’s equation of geometrical optics describes the evolution of the basic luminance on phase space. At an optical interface, the laws of optics describe non-local boundary conditions for the basic luminance. A discontinuous Galerkin method is employed to solve Liouville’s equation for a dielectric total internal reﬂection concentrator.


Mathematical model
In the design of optical systems for illumination purposes, ray tracing is typically used to provide insight into the effectiveness of the system.A large number of rays needs to be traced to get enough information.An alternative approach is based on a phase space description of light, where the evolution of a light ray can be described in terms of a Hamiltonian [1].A light ray is described by its position (q, z) ∈ R d+1 and momentum vector (p, p z ) ∈ R d+1 (d = 1, 2), where the momentum vector is simply the unit direction vector multiplied by the refractive index n.Assume that light is travelling in the positive z-direction, such that we can parametrise the light ray in terms of the z-coordinate.Furthermore, this allows us to write p z = n 2 − | p| 2 .The coordinates q and p describe a point in the 2d-dimensional phase space, and they evolve according to Hamilton's equations [1].
In the absence of Fresnel reflections, scattering and other losses, the energy (per time) of a beam of light remains constant when propagating along the z-axis.This energy is the luminous flux Φ.An infinitesimal element of luminous flux dΦ can be related to an energy density ρ on phase space by dΦ = ρ dU.Here, dU = dq dp represents an infinitesimal volume element on phase space, which in optics is known as étendue and the energy density ρ is called basic luminance.Conservation of étendue and luminous flux leads to the basic luminance invariance of ρ which can be expressed as The basic luminance invariance even holds when light is reflected or refracted [2,3].If we need to consider both forward and backward propagating light, then relation (1) needs to be augmented with the direction of the light, see [4].* e-mail: r.a.m.v.gestel@tue.nlFrom relation (1) one can derive Liouville's equation, which reads where H(z, q, p) = − n(z, q) 2 − | p| 2 is the Hamiltonian.At an optical interface (that is: a discontinuity in the refractive index field n) the momentum of a light ray changes discontinuously.For forward propagating light and no losses, the basic luminance remains invariant so that at an optical interface the following relation holds where the ± denote one-sided limits towards the optical interface that correspond to incident and outgoing light for − and +, respectively.Moreover, p(z + ) is computed from p(z − ) with the vectorial laws of reflection and refraction.Relation (3) states that at an optical interface different parts of phase space are in contact with each other, i.e., at an optical interface we have non-local boundary conditions.The change in momentum depends on the physical surface unit normal of an optical interface.Consequently, solving Liouville's equation for curved optical interfaces, where the normal changes as a function of the position, is in particular very challenging.

Numerical solution to Liouville's equation
In [4] we describe a discontinuous Galerkin method to solve Liouville's equation numerically for twodimensional optics (d = 1).The method is energypreserving and allows for control over the accuracy of the numerical solution by choosing discretisation parameters.
The method is computationally efficient as it can converge  rapidly to exact solutions for model problems, see [4].Furthermore, the illuminance and luminous intensity can be computed from ρ.
We apply the method to the dielectric total internal reflection concentrator (DTIRC) shown in Figure 2.For details about the design of such a device, we refer the reader to [2,5].The source is located at z = 0 and the target/exit of the device is at z = Z t = 2.648668.In Figure 1 we show the numerical solutions, where at z = 1 2 Z t one can see that light has been refracted into the medium n 1 and at z = Z t some of the light has been reflected via total internal reflection yielding the top and bottom patches.The numerical solutions are of high accuracy, as indicated by the good resolution of high gradients present in the solution.Moreover, the numerical solution can capture the discontinuity across the optical interfaces perfectly.In the solution process we naturally compute ρ at certain z-values when stepping towards the target.Thus, in addition to the solution at the target we also know intermediate distributions, hence, providing a complete picture of light propagation through the optical system.

Figure 1 :
Figure 1: Distributions of ρ for the DTIRC

Figure 2 :
Figure 2: A DTIRC and a couple of light rays.The gray colour represents a refractive index n 1 = 1.5 and the white colour the background medium with n 0 = 1