An Energy Conservative hp-method for Liouville’s Equation of Geometrical Optics

Liouville’s equation on phase space in geometrical optics describes the evolution of an energy distribution through an optical system, which is discontinuous across optical interfaces. The discontinuous Galerkin spectral element method is conservative and can achieve higher order of convergence locally, making it a suitable method for this equation. When dealing with optical interfaces in phase space, non-local boundary conditions arise. Besides being a difficulty in itself, these non-local boundary conditions must also satisfy energy conservation constraints. To this end, we introduce an energy conservative treatment of optical interfaces. Numerical experiments are performed to prove that the method obeys energy conservation. Furthermore, the method is compared to the industry standard ray tracing. The numerical experiments show that the discontinuous Galerkin spectral element method outperforms ray tracing by reducing the computation time by up to three orders of magnitude for an error of 10-6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-6}$$\end{document}.


Introduction
Illumination optics deals with the design of optical components for various applications, like LED lighting [28] and automotive headlamps [10,36]. The design of these optical components requires a different approach than used in imaging optics, as imaging effects are highly undesirable [9]. Light propagation through an optical system is usually computed using ray tracing [14]. This means computing the evolution of many rays through an optical Ghattas. Across an optical interface the numerical fluxes are discontinuous and therefore we have to take a different approach. Inspired by [23], we present a method that directly incorporates the laws of optics and obeys energy conservation.
The article is outlined as follows: in Sect. 2 we discuss the conserved quantities in an optical system and Liouville's equation, and in Sect. 3 we discuss the DGSEM. In Sect. 4 we discuss the energy conservative treatment of the optical interfaces, and in Sect. 5 we present numerical experiments proving energy conservation for two examples. The first example features a smooth refractive index field, while the second example is a test case featuring a discontinuity in the refractive index described by van Lith et al. in [33]. In the latter example, we compare the DGSEM for solving Liouville's equation to quasi-Monte Carlo ray tracing for obtaining the illuminance. Finally we present our conclusions in Sect. 6.

Conserved Quantities and Liouville's Equation
In optics we consider the transfer of luminous flux between surfaces. A source emits a beam of radiation or light, carrying a finite amount of luminous flux denoted by Φ. In the absence of losses by absorption or scattering in an optical system, the total flux Φ is conserved, i.e., energy throughout the optical system is conserved. A related quantity is the luminance denoted by ρ * , which is defined as [9,25] where dΦ is an infinitesimal amount of flux carried by an infinitesimal beam, dω is an element of solid angle subtended at the center of the source by the area at the detector d A, and d A cos θ the projected area perpendicular to the beam, i.e., d A is an element of the surface area of the detector and θ describes the angle between the normal of the detector and the beam. The solid angle describes a cone on the unit sphere with the center of the source as its vertex and dω the area on the unit sphere subtended by the cone. Another important quantity that is also conserved in an optical system is étendue, which is defined by [9] dU := n 2 d A cos θ dω, (2) where n is the refractive index in which the beam is immersed. This allows us to write the luminance as ρ * = n 2 dΦ dU .
When a beam is propagating through a homogeneous medium the luminance ρ * is conserved, as is implied by conservation of energy and conservation of étendue. When a beam of light strikes an optical interface, e.g., a lens or a mirror, the beam is subject to Snell's law of refraction or the law of specular reflection. In the case where the beam is refracted, e.g., a transition from a medium with refractive index n 1 to a medium with refractive index n 2 , the luminance is not conserved. Assuming no Fresnel reflections, and applying conservation of energy and étendue, we obtain [25,26] where ρ * i and ρ * t describe the incident and transmitted luminance, respectively. The quantity ρ * /n 2 , known as basic luminance is conserved for refractions, cf. (4). A similar result can be derived for reflections, where the refractive indices are equal. Relation (4) will be referred to as basic luminance invariance. For a complete derivation including Fresnel reflections see [9,25,26].
The definitions of luminance, étendue and basic luminance invariance described above hold for three-dimensional optics, whereas in two-dimensional optics the definitions are slightly altered. For more details, see [9]. In summary, we denote the basic luminance for both two-and three-dimensional optics by ρ, which is defined by where the étendue dU for two-and three-dimensional systems reads [9] dU = n 2 d A cos θ dω for 3D optics, n dl cos θ dθ for 2D optics, where for 2D optics θ denotes the angle between the normal of the detector and the beam, and dθ an element of angle subtended at the center of the source by the infinitesimal line segment at the detector dl. The basic luminance is related to the luminance by ρ = ρ * /n 2 for three-dimensional optics, whereas for two-dimensional optics ρ = ρ * /n.

Liouville's Equation
In geometrical optics the evolution of light rays in a beam of light can be cast in a Hamiltonian system, where we denote with q ∈ R d the position and p ∈ R d the momentum coordinates [35]. For two-dimensional optics d = 1, while for three-dimensional optics d = 2. Both terms together describe a point (q, p) in phase space, where the phase space P is defined as the collection of all positions q and momenta p at a certain position along the optical axis denoted by the z-coordinate. A point in phase space evolves when we move along the optical axis. The momentum p = ( p, p z ) ∈ R 3 is restricted to Descartes' sphere | p| = n(z, q) where n is the refractive index field as a function of the three-dimensional position coordinates q = (q, z) [35]. This restriction invites us to use spherical coordinates to represent the momentum vector p as p = ( p, p z ) = n (sin θ cos ϕ, sin θ sin ϕ, cos θ ) , (7) where θ represents the polar angle, describing the angle between the z-axis and p measured from the z-axis, and ϕ the azimuthal angle for describing the direction in the q-plane. Therefore, at a given position z 0 along the optical axis, one can visualise the phase space coordinates on the screen that is perpendicular to the z-axis and intersects the z-axis at z 0 , where q is the position on the screen and p describes the projection of p on the screen [32]. The restriction of the momentum for physical rays p to Descartes' sphere also implies that the two-dimensional momentum vector is restricted by | p| ≤ n, describing a region known as Descartes' disc [35]. The phase space coordinates of a light ray evolve as a function of the distance along the optical axis according to Hamilton's equations, which read with h = h(z, q, p) the optical Hamiltonian given by Here, σ ∈ {−1, 0, +1} denotes the direction of the light ray travelling along the optical axis, with σ = 0 being marginal rays that travel perpendicular to the optical axis [35]. For simplicity, we assume that all rays travel in the positive z-direction given by σ = +1.
Hamilton's Eqs. (8) hold for a single light ray, however, this may be generalised to a beam of light carrying a finite amount of energy in terms of luminous flux. The flow generated by Hamilton's equations describes canonical transformations, otherwise known as symplectic transformations, on phase space. These transformations preserve the symplectic structure of phase space [1]. In other words the phase space volume element dq 1 dq 2 d p 1 d p 2 is constant. In the context of optics this has the equivalent meaning of étendue conservation. In fact dU = dq 1 dq 2 d p 1 d p 2 , which can be obtained from the first two components of the threedimensional momentum vector described by expression (7). The Jacobian determinant of p with respect to the polar and azimuthal angles θ and ϕ, can be computed as with dω = sin θ dθ dϕ an element of solid angle. Noting that the differential area on a screen can be written as d A = dq 1 dq 2 and substituting (10) into relation (6) for 3D optics, we obtain The basic luminance invariance (4) implies that ρ remains constant if we move an arbitrary distance Δz along the optical axis, i.e., ρ (z + Δz, q(z + Δz), p(z + Δz)) = ρ(z, q(z), p(z)). (12) Note that this relation also holds when a beam of light is reflected or refracted. If the solution is sufficiently smooth, one can derive Liouville's equation by subtracting the right-hand side of (12) from its left-hand side and dividing by Δz and subsequently taking the limit Δz → 0, resulting in ∂ρ ∂z Here, we have already applied Hamilton's equations (8). The advective form of Liouville's equation (13) may be written in conservative form by assuming that h is twice differentiable, upon which we obtain where we have used that the velocity field u is divergence-free, and the superscript T denotes transpose. Note that an optical interface causes the Hamiltonian to be discontinuous. Therefore, at an optical interface, Liouville's equation does not hold and we must apply (12) together with Snell's law and/or the law of specular reflection, in the limit Δz → 0.
Solving Liouville's equation on phase space at any point z along the optical axis tells us how the basic luminance changes when light propagates through an optical system, allowing us to compute at each z-coordinate the related integral quantities such as luminous flux on the screen. The total luminous flux Φ in the optical system at z = const reads Here, the phase space dependence on the z-coordinate is denoted explicitly, since the momentum domain is restricted according to Descartes' disc. Assuming the optical system is lossless the total luminous flux should be constant, i.e., Φ(z) = Φ(0). An infinitesimal element of illuminance E is defined by Applying definition (5) for the basic luminance and relation (11), dE (16) can be rewritten as Next, integrating over momentum space results in the illuminance E(z, q), i.e., where p = ( p 1 , p 2 ) ∈ P(z) in which P(z) denotes the two-dimensional momentum space at a certain position z along the optical axis. Alternatively, an infinitesimal element of luminous intensity I is defined by Applying again definition (5) for the basic luminance and definition (6) for the étendue, dI (18) can be written as dI = ρ n 2 cos θ d A.
Subsequently using the relation for p z defined in (7) and dA = dq 1 dq 2 , followed by integration over the position coordinates on the screen, denoted by q = (q 1 , q 2 ) ∈ Q(z), we obtain With these definitions, the main quantities of interest in optics can thus be easily computed from the basic luminance, satisfying Liouville's equation. In the next section, we explore a method for solving Liouville's equation.

Derivation of DGSEM
In what follows, we restrict ourselves to two-dimensional optical systems, hence reducing the complexity from a four-dimensional phase space to a two-dimensional phase space with position coordinate q and momentum coordinate p, and the distance along the optical axis denoted by z. Note that the position and momentum are now scalar quantities, therefore we omit the bold-face notation for these quantities. Next, we outline a spatial semi-discretisation of Liouville's equation, leaving only z continuous. For the semi-discretisation we apply the discontinuous Galerkin spectral element method (DGSEM) described by Kopriva in [21], to the two-dimensional Liouville equation for ρ = ρ(z, q, p) in conservative form and the flux vector f now reads The Hamiltonian h for two-dimensional optics reduces to and consequently the velocity u reads u = 1 For phase space discretisation, the two-dimensional phase space domain P is covered with straight-sided quadrilaterals Ω k ⊂ P with k the index of the element. In a more general discretisation, the boundaries of quadrilaterals are allowed to be curved, such that curved boundaries from physical constraints can be modelled appropriately. In fact, when the refractive index field changes continuously as a function of q, then the maximum allowed momentum varies as a function of q due to the restriction of p to Descartes' sphere. This restriction can be accommodated by curved boundaries when solving Liouville's equation, see [31]. For a discussion on DGSEM with curved quadrilateral elements, see for example [8,18,21,22]. In this paper we only consider straight-sided quadrilaterals.
Each quadrilateral Ω k has four vertices {x 1 , x 2 , x 3 , x 4 } labelled in counter-clockwise direction where x = (q, p) T and we have omitted the element index (superscript k), see Fig. 1. For ease of computation, the reference square χ = [−1, 1] 2 is mapped to each quadrilateral Ω k , transforming a point in the reference domain (ξ, η) ∈ χ to a point in physical space x(ξ, η) ∈ P using the following bilinear transformation The Jacobian of the transformation is given by The divergence term in (20a) can be rewritten by applying the chain rule resulting in and f is an auxiliary flux defined by the product of the adjoint Jacobian matrix and the flux f , i.e., Applying the transformation (25) to Liouville's Eq. (20a), we obtain where ρ = ρ(z, ξ, η). The weak formulation of Liouville's equation is obtained by first multiplying the PDE (27) by the Jacobian determinant J and by a smooth test function φ, and subsequently integrating over the reference domain χ. This results in The second term is rewritten by applying the product rule and Gauss's theorem, so that wheren is the outward unit normal on ∂χ and the orientation of the closed curve ∂χ is counter-clockwise. Using this, we obtain the weak formulation of Liouville's equation on the reference domain Note that for strong solutions we require the flux to be differentiable, hence, h(z, q, p) should be twice differentiable. However, the DGSEM uses the weak form of the solution and only requires the flux to be continuous, therefore, h(z, q, p) being once continuously differentiable is sufficient. For typical optical interfaces this is not sufficient since the refractive index field is discontinuous and, therefore, h(z, q, p) and also the flux are discontinuous. In particular, for these interfaces we require a special treatment of the fluxes which we will discuss in Sect. 4.

Tools for Approximating the Solution
The solution ρ in Eq. (29) is approximated by an expansion in basis functions [21]. We choose one-dimensional basis-functions ϕ i , i = 0, ..., N , for which ϕ i (ξ j ) = δ i j holds for chosen points ξ j . Moreover, we require that the basis-functions form an orthogonal basis with respect to the standard L 2 -inner product. A suitable choice are the Lagrange polynomials defined on Gauss-Legendre nodes. In the following, we will replace ρ by the approximation where ρ i j (z) are the expansion coefficients for the chosen basis. In this paper we will restrict ourselves to using the same basis functions in both directions with an equal number of expansion coefficients, although in general this restriction is not necessary. The quadrature rule defined by Gauss-Legendre nodes {ξ i } N i=0 and corresponding weights {w i } N i=0 allows us to approximate the integral of any function g, i.e., with −1 < ξ i < 1 and w i > 0. The quadrature rule on the reference domain has a tensor product structure, hence, Thus, we place nodes inside the reference domain at the Gauss-Legendre nodes (ξ i , η j ). Furthermore, for one-dimensional integrals the Gauss-Legendre quadrature gives exact integration for at least all polynomials of degree 2N + 1. Focusing on one dimension, the Lagrange polynomials on the Gauss-Legendre nodes read which satisfy the Kronecker property, i.e., It can be readily verified, that the Lagrange polynomials defined on the Gauss-Legendre nodes are orthogonal with respect to the standard L 2 -inner product, i.e., where the integration is exact since i (ξ ) j (ξ ) is a polynomial of degree at most 2N .
A very useful property of Lagrange polynomials is that the polynomial interpolation of any function g is rather easy, i.e., where I N is the polynomial interpolation operator using N + 1 nodes. The approximation of the derivative of g is defined as the derivative of the interpolant, i.e., Note that if g is a polynomial of degree N or less, both the interpolant and the derivative are exact. In what follows, we require the derivative at the nodes, i.e., where The elements of the differentiation matrix D = (D i j ) can be found by differentiating the Lagrange polynomial followed by evaluation at the node, and thus read and the diagonal elements read due to the fact that the derivative of a constant function vanishes.

Approximating the Solution
To derive an approximation of the solution, we expand both the solution and the flux in Lagrange polynomials. The expansions read The coefficients, indicated by the index-subscript i j, in each expansion are related to the position of an element's interior node (q i j , The auxiliary flux coefficients f i j (z) are related to ρ by with u i j the transformed velocity, similarly defined to (26). Here the velocity u i j (z) = u(z, q i j , p i j ) depends on z if the refractive index n depends on z. In the following, we omit f i j 's dependence on z for ease of notation. Next, we have to approximate the integrals in Eq. (29). The test function φ is chosen to be in the same basis as the solution ρ, resulting in a Galerkin method. Therefore, taking allows us to derive (N + 1) 2 equations for the (N + 1) 2 coefficients ρ i j . Combining this together with the approximations (40a) and (40b) for ρ and f we can approximate the integrals using the Gauss-Legendre quadrature rules. Therefore, substituting the approximation (40a) in the first term of (29), we obtain Applying the Kronecker property (34) of the Lagrange polynomials, the sums reduce to Note that the integral is exact for the given combination of a bilinear mapping x(ξ, η) and Lagrangian polynomials, since then the integrand is a polynomial of degree 2N + 1 in ξ and in η. The Gauss-Legendre quadrature rule is exact for this bivariate polynomial.
For the third term in (29), we substitute the approximation (40b) and denote f = ( f , g), resulting in where we have used the definition of the differentiation matrix (39). Furthermore, we introduce the following auxiliary matrix for ease of computation. The third term then reads In what follows, we will replace the flux appearing in the boundary integral from Eq. (29) with a numerical flux F = F, G . The boundary integral can be split into four parts and For the bottom part, with η = −1, the integral can be exactly evaluated using Gauss-Legendre quadrature, such that we obtain Similarly, we can compute the other components and the result for the full boundary integral reads In the discontinuous Galerkin spectral element method the elements communicate by fluxes through the faces of each element. The solution on the boundary between two elements is allowed to be discontinuous, thus the limit towards the boundary of an element can have two values, one for each element it touches. The flux on the boundary must be replaced by a numerical flux so that the neighbouring elements can communicate. The numerical flux depends on the values of ρ just left and right of the boundary, i.e., ρ L and ρ R , which are computed by evaluating the interior solution (40a) at the boundary. For the numerical flux we take the upwind flux. Due to the transformation (26) of the flux to the reference domain, the physical upwind flux F is scaled at an edge by Δl/2 with Δl the length of the edge, such that the upwind flux F = 1 2 Δl F over an edge reads wheren denotes the outward normal vector w.r.t. the element left of the boundary.
Next, we substitute expressions (43), (45) and (47) in equation (29), so that we obtain the semi-discrete ODE system for the expansion coefficients ρ i j (z): with the numerical fluxes F = F, G given by (48). This ODE system can be solved using any numerical time integrator, e.g., the classical fourth order Runge-Kutta method. Other popular choices in the literature are explicit low-storage Runge-Kutta methods, see [6,19,34]. The discontinuous Galerkin spectral element method approximates the exact solution by an N th degree polynomial, so the global spatial error e for a typical mesh size Δx behaves as Furthermore, the scheme is restricted by stability in terms of a CFL condition. For discontinuous Galerkin methods on quadrilaterals there is no direct known bound for the CFL condition. For triangular grids the relation between the Courant number and the shape of the triangles is studied in [7,30].

Optical Interfaces
In the phase space representation, the flow of ρ describes a beam of light propagating through an optical system. When the beam hits an optical interface, the momentum p changes discontinuously according to the law of specular reflection or Snell's law of refraction. Furthermore, from the discussion in Sect. 2 we know that the total luminous flux should remain constant throughout the optical system. The numerical solution should respect the actual physics, therefore, the discontinuity at optical interfaces coupled with conservation of energy should be incorporated into the DGSEM when we solve Liouville's equation.
In the DGSEM the solution is allowed to be discontinuous across the boundary connecting two or multiple elements. Therefore, the mesh in phase space is aligned such that elements adjacent to the interface have edges that coincide with the optical interface, across which the solution is discontinuous [31]. The elements in the DGSEM communicate through numerical fluxes, hence, we have to incorporate both Snell's law and the energy conservation constraint in the numerical flux when integrating (49). In particular, for a beam of light moving towards an optical interface, i.e., the velocity is directed towards the interface, we have to leave ρ free, whilst for a beam moving away from an optical interface, i.e., the velocity is directed away from the interface, we have a Dirichlet boundary condition for ρ due to Snell's law or the law of specular reflection [31].
Refraction or reflection causes the elements to be connected in a non-trivial manner at the optical interface. For example, one single element, on the side where light is moving towards the interface, can contribute to multiple elements on the other side. This occurs because both Snell's law and the law of reflection are non-linear in the momentum p. Therefore, this requires special treatment of the numerical fluxes to ensure that the scheme obeys energy conservation.
Van Lith et al. present Snell's function in [33], which is an explicit version of Snell's law and the law of specular reflection combined on phase space. Snell's function S relates the momentum p of an incident ray to the outgoing momentump of the ray. Let n 1 be the refractive index of the incident medium and n 2 the index of the transmitted medium. Then for a generic two-dimensional optical interface in the (q, z)-plane with surface unit normal ν = (ν q , ν z ) directed towards the incident medium, Snell's function reads with the auxiliary variables δ and ψ defined by In the expression for ψ the plus sign should be taken for rays that propagate in the positive z-direction, while the minus sign should be taken for rays that propagate in the negative z-direction. Furthermore, the sign of n 2 , i.e., sgn(n 2 ), in the first case of (51a) can be used to accommodate embedded mirrors in a medium of refractive index n 1 ≥ 1 by taking n 2 = −n 1 , see [32]. Note that there is a so-called critical momentum p c when n 1 ≥ n 2 so that δ = 0. For δ < 0 all light will be reflected, referred to as total internal reflection (TIR), while for δ ≥ 0 light will be refracted. The outgoing momentum is computed asp = S( p; n 1 , n 2 , ν), for which we will frequently use the shorthand notationp = S( p) and take the other parameters as given. Furthermore, the inverse of Snell's function will also be frequently used, for which we will use the shorthand notation p = S −1 (p). This means, find the momentum p such thatp = S( p). For example, for refraction the inverse reads [33] Snell's function (51) combined with (12) results in [33] ρ where p + = S( p − ; n 1 , n 2 , ν) and the ± denote one-sided limits towards the optical interface. This relation allows us to relate the basic luminance ρ on both sides of the interface.
To elaborate the energy conservation constraint, we consider the following flat optical interface parallel to the z-axis Note that the optical interface in phase space is represented by a line parallel to the p-axis.
The optical interface has two sides where on one side the normal in phase space is directed towards q < q 0 and describes the part with refractive index n 1 , whereas on the other side the normal is directed towards q > q 0 and describes the part with refractive index n 2 . The normal in phase space is given byn = (±1, 0) with the plus sign for the direction towards q > q 0 . Since the optical interface is represented by a line parallel to the p-axis at some constant q-value, only the q-component of the flux (20b) needs to be considered, i.e., cf. (22).
In what follows, we assume that light is initially in the medium with refractive index n 1 . We partition the optical interface, represented in phase space, into line segments for both sides of the optical interface. The partitioning is based on whether light is moving towards or away from the optical interface.
The line segment on the side of the optical interface with velocity directed towards the optical interface is denoted L. The line segment L describes the incoming momentum of light from the medium with refractive index n 1 , due to the assumption of light being initially in this medium.
The line segments just on the optical interface with velocity directed away from the optical interface is denoted R. The line segments R describe the outgoing momentum of light, and is further split into two parts, i.e., R = R R ∪ R T . The line segment denoted R R represents the momentum of light after total internal reflection, and is part of the optical interface where the refractive index is n 1 , while the line segment denoted R T is in the medium with refractive index n 2 , and represents the momentum after transmission.
To distinguish the momentum taken from either side of the optical interface, we write p ∈ L andp ∈ R. First, consider the integral of the flux entering an arbitrary momentum interval [p 1 ,p 2 ] ⊆ R T . The integral reads where q + 0 denotes the limit towards the optical interface from the line segment R T . Relation (53) implies where q − 0 denotes the limit towards the optical interface from the line segment L. Subsequently, we transform the integral usingp = S( p) = S( p; n 1 , n 2 , ν) resulting in The relation for reflection can be derived similarly by considering the integral of the flux entering an arbitrary momentum interval [p 3 ,p 4 ] ⊆ R R . We obtain the relation The relations (58) and (59) describe how the fluxes leaving L are related to the fluxes entering R T or R R , respectively. Henceforth, they are known as energy conservation constraints. For the particular optical interface (54) the constraints can be simplified. Since, we assumed that light was initially in the medium with refractive index n 1 , the optical interface normal on the full position space is equal to ν = (ν q , ν z ) = (−1, 0). Snell's function (51) reduces for this flat interface tō with p c = n 2 1 − n 2 2 . Then, the energy conservation constraint for R T , given by (58), can be simplified by noting that for refraction the following relations hold Hence, we obtain and similarly for reflection the constraint for R R , given by (59), reduces to The balances (61) have to be combined with relation (53) to ensure the scheme conserves energy. However, the coupling between the line segments L and R is not straightforward, which will be discussed in the next section.

Conservative Handling of Fluxes
First, consider only the refractive part of the optical interface. Elements adjacent to the optical interface in phase space are shown in Fig. 3a. These elements have edges on the optical interface and these edges are denoted by L i (i = 1, 2) and R j ( j = 1, 2, 3, 4). Due to refraction, the value of ρ in the elements that contain R j as an edge is determined by the flow through the elements that contain L 1 and L 2 . In fact, taking a closer look at how Snell's function connects the line segments from L to R in momentum space at the optical interface, we obtain for example Fig. 3b. In Fig. 3b L i and R j denote line segments along the optical interface. The basic luminance ρ along the line segments L i and R j are represented by their inner-element solution evaluated at the optical interface. To simplify notation, we denote these polynomials along the optical interface by ρ L i ( p) with i = 1, 2 and ρ R j ( p) with j = 1, 2, 3, 4. For example: where ζ = ζ( p) denotes the line segment's local reference coordinate along the interface. In Fig. 3b Hence, the endpoints of these line segments are found applying Snell's function to the endpoints of L i , i.e.,p L i = S( p L i ). Note that due to Snell's function, the line segments L i are stretched or compressed in the momentum direction. Computing the endpointsp L i allows us to determine which line segments before the optical interface contribute to a single line segment after the optical interface. From the figure we see that part of L 1 contributes to R 2 (the blue coloured region). Therefore, a relation connecting ρ L 1 ( p) and ρ R 2 ( p) on opposite sides of the optical interface must be found. Hence, as a first step applying relation (53) to a polynomial on L i , allows us to find the corresponding ρ onL i , i.e., The coupling between line segments that do not exactly match, as shown in Fig. 3b, is similar to what is known as a geometrically non-conforming mesh [20,23]. In [23] the authors describe a discontinuous Galerkin method for non-conforming meshes, applied to Maxwell's equations that form a hyperbolic system of PDEs. In their approach to treating non-conforming interfaces the solutions are first transferred to an intermediate construct called a 'mortar', and on this mortar the numerical fluxes are computed and transferred back to the corresponding elements. The transfer of the solutions and numerical fluxes is done using a least-squares matching, with integrals evaluated using Gauss-Legendre quadrature [5].
We will take a slightly different approach since in Liouville's equation for optics the flux f is discontinuous across an optical interface. Instead, relation (53) and Snell's function together describe how ρ transforms across an optical interface, cf. (64), therefore, a leastsquares matching of the polynomials describing ρ along either side of the interface is used with Snell's function directly incorporated and an additional constraint is used to satisfy energy conservation.
For the reflective part of a flat interface that is parallel to the z-axis, Snell's function reduces top = −p, see (60). The conservative treatment of these types of optical interfaces is easily accommodated by choosing a mesh such that the elements and nodes are symmetric with respect to the line p = 0 and, therefore, the constraint (61b) is easily satisfied. Due to this choice of mesh each nodep j ∈ R R will exactly correspond to −p j = p j ∈ L and a point-by-point transfer of ρ can be made. Henceforth, the following exposition of the method describes the method considering only refraction.

Contribution from One Element
From Fig. 3b we see that the line segment R 2 only depends on the solution in L 1 . The polynomial ρ R 2 must thus be computed from the polynomial ρ L 1 with the additional constraint of energy conservation. That is, the integral of the flux within the blue interval on either side of the optical interface should be equal analogous to equation (61a). Therefore, the constrained least-squares approximation reads and the momenta on both sides are related by p R i := S −1 (p R i ), see Fig. 3b. Furthermore, the numerical fluxes are defined as expansions in the Lagrange polynomial basis on Gauss-Legendre nodes, similar to (62), with flux coefficients F j := u j ρ j . The minimisation of the integral in (65a) requires finding a polynomial that matches in the least-squares sense, while the constraint (65b) ensures that the scheme conserves energy.
The integrals in the constrained minimisation problem (65) are transformed to reference line segments, so that we can compute the integrals using Gauss-Legendre quadrature. To be more specific, the integral on the LHS of (65b) and the integral in (65a) are transformed to the reference line segment along R 2 , while the integral on the RHS of (65b) is transformed to the reference line segment along L 1 . Omitting the element's subscripts, applying relation (64) and introducing an auxiliary function Ξ for ease of notation, we obtain where Δp R :=p R 3 −p R 2 and Δp L := p L 2 − p L 1 . Furthermore, the coefficients σ L ∈ [−1, 1] and λ L ∈ [0, 2] denote the offset and scaling in L 1 's reference frame, such that p(σ L ) = p R 2 and p(σ L + λ L ) = p R 3 in L 1 . Finally, the auxiliary function Ξ reads This function relates the reference frame coordinates for a momentum interval [p R ,p R + Δp R ] past the optical interface to the reference frame coordinates on a momentum interval [ p L , p L + Δp L ] before the optical interface.
Next, we write the constrained minimisation problem (66) in terms of a Lagrange function L with a Lagrange multiplier μ for the energy conservation constraint, i.e., The coefficients ρ R j for the polynomial ρ R ∈ P N can then be computed by solving Recalling that both F L and F R are written as expansions in the Lagrange polynomial basis on Gauss-Legendre nodes, we obtain for the energy conservation constraint (69b) where we have used F j = u j ρ j . To evaluate the second integral, we transform it to the reference interval [−1, 1] using ζ(ξ) = σ L + 1 2 λ L (1 + ξ ). Next, we replace both integrals with Gauss-Legendre quadrature to find the exact values, since both integrands are at most N th degree polynomials, therefore we obtain with ξ k and w k the Gauss-Legendre nodes and weights, respectively. Recalling again that the polynomials ρ L and ρ R are written as an expansion in a Lagrange polynomial basis, cf. (62), we can rewrite the equations (69a) to (72) For the evaluation of the first integral, we introduce a generic auxiliary variable S i j , given by with and S i j given by (73). The integral M i j describes the orthogonality of the Lagrange polynomials on Gauss-Legendre nodes, and therefore, is easily evaluated, cf. (35), to be The coefficients M i j and S i j are elements of the matrices M, S ∈ R (N +1)×(N +1) . The matrix M is simply a diagonal matrix containing the Gauss-Legendre quadrature weights, i.e., M = diag(w) with w = (w 0 , w 1 , ..., w N ) T . By defining as the components of the vectors α R and β L , we can write the linear system given by (74) and (71) for ρ R = (ρ R 0 , ρ R 1 , ..., ρ R N ) T and μ compactly in matrix-vector form: where we take the arguments for S as understood. Furthermore, • denotes the Hadamard product for vectors a = (a 0 , a 1 , ..., a N ) T and b = (b 0 , b 1 , ..., b N ) T , i.e., a • b = (a 0 b 0 , a 1 b 1 , ..., a N b N ) T . Let us denote the matrix on the LHS by A: The determinant of this matrix reads see Appendix A for details. Therefore, the matrix is only singular if all coefficients α R i = 0, or equivalent if all velocities u i = 0, which would mean no flux can enter the element from that side. Hence, we can safely assume that the matrix A is regular.
An analytical inverse for the matrix A is derived in Appendix A, and reads where the coefficients of the matrix B read Now, we can directly obtain an expression for the Dirichlet boundary condition values ρ R in terms of ρ L , i.e., Note that for problems where the refractive index n does not depend on z, the coefficient matrix C relating ρ R and ρ L can be pre-computed and re-used during integration along the z-axis.

Contributions from Multiple Elements
From Fig. 3b we see that the element R 3 depends on bothL 1 andL 2 . The idea remains the same, i.e., to use a least-squares matching with a constraint to ensure that the scheme is energy conservative. The constrained least-squares problem for R 3 reads where p R 34 := S −1 (p R 34 ) andp R 34 is the momentum value where the intervalsL 1 andL 2 meet, see Fig. 3b. Furthermore, ρL contains the contributions from ρ L 1 and ρ L 2 , and is defined by The integrals in (83) are transformed to their respective line segments, e.g., the integral on the LHS of (83b) and the integral in (83a) are transformed to the reference interval [−1, 1] along R 3 , such that we obtain where we write R instead of R 3 for brevity, and κ is defined such that p(κ) =p R 34 in R 3 and Δp R :=p R 4 −p R 3 , Δp L 1 := p L 2 − p L 1 and Δp L 2 := p L 3 − p L 2 . Note that σ L 1 + λ L 1 = 1 and σ L 2 = −1, however, for illustration purposes we will keep using the variables rather than these values. The Lagrange function L for this constrained minimisation problem reads The coefficients ρ R j for the polynomial ρ R ∈ P N can be found by solving Following the same steps as in Sect. 4.2, we obtain the following system of equations with (87d) The linear system described by (87) can once again be assembled into a matrix-vector form: where we have used the shorthand notation S L 1 = S i j (−1, 1+κ) and S L 2 = S i j (κ, 1−κ) . Note that the matrix on the LHS is exactly the same as the matrix obtained in the previous section, except for possibly different values for α R j . Therefore, we can again solve the linear system explicitly for the Dirichlet boundary condition values ρ R , resulting in cf. (82). This result can of course be generalised to K elements contributing to ρ R , resulting in

Overview
To summarise, during a z-step the numerical fluxes over the optical interface are evaluated as follows. First, the elements are identified that have an edge on the optical interface. Those elements are separated into elements with velocities directed towards the optical interface, denoted L, and elements with velocities directed away from the optical interface, denoted R. For the elements from L the solution is evaluated at edges on the optical interface. The numerical flux over the edges for the elements L can be directly computed as there is no constraint on ρ. For each element from R there is a Dirichlet boundary condition on the edge at the optical interface given by (53), that is incorporated into the numerical flux. The value for the Dirichlet boundary condition is determined from the elements L, as follows. To determine which elements from L contribute to a single R element, Snell's function is applied to the momentum boundaries of the elements L. Subsequently, the geometric quantities relating the element sizes are computed. Next, the momenta p at the quadrature nodes, for evaluation of the integral S i j , are determined. Subsequently, we apply S −1 to these nodes, and compute Ξ using (67). Hereafter, the integrals S i j are evaluated and the coefficients α L j , α R j are computed. Finally, the values for the Dirichlet boundary condition can be found from their contributing L-elements by applying (90).

Results
Numerical experiments were performed for two examples. The first example features light propagating through a gradient-index medium. The smooth refractive index field of the medium fits naturally into the DGSEM for solving Liouville's equation. For such optical systems ray-tracers usually have to resort to difficult to obtain closed-form expressions for the trajectory of the rays [2], or use symplectic integrators to solve Hamilton's equations for every ray [27]. Solving Liouville's equation with the DGSEM provides directly the energy distribution, i.e., the basic luminance ρ for the optical system. Furthermore, the method conserves energy by design.
The second example features a single optical interface. The problem exhibits both total internal reflection and refraction. At the optical interface we apply the strategy outlined in Sect. 4. Furthermore, a comparison is made between solving Liouville's equation using the DGSEM and applying quasi-Monte Carlo ray tracing [14]. The illuminance is solved using both methods and the performance of both methods is tested.

Elliptic Waveguide
As a first example, we consider the elliptic waveguide [35] which features a smooth refractive index field given by The parameters n 0 and κ are taken to be n 0 = 1.4 and κ = n 2 0 − 1. The refractive index field and several rays are shown in Fig. 4. We observe that the elliptic waveguide contains light much like an optical fibre. Hamilton's equations (8) for rays inside the elliptic waveguide The result using a 6th degree polynomial (N = 6) and K = 16 × 16 = 256 elements is shown in Fig. 6, together with the initial condition. The numerical solution at z = Z has roughly the same phase space area and is approximately a rotation of the initial condition. The scheme is also energy conservative up to machine precision, as can be seen in the plot of Fig. 7. Here, the relative error in the total luminous flux is plotted as a function of z. The total luminous flux is computed according to its definition (15) by applying a quadrature rule in agreement with the chosen polynomial degree. Furthermore, a convergence test is performed for this example by changing the number of elements K and varying the polynomial degree from N = 1, 2, ..., 6. The numerical solution is compared to the exact solution, which can be found from the trajectory of the rays given by expressions (93). The expressions describe the evolution of a ray, given the initial conditions of the ray. For the analytical solution to Liouville's equation at an arbitrary z we want to know where the ray originated from, since the phase space coordinates on the mesh are known. This amounts to tracing the ray backwards starting from an arbitrary z to z = 0, resulting in with q(z) and p(z) given in (93).
Using the exact solution (97) we can evaluate the discretisation error for which we take the L 1 -norm, i.e., where ρ DG denotes the numerical solution and ρ denotes the exact solution (97). The integrals in (98) are evaluated using Gauss-Legendre quadrature with N + 3 nodes. The convergence order γ DG is estimated from the empirical relation with C DG > 0 an arbitrary constant. The convergence data is shown in Table 1. The spatial discretisation is done using an N th degree polynomial, and therefore the spatial order of accuracy is N + 1. The temporal discretisation is done using a 4th order explicit Runge-Kutta method, where we choose Δz to be the maximum allowable step such that the temporal integration is stable. Furthermore, a uniform rectangular mesh is used, where upon mesh refinement the mesh size in each direction is halved and similarly Δz is halved to ensure stability.
The global error depends on whether the spatial or temporal discretisation errors dominate. From Table 1, we observe that the spatial discretisation error dominates for the polynomial degrees N = 1 to N = 6. Choosing a smaller Δz-step in the numerical experiments did not influence the discretisation error. The results show that we obtain the expected N + 1 order of convergence.

Bucket of Water
To illustrate that the strategy outlined in Sect 4 for handling optical interfaces is energy conservative, we apply it to a test case. The test case 'bucket of water' introduced by van Lith et al. [32,33] is a suitable choice. The refractive index for this problem is given by where we take n 1 = 1.4 and n 2 = 1. Using an initial basic luminance ρ 0 that has support D where q < 0 and p > 0 for all (q, p) ∈ D, the solution features both refraction and total internal reflection in two separate quadrants of phase space. The exact solution reads [33] ρ(z, q, p) = where p c = n 2 1 − n 2 2 ,p = −S(− p; n 2 , n 1 , − ν) with ν = (−1, 0), and The region described by {q < 0, p ≥ 0} features propagation through the medium with refractive index n 1 . The region {q < 0, − p c < p < 0} describes light that was reflected at the optical interface, and the region {q > 0, p ≥ 0} describes light that was refracted.
As an initial condition we use with ϕ m defined in (95) and on the part of the boundary of the domain that is not on the optical interface, we prescribe ρ = 0 whenever the velocity field is pointing into the domain, otherwise we leave ρ free. Since the q-position is restricted to q ∈ [−1, 1], this means that at q = ±1 we place virtual detectors that capture any luminous flux leaving the system. For the parameters in (102), we take q 0 = −0.35, σ q = 0.25, p 0 = 0.45, σ p,0 = 0.45, p 1 = 1 2 (1.3 + p c ) and σ p,1 = 1.3 − p 1 . Furthermore, we take m = 7 unless specified otherwise.
Again, the explicit 4th order RK method from the previous example is used with a constant Δz-step as determined by the stability of the temporal integration. The numerical solution is integrated from z = 0 to z = Z = 0.7 and z = 2Z , and is shown in Fig. 8, together with the used initial condition. The result was obtained using a 6th degree polynomial (N = 6) and K = 480 elements. The mesh uses only rectangular elements and is almost uniform. To easily treat the optical interface, we have shifted the elements below and above the critical momentum p c ≈ 0.98 in the p-direction such that the critical momentum is aligned with the edges of these elements. The mesh spacings for K = 480 are Δq = 0.1 and Δp ≈ 0.1. In Fig. 8 the quadrants featuring reflection and refraction can be clearly distinguished, while the solution is, as expected, perfectly discontinuous along the optical interface. Furthermore, at z = 2Z some of the solution has passed q = 1, meaning some energy has hit the detectors. We observe that a total 7.5 % of the initial luminous flux has hit the detectors at z = 2Z . Taking into account the luminous flux on the detectors, we compute the relative error in the total luminous flux as a function of z which is plotted in Fig. 9a. The plot shows that the method obeys energy conservation up to machine precision.
Furthermore, to show that the optical interface treatment does not incur any penalty on the convergence order, we compute the discretisation error for this example defined in (98). The convergence data for N = 1, 2, ..., 6 is shown in Table 2. Also for this example, we observe that the spatial discretisation error is dominant and choosing smaller Δz-steps did not result in different discretisation errors. Moreover, the expected spatial order of convergence N + 1 is obtained.
Next, we verify the exponential convergence of DGSEM by increasing the polynomial degree, whilst keeping the number of elements fixed to K = 1920 and choosing m = 28 in (102). For temporal integration a fixed number of 2 · 10 4 z-steps are performed, chosen such that the temporal integration error does not interfere with the convergence test. The result is shown in Fig. 9b and exponential convergence is observed.

Comparison with Ray Tracing
We compare forward ray tracing with bin-counting to solving Liouville's equation using the DGSEM. Solving Liouville's equation already has two advantages, i.e., it conserves energy and provides a more complete picture due to solving the luminance instead of its integrated quantities, the illuminance or luminous intensity. The latter advantage also comes at a price of having to solve a two-dimensional problem in phase space followed by integration to compute these quantities. Ray tracing on the other hand can directly use bins on a one-dimensional grid to compute either the illuminance or luminous intensity. For a fair comparison, we compute the illuminance E, defined by (17), for this test case using both quasi-Monte Carlo ray tracing and the DGSEM. For quasi-Monte Carlo ray tracing we fix the number of bins to B = 1000 and employ a uniform grid on q ∈ [−1, 1], i.e., with Δq = 2 B . The jth bin is defined by [Q j , Q j+1 ] with midpoint q j = 1 2 (Q j + Q j+1 ). The global error for quasi-Monte Carlo integration using a 2D Sobol sequence behaves as  [11]. The 2D points are in our case the initial phase space coordinates (q i , p i ) ∈ P of each ray. For more details on quasi-Monte Carlo integration, see [24]. In the bucket of water example M = N RT denotes the number of rays and we use a fixed number of bins.
For the DGSEM we compute the luminance followed by integration such that we obtain the illuminance. Ray tracing defines an average illuminance on each bin, hence, for a fair comparison we also average the illuminance for the DGSEM when computing the discretisation error. For the discretisation error we take the L 1 -norm and compare the numerical solution to the exact illuminance, which is computed by integrating the exact luminance (101) numerically up to machine precision.
Once again we take the initial condition (102) and (95) with m = 7. The illuminance computed using ray tracing with N RT = 0.64 · 10 6 rays and the illuminance obtained with DGSEM on a mesh with K = 480 elements and N = 4 are shown in Fig. 10a, together with the exact solution. The ray tracing (RT) solution is noisy, which is inherent in the method due to the quasi-random Monte Carlo process, while the DGSEM solution is almost indistinguishable from the exact solution.
The discretisation error for ray tracing for increasing number of rays is shown in Table 3, while the results for the DGSEM with increasing number of elements K is shown in Table 4. In the tables e RT and e DG denote the errors for ray tracing and solving Liouville's equation using the DGSEM, respectively, while t RT and t DG denote their respective computation times using only a single core. Furthermore, γ RT is estimated from the empirical relation while γ DG is estimated from the empirical relation (99). From the tables we observe that ray tracing uses 2.62 · 10 9 rays taking almost an hour and a half, while the DGSEM achieves roughly the same accuracy in only 8.0 seconds when using 1920 elements. Varying the polynomial degrees results in the performance graph shown in Fig. 10b. It can be observed that the DGSEM always achieves a better accuracy for N ≥ 1, compared to quasi-Monte Carlo ray tracing in the same amount of time. The DGSEM significantly outperforms ray tracing and, moreover, can achieve high accuracies in reasonable time.

Conclusions
We have solved Liouville's equation for geometrical optics, using the discontinuous Galerkin spectral element method. For smooth refractive index fields the scheme obeys energy conservation by design. At optical interfaces Snell's law of refraction or the law of specular reflection needs to be applied. Together with the basic luminance invariance, this corresponds to non-local boundary conditions in phase space. A method was presented to treat the non-local boundary conditions along the optical interface such that the scheme remains energy conservative in the presence of optical interfaces. Energy conservation is verified in an example. Moreover, in the same example the scheme was compared with forward ray tracing when computing the illuminance. Ray tracing uses bins on a one-dimensional grid to compute the illuminance, while the DGSEM has to solve a two-dimensional problem followed by integration. This still resulted in a better performance compared to ray tracing. In particular, for a fourth degree polynomial, the DGSEM has a computation time of 8.0 seconds, while ray tracing took 1 hour and 26 minutes to achieve almost the same accuracy.
At the moment Fresnel reflections [15,16] were not included. Therefore, an obvious next step will be to include this in the method by modifying the basic luminance invariance over an optical interface (53), see [26]. Moreover, only two-dimensional optics was considered in this paper. Hence, for future research we intend to extend the method to a three-dimensional optics settings. This requires a four-dimensional phase space together with the propagation along the z-coordinate, making it a five-dimensional problem. Despite the increased computational complexity due to the high dimensionality of the problem, the high convergence rates of DGSEM will still be an advantage over ray tracing, where in theory the DGSEM might achieve a better performance [31].
The second term on the RHS can be easily evaluated using a cofactor expansion along the first column, since it has all zeros except for α 0 w 0 and the remaining minor is the determinant of an identity matrix. Therefore, we obtain Repeating these steps we obtain the expression for r , defined in (80), i.e., so that