A Novel Scheme for Liouville’s Equation with a Discontinuous Hamiltonian and Applications to Geometrical Optics

A novel scheme is developed that computes numerical solutions of Liouville’s equation with a discontinuous Hamiltonian. It is assumed that the underlying Hamiltonian system has well-defined behaviour even when the Hamiltonian is discontinuous. In the case of geometrical optics such a discontinuity yields the familiar Snell’s law or the law of specular reflection. Solutions to Liouville’s equation should be constant along curves defined by the Hamiltonian system when the right-hand side is zero, i.e., no absorption or collisions. This consideration allows us to derive a new jump condition, enabling us to construct a first-order accurate scheme. Essentially, the correct physics is built into the solver. The scheme is tested in a two-dimensional optical setting with two test cases, the first using a single jump in the refractive index and the second a compound parabolic concentrator. For these two situations, the scheme outperforms the more conventional method of Monte Carlo ray tracing.


Introduction
Lighting is one of the most important, if not one of the most neglected, parts of our daily lives. Everywhere around us we apply lighting, our homes and our offices. Perhaps surprisingly, a good lighting design uses complicated optical elements. These optical elements are designed using the rules of geometrical optics. The equations that govern light rays turn out to constitute a Hamiltonian system [1]. We wish to exploit the Hamiltonian nature to come up with more efficient solvers for the numerical simulation of optical systems.
Hamiltonian systems are ubiquitous in mathematical physics. Their application ranges from classical mechanics and geometrical optics to statistical mechanics [2][3][4]. Also quantum mechanics uses the vernacular and structure of Hamiltonian systems, Schrödinger's equation being completely analogous to the Hamilton-Jacobi equation [5]. Hamiltonian systems have a rich mathematical background in Lie group theory [6]. The evolution operators involved have special properties and they are called symplectic or canonical transformations.
Symplectic methods are relevant in many fields, such as particle accelerators [7,8] and illumination optics [9,10]. There are specialized numerical integrators which are symplectic and as such conserve any quadratic form [11][12][13]. Such integrators play an important role in for instance astrophysics [14,15]. One remarkable application is the usage of a symplectic integrator to obtain an energy-preserving computational fluid dynamics scheme [16].
The main idea of the Hamiltonian formulation of mechanics is to have a single functional which generates the motion, the well-known principle of least action [17]. The principle is completely equivalent with Newton's laws [18]. In geometrical optics, the governing principle is Fermat's principle, which roughly states that light travels between two points in the shortest time. These principles can then be used to derive the equations of motion, or equivalently, Hamilton's equations.
In geometrical optics, the velocity is fixed for a light ray, which leads to the length of the momentum vector being fixed. This reduces the dimension of the momentum space. Moreover, the origin of the ray is also arbitrary, which reduces the dimension of the position space. We thus end up with a two-dimensional momentum space and a two-dimensional position space, combining to a four-dimensional phase space.

Lagrangian Versus Eulerian Picture
One way of obtaining information about a Hamiltonian system is simply to integrate the equations for a given set of initial conditions. This method requires special time integrators called symplectic integrators [19]. This approach is known as the Lagrangian flow field specification, where one follows a single trajectory in phase space. By doing this for many different initial conditions, one can characterise the flow generated by the Hamiltonian.
Another approach is to look at phase space as a whole and characterise the flow as a function of fixed coordinates and time. This approach is known as the Eulerian flow field specification and likewise, it can characterise the flow generated by the Hamiltonian [20]. Typical Eulerian problems include travel time problems which satisfy the Hamilton-Jacobi equation [21,22].
In a geometrical optics setting, the Lagrangian method is represented by a technique called ray tracing. One follows a light ray through an optical system by integrating the Hamilton equations. When a discontinuity in the refractive index is encountered, Snell's law is applied. To obtain large scale information, as is necessary in illumination optics, many rays have to be traced, on the order of several millions.
We propose another approach, which represents an Eulerian method, that involves solving Liouville's equation. This supplies us with global information about the system at once, since we only have to integrate Liouville's equation once. Liouville's equation is a hyperbolic partial differential equation, and as such there are explicit numerical methods available, which is another advantage of our approach.

Outline of the Paper
This paper proposes a new numerical method for solving Liouville's equation, in particular when encountering jumps in the Hamiltonian. To the author's knowledge, such methods are in the earliest stages of development. Jin and Wen developed a numerical solver for a time-dependent description of geometrical optics [23]. However, we shall develop a general framework for the treatment of Liouville's equation with discontinuous Hamiltonians.
In Sect. 2, we shall start with a brief outline of Hamiltonian dynamics. We will give a short introduction to canonical transformations and their implications. We introduce Liouville's equation, which is a consequence of the incompressibility of Hamiltonian flows. In the case of discontinuous changes of the Hamiltonian, the Hamilton equations imply a discrete canonical transformation on phase space.
We then switch to optics, where such discrete transformations can be easily realised in the form of a discontinuous change in refractive index, for example a lens or mirror. A short derivation of the explicit form of Snell's law in vector form is given in Sect. 3. It is applied in a special case where Liouville's equation can be solved analytically in Sect. 4, using the method of characteristics (MOC). The MOC leads to a jump condition for Liouville's equation, which we will derive in Sect. 5. The jump condition is independent of the evolution parameter or any particular form of the Hamiltonian. In Sect. 6 we derive the numerical scheme, which we illustrate for geometrical optics using a piecewise constant refractive index for simplicity. The explicit form of Snell's law together with the jump condition are at the heart of the scheme. We apply the scheme to a test case, the same as used in Sect. 4, and we show a good comparison between the analytical and numerical solutions. Next, we apply our scheme to the compound parabolic concentrator (CPC) to show that it can also handle more complex geometries. For the CPC, we will show that our scheme is more efficient in terms of computational time compared to Monte Carlo ray tracing.

Canonical Transformations and Conserved Quantities
Let us define phase space P = Q×P as the collection of positions q ∈ Q and momenta p ∈ P. At the moment, we do not specify q and p just yet, since they have different dimensions in optics and mechanics. In mechanics, phase space is conventionally a six-dimensional space, while in geometrical optics, phase space is usually four-dimensional. The Hamilton equations are given byq where h : R + × P → R is the Hamiltonian. The dot means differentiation with respect to an evolution coordinate. In mechanics, it is customary to take time as the evolution coordinate. In geometrical optics, the length down the optical axis or the arc length of a ray takes the role of time. The origin is arbitrary, but it is customary to put it at the light source. We denote this "time" by z, the dot in that case becomes d dz . Each fixed z defines a plane in physical space which we call a screen, it is perpendicular to the optical axis and intersecting at z. It becomes clear that for fixed z we can characterise a ray by a two-dimensional position vector and a two-dimensional direction vector projected on the screen. The momentum vector is proportional to the direction vector multiplied by the local refractive index. The complete ray is then represented by the position and momentum vectors as a function of z.
The optical Hamiltonian is given by where n : R + × R 2 → R is the refractive index field and σ ∈ {1, 0, −1} is the direction index, indicating in which direction a ray travels along the optical axis [3]. It is important to note that (2) represents the Hamiltonian of a single ray. To ensure physical solutions, we must restrict the momentum to |p| ≤ n(z, q). The upshot is that the momentum space consists of two disks, corresponding to σ = ±1, where one can view σ = 0 as a limiting case on either disk. We restrict ourselves to forward travelling rays, meaning we look only to rays in the disk labelled by σ = 1. Let us denote R 3 -vectors with an arrow, to distinguish them from the phase space quantities in bold face. To compute the momentum in R 3 , we need to use T is the unit direction vector of the ray in R 3 . Hence, a trajectory in phase space contains the same information as the path of a ray in real space. Likewise, the trajectory of a mechanical particle in phase space is equivalent to its trajectory in real space.
All cases of Hamiltonian evolution can be considered actions on phase space. These actions are of course governed by the Hamilton equations, and are known as canonical transformations, otherwise known as symplectic transformations or symplectic mappings. Discontinuous changes in the Hamiltonian can also constitute a canonical transformation on phase space [24].
One important theorem in Hamiltonian mechanics is Liouville's Theorem [25], the application of which asserts that volume in phase space is preserved. Hamiltonian flow is incompressible, and as such, we can derive Liouville's equation, whenever the Hamiltonian is smooth i.e., ∂ρ ∂z where ρ = ρ(z, q, p) is the phase space density, which is a particle density in mechanical settings and an energy or power density in optical settings. However, contrary to mechanics, ρ has a direct practical application in geometrical optics: it is the radiance in terms of radiometric quantities or the luminance in terms of photometric quantities. Other practical information can be obtained by integrating over a particular dimension, for instance the intensity (radiant or luminous) can be found by integrating over all positions. In fact, optical engineers are more and more working directly in phase space, not even bothering to convert the momentum to angles, see for instance Rausch et al. [9]. The derivation of Liouville's equation is based on the observation that energy is transported along rays combined with Reynolds' Transport Theorem and the incompressibility of the flow [26,27]. An alternative derivation and application to level-set methods have been performed by Qian et al. [28,29]. Hauray has shown existence and uniqueness of solutions for Liouville's equation with a mechanical Hamiltonian which has a force field of bounded variation, [30]. However, we are mainly interested in optics, particularly when the refractive index exhibits discontinuous changes, and as such (4) may not be well-defined.
Equation (4) is a hyperbolic PDE and whenever h is not smooth, the coefficients can be discontinuous or measure-valued. In general, such equations are rather difficult to deal with [31][32][33][34]. In the case when the Hamiltonian is not smooth, we can still find a canonical transformation that represents the discontinuous change [35,36]. Liouville's equation is not valid locally, but we must employ the fact that energy is transported along rays, consequently where pluses denote limits from one side of the interface and minuses the other. We will rely heavily on (5) to derive our numerical scheme. Both Liouville's equation and the jump condition (5) express continuity of ρ along rays. The difference is that Liouville's equation also needs differentiability, whereas (5) does not. Note that it characterizes the physically correct solution. In optics, energy is transported along rays, a statement which follows directly from Maxwell's equations. The quantity ρ being an energy density, conservation of energy gives us that ρ should be constant along rays, no matter if they are refracted or not. Numerical schemes aimed at recovering certain physical properties of the solution are sometimes referred to as well-balanced schemes [37][38][39].
It is important to note that the previous discussion and the jump condition (5) are also valid in a mechanical context. In fact, there is a one-to-one correspondence between rays and particle trajectories. However, energy transported by rays is a more tangible concept, as opposed to phase space density transported by particles.
For the remainder of this work, we shall be dealing with geometrical optics, where discontinuities in the Hamiltonian are more easily imagined. We shall refer to such a discontinuity in refractive index as an interface, examples of which are mirrors and lenses. Below we give a short treatment of Snell's law, which is the law governing rays at interfaces. It will allow us to apply (5), since p(z + ) and p(z − ) are related through Snell's law.

Snell's Function
We derive the explicit form of Snell's law, which is the key to obtaining exact as well as numerical solutions. Consider a ray of light incident on an interface where the refractive index jumps from n 1 to n 2 , for instance a ray travelling in air and hitting a piece of glass. The tangential momentum at the point of impact must remain constant. From this fact, we can derive laws governing the refraction of light, commonly stated as Snell's law, i.e., where 0 ≤ θ i ≤ π 2 is the angle of incidence as measured from the surface normal → ν , and 0 ≤ θ t ≤ π 2 is the angle of refraction, see Fig. 1. The convention of the angles in Snell's law is to measure them with respect the normal that points from the surface into the medium in which the ray is propagating. Furthermore, the incident angle is always positive. This means the incident and reflected ray are measured with respect to the normal → ν , but the refracted ray is measured with respect to surface normal − → ν . As the incident angle is always positive, the reflected ray angle is always negative. With this convention in place, we see that using n 2 = −n 1 in (6) gives the law of specular reflection, i.e., θ r = −θ i . In addition to (6), the incident ray, surface normal and refracted ray all lie in one plane called the plane of incidence. This statement together with (6) is enough to uniquely specify the transmitted ray. This law has been known for over a millennium [40]. We shall write the incident unit direction vector as → i , similarly for the reflected ray direction we use → r and for the transmitted ray direction we shall write → t . We can exploit the fact that → t is in the same plane as → ν and → i by constructing an orthonormal basis, see Appendix 1 for the details. Note that Snell's law only uses local information, which is to say only the surface normal at the point of impact of the incident ray is needed. Hence, for a given ray, Snell's law does not depend on the curvature of the surface, provided we use the surface normal at the intersection point of the ray and the surface. The conclusion of all these considerations is presented in the following theorem.
Theorem 1 Consider a ray travelling in a medium with refractive index n 1 with momentum p incident on an interface with (local) surface normal → ν ∈ R 3 (with → ν = 1), where the refractive index changes discontinuously to n 2 . Then, the momentum after encountering the interface is given by p = S(p; n 1 , n 2 , ν), with S defined as, where δ := n 2 2 − n 2 1 + ψ 2 and ψ := with → p and → ν being the R 3 -vectors and p and ν being their first two components, respectively.
Remark We use ν as an input parameter for S instead of → ν . The first two components of → ν provide us with enough information, since → ν = 1. In particular, we have where we have to choose the sign such that ψ ≤ 0, which follows from the angle convention discussed earlier.
The proof of Theorem 1 can be found in Appendix 1. Analogous reflection laws have been derived by Cockburn et al. [41] for a level-set approach. We will refer to S defined in (7a) as Snell's function. It can accommodate mirrors embedded in a medium of refractive index n 1 by choosing n 2 = −n 1 , since then √ δ = |ψ|, while ψ ≤ 0. Thus, even though δ ≥ 0 in this case, the refraction part of Snell's function is then equal to the reflection part. Note that when n 2 < n 1 , there is a range of angles such that δ < 0. Such rays suffer total internal reflection (TIR) and are reflected. The critical angle θ c can be found by setting n 1 sin θ c = n 2 , or equivalently δ = 0, leading to θ c := arcsin n 2 n 1 . Rays incident at angles larger than θ c suffer TIR.
Using Snell's function, we can now also tackle the reverse problem. Given a ray with momentum p in a medium with index n 2 , find a ray with momentum p in a medium with index n 1 such that, when refracted, S(p; n 1 , n 2 , ν) = p . Corollary 1 Given a ray with momentum p , we can find a ray with momentum p such that, when refracted, will end up with momentum p . The momentum p is then given by Proof We apply the Helmholtz reciprocity principle [26], also known as a backward ray, to find that we can reverse ray directions with impunity. Recalling the angle convention and applying Snell's function gives −p = S −p ; n 2 , n 1 , −ν .

Method of characteristics
One very useful tool in the analysis of hyperbolic PDEs is the method of characteristics (MOC) [42]. The MOC effectively turns a hyperbolic PDE into a set of ordinary differential equations. The idea is as follows, we introduce time dependent coordinates given by (1). We investigate the time derivative of ρ (z) := ρ (z, q(z), p(z)), given by, where the last equality comes from Liouville's equation (4). This implies ρ (z) = const. The curves in R + × P defined by the Hamilton equations are called the characteristics of Liouville's equation, as they reduce Liouville's equation, a PDE, to a set of ODEs. Hence, in geometrical optics, the characteristics of Liouville's equation are rays. Likewise, in mechanics, the characteristics are particle trajectories. We therefore see that whenever h is smooth, we can trace the characteristics all the way back to z = 0, giving ρ (z, q(z), p(z)) = ρ (0, q(0), p(0)) = ρ 0 (q(0), p(0)) .
Let us define an interface as a surface where the refractive index changes discontinuously. Furthermore, the surface must have a well-defined normal at every point. For instance, the surface may be given by q = Q(z), with Q : R + → R 2 differentiable. Thus, if we allow for an interface, we see that Snell's law (7a) combined with (5) allows us to connect two regions where the Hamiltonian is smooth. Note that q is continuous, where its z derivative is discontinuous. Furthermore, since z is the evolution coordinate, it is always continuous. These considerations lead us to the following definition.

Definition 1 (Base) Characterstics of Liouville's equation
Wherever the refractive index is sufficiently smooth, the characteristics of Liouville's equation are curves defined on R + ×P given by solutions to the Hamilton equations (1). Wherever the refractive index is discontinuous, and we have an interface, the characteristics change discontinuously according to Snell's law (7a). A base characteristic is a characteristic projected onto phase space. Thus, while a characteristic is a curve in R + × P, a base characteristic is a curve in P.
Note that in the geometrical optics setting, we define the characteristics of Liouville's equation as physical light rays. On both sides of the interface we have smooth characterstics and (11) holds, while the two sides are connected by (5). Hence, the characteristics are piecewise smooth and only discontinuous at the interfaces. Therefore, (11) holds even when interfaces are present. Jin and Wen showed that defining the analytical solution in this manner leads to a well posed problem [23].
For an application of the MOC as a solution method, see Appendix 2, where we solve a two-dimensional optical problem of a single transition given by We will also refer to this problem as the Bucket of Water Problem, since setting n 1 = 1.4 and n 2 = 1 is the situation of a water-air transition.

Gradient Jump Conditions
The method of characteristics is needed to find jump conditions on the gradient of ρ, which we will need in order to derive the numerical scheme. Through (11), we find that a linear combination of values of ρ along several characteristics is also constant. Hence, we can advance as follows to derive the jump condition we need. First, take several points in phase space close to each other and on one side of the interface. Next, we move along characteristics by a small amount Δz such that all the points cross the interface. By applying (11), we find that the finite differences we can construct using the points must be constant in z. Taking limits allows us to relate gradients on the two sides of the interface; see Fig. 2. This reasoning is simplified greatly by using a flat surface. However, since Snell's law only depends on the local gradient, the treatment of a curved surface can be derived from that of a flat surface through a suitable coordinate transform.
Theorem 2 Let h : P → R be piecewise smooth, let the interface be flat and the initial condition in (11) differentiable, i.e., ρ 0 ∈ C 1 (P). Then (5) implies where ·| ± is shorthand for evaluation at z ± , q(z ± ), p(z ± ) . Furthermore, Snell's law implies Proof 1. Let z → (q 1 (z), p 1 (z)) be a base characteristic that intersects the interface at z . Let us denote one-sided limits towards z with a superscript plus or minus. Thus, and similarly for q ± and z ± , see Fig. 2. We use as an initial condition for this characteristic For a second base characteristic, let z → (q 2 (z), p 2 (z)) be such that it has the same intersection point in Q with the interface. Both base characteristics should intersect the surface at the same point q − , but with slightly different momenta and at slightly different z. Therefore, let ε > 0 be some small number, then the second characteristic should satisfy where we use ·| − as shorthand for evaluation at (z − , q − , p − ). If we advance along z by an amount ε, we see that q 2 (z − + ε) = q − to first order. In Fig. 2, (q 1 , p 1 ) is marked by a bullet, while (q 2 , p 2 ) is marked by a square. Next, let us define the values of ρ for these two characteristics as which are both constants due to (11). Note that due to these definitions, we have, after a Taylor series expansion in ε, which also holds for all z < z , since ρ 1 and ρ 2 are constant along characteristics, see (11). Note that we can also evaluate ρ 1 and ρ 2 at z = 0, so that the right hand side of ( * ) is well-defined. 2. According to (5), ρ 1 and ρ 2 are also constant when encountering an interface. The characteristics change discontinuously, but the value of ρ 1 and ρ 2 will not change. Thus, we see that ( * ) is also true across an interface. Therefore, we can advance z by some small amount Δz such that both characteristics have crossed the interface. Choosing Δz = 2ε, we end up with Note that in a geometrical optics setting, we have that q is continuous and therefore q − = q + . Denoting p + = S(p − ; n q − , n q + , ν), we obtain for the momenta The evolution of the second characteristic is sketched in Fig. 3. Roughly speaking, during the first ε of propagation, it reaches the interface and changes discontinuously according to (7a). During the second ε of propagation, it travels through the second medium. From now on, we will assume the parameters of Snell's function understood and suppress the notation. We use a Taylor expansion on Snell's function to obtain where S(p − ) = p + and ∂S ∂p p − is essentially the Jacobi matrix of Snell's law. We take again the finite difference similar to ( * ), but now evaluated at z + Δz, giving However, since ρ 1 and ρ 2 are both constant, it follows that this must be equal to ( * ), yielding ∂h ∂p where taking the limit of ε → 0 removes the O(ε) terms. 3. Next, we find a third characteristic that passes through the same intersection point q − and at the same z-coordinate z − , but with a slightly different momentum. Hence, the third characteristic should have as an initial condition, where k ∈ P is an arbitrary vector satisfying |p − − εk| ≤ n(q − ), which ensures that the momentum p 3 is physical. We also define the value of ρ along the third characteristic, where we once again point out that ρ 3 is in fact a constant. Furthermore, analogous to ( * ), we have However, since ρ 1 and ρ 3 are constant, we find that this expression must be equal to (#), yielding Furthermore, this equality must hold for all admissible k, and letting ε → 0 yields (14). 4. We can rewrite ( * * ) into the form where applying (14) gives (13).
Remark Using the definition of the Poisson bracket, we can rewrite the jump condition (13) into a very concise form, Corollary 2 When the interface is curved, we must adjust the result of Theorem 2 as follows, where Q : R + → Q differentiable and Q(z ) gives the intersection point of a characteristic with the surface. Furthermore, (14) is also valid.
Proof A curved interface is represented by a plane in phase space moving with velocity Q (z). Hence, we perform a coordinate transform affecting the q-direction only, defined as which results in a coordinate system where the surface is standing still. Thus, in the new coordinate system, the surface is flat and we can apply Theorem 2. To finish, we note that Theorem 2 can be read as dρ dz − ∂ρ ∂z Hence, we obtain that we must apply ( * ) to ρ(z, q(z) − Q(z), p(z)), yielding (17).
Remark We may alternatively view the coordinate transform as a canonical transformation, with the new Hamiltonian given bỹ It is important to note that Theorem 2 is true in a very general setting. Its formulation is independent of the particular form of the Hamiltonian or the coordinate system. It is therefore true in mechanics as well as optics, it is even valid independent of the evolution coordinate. We can, for instance, find a Hamiltonian system describing geometrical optics in terms of real time or arc length.
The jump condition (13) is a consequence of the transformation on phase space from z − to z + being symplectic [36]. From (11), we see that ρ is constant along the characteristics. Therefore, if the influence of an interface is to move characteristics closer to each other in the q-direction, they must get further apart in the p-direction, which is reflected in the gradient of ρ through Theorem 2.

Derivation of the Scheme
In this section, we aim to develop a scheme to solve Liouville's equation numerically. The setting is two-dimensional to avoid large and cumbersome expressions. The numerical scheme for higher-dimensional systems is indeed very similar. Note that the position q and momentum p are now scalar quantities. We define the advection speeds on the screen, given by We look for approximate solutions to under the additional condition that whenever n changes discontinuously, we have with p − and p + related through Snell's law and q − = q + , z − = z + . For completeness, we quote Snell's law in two-dimensional form. Note that ν ∈ R such that |ν| ≤ 1, and we have where δ := n 2 2 − n 2 1 + ψ 2 and ψ := pν ± n 2 where the sign is to be taken such that ψ ≤ 0, which follows from the angle convention of Snell's law. We apply a grid on phase space such that we have {q i } N i=1 for the positions and p j M j=1 for the momenta. We rescale the position space such that we have q ∈ [0, 1], giving where i ∈ {1, . . . , N }. Note that the advection speeds, a and b, go to infinity for p close to n(q, z). Therefore, in many cases it is practical to choose a maximum allowed momentum in the system, p max . This corresponds roughly to setting a maximum angular aperture. The discretization of p is thus defined as for j = {1, . . . , M}. Finally, we discretize z as for t = {1, . . . , T } and z max is the total length of optical axis along which we integrate Liouville's equation. The numerical approximation of the solution is then denoted by ρ t i j ≈ ρ(z t , q i , p j ).
Let us define the positive and negative part of a number c ∈ R as follows, Whenever the refractive index is differentiable, (20) has a classical solution and the upwind scheme is straightforwardly found as where a t i j := a z t , q i , p j , b t i j := b z t , q i , p j .
As one can see, the expression in a two-dimensional optical system is already quite cumbersome, though not complicated. A three-dimensional optical system will just add four more upwind difference terms. We now wish to find a scheme that gives us the correct physical solution whenever we allow n to have discontinuities. The correction to the upwind scheme applies only locally around the interface, away from the interface the scheme will be given by (27). Thus, to illustrate our method, we use the simplest case of a piecewise constant refractive index. We choose a system that has 0 ≤ q ≤ 1, fix 1 < k < N , and let us place the interface at q k+ 1 2 = q k + 1 2 Δq, thus we have It is clear that ∂n ∂q = 0 almost everywhere, thus b(z, q, p) = 0 at all the grid points. Away from the interface, i ≥ k + 2 or i ≤ k − 1, we have a smooth refractive index, resulting in the following scheme, where now a does not depend on z so that a i j = a(q i , p j ). Even close to the interface, this scheme works as long as the upwind grid point is not on the other side of the interface. Hence, since the sign of a i j is equal to the sign of p j . Let us now consider the collection of grid points (q k , p j ) with p j < 0 and (q k+1 , p j ) with p j > 0. These are grid points that have their upwind grid point on the other side of the interface. Our scheme has to be different here since the characteristics have a jump in momentum when crossing the interface. We wish to approximate ∂ρ ∂q at the grid point (q k , p j ), which we can do by utilizing Theorem 2. The idea is straightforward, we make a Taylor expansion from both sides of the interface towards the interface. However, we keep Snell's law in mind and make our expansions about the point where the characteristic is discontinuous, see Fig. 4. This allows us to approximate the gradient at (q k , p j ) using information from the other side of the interface. The resulting scheme is summarized in the following theorem.

Theorem 3
Consider the collection of grid points such that q k , p j < 0 . Let us denote p = −S(− p j ; n 2 , n 1 , −ν), and δ ≥ 0 in (22a), where δ is given by (22b). The scheme is then given by whereã and a = p where θ = ( p − p r −1 )/Δp, with r such that p r −1 < p ≤ p r . In case when, δ < 0, reflection occurs and we have to use (32a), now with and a = p Remark The case for q k+1 , p j > 0 is similar.
Proof 1. We start by using the method of lines (MOL) approach, which is to say, we discretise space, but leave z continuous. Let us assume that j is such that δ ≥ 0, thus we have refraction. Let us further assume that the initial conditions are smooth. Furthermore, define p = −S(− p j ; n 2 , n 1 , −ν), which we shall shorten to p = −S(− p j ), such that p is the momentum that becomes p j if the characteristic traverses the interface. We define r as the unique index such that p r −1 < p ≤ p r . 2. Performing a Taylor expansion about the relevant grid points close to the interface on the left side reveals, and similarly on the right side, where h.o.t. represents (mixed) higher order terms. Next, we perform another Taylor expansion, this time for the derivative, i.e., Now we apply (13) where again, h.o.t represents mixed higher order terms. 3. We will now take a linear combination of the values of ρ at the three grid points given by ( * ) and ( * * ). We wish to find a linear combination which approximates the q-derivative of ρ at the grid point (q k , p j ). Assume that λ l ∈ R, for l = 1, 2, 3, then the upwind finite difference should satisfy where the approximation should be correct up to first-order. Using (21), we find the following system of equations, where one should note that the determinant of this linear system is non-zero. Definingã as the harmonic mean of a k j and a , see (32b), we find the solution of this system as Δz , whence we find the scheme Recognizing that ρ := θρ t k+1,r + (1 − θ)ρ t k+1,r −1 is nothing but a linear interpolation approximating ρ(z t , q k+1 , p ), we find (32a).

4.
When j is such that δ < 0, we have reflection and the derivation is largely the same. The only difference being that the upwind grid points are now given by (q k , p r ) and (q k , p r −1 ) instead of (q k+1 , p r ) and (q k+1 , p r −1 ). After doing completely similar operations, we find (32a) but now with ρ given by (33b) and a given by (33a).
Remark When the surface is curved we find that we must furthermore replace a k j with a k j − Q (z ) and a with a − Q (z ).
The correction for the advection speed enforces that the mapping from z − to z + is symplectic. In going from a high refractive index to a lower one, the transmitted light is stretched in the p-direction. The slightly higher advection speed near the interface can be interpreted as a shrinking of the q-direction, thereby preserving area in phase space.
From Theorem 3, we see that the scheme close to an interface is still an upwind scheme and in fact, (32a) has completely the same form as (30). However, we must replace both the advection speed and the value of ρ at the upwind grid point. Fortunately, the replacement values can be explicitly computed by invoking Snell's law. The scheme is in essence a first order accurate upwind scheme, with the correction only occurring near the interface. We thus expect the scheme to be globally first order accurate.
When considering variable refractive index fields together with interfaces, we can again employ Theorem 3. However, away from the interface we must now use (27) and only correct the terms approximating ∂ρ ∂q . The interface position does not depend on p and therefore the flow across an interface may only be driven by ∂h ∂ p . Thus, the flow across an interface cannot depend on the gradient of the refractive index field, provided we use a standard Cartesian discretisation. It is for this reason that we shall not consider variable index fields in the following.

Implementation of Curved Interfaces
Most optical systems have a curved interface, like any lens or focussing mirror, and consequently do not have an axis along which the refractive index field is constant. Thus, as we move along the optical axis, the refractive index field will change. The interfaces will therefore also move, in position space Q, as a function of z. Let us assume that we can describe the position of the interface with a differentiable function Q : R + → R d . Then the position of the interface defines a moving plane in phase space, or line in the two-dimensional case.
In such a case, it might occur that at one z-level, the surface is on one side of a column of grid points and at the next it is on the other side. Thus, fix some l ∈ {1, . . . , N } and consider the maximum z-level τ such that q l ≤ Q(z τ ) < q l+1 . Then, by definition, we have that Q(z τ +1 ) ≥ q l+1 and it may happen that ρ τ +1 l+1,m = 0 for all m = {1, . . . , M}, see Fig. 5a. This situation is especially likely to happen whenever the surface is a reflector.
The situation described is of course non-physical and completely due to the discretisation of phase space. However, we can correct the resulting error by interpolation across base characteristics. Hence, we keep track of when the interface moves across a column of grid points and check whether the solution needs correction. If it does need correction, we perform an interpolation across the base characteristics. Let us fix j and assume that some base characteristic undergoes a momentum change from p j to p and p r −1 < p ≤ p r . Then we have, possibly, that ρ τ +1 l,r −1 = 0, ρ τ +1 l,r = 0 and ρ τ +1 l, j = 0. We interpolate between these three points to correct ρ τ +1 l+1, j . In this way, whenever the solution is zero somewhere due to an artefact of the discretisation, we may correct it. After the interpolation in the p-direction to obtain values for ρ on q l and p , we have essentially a one-dimensional case across the base characteristic.

Results
As argued previously, a continuously varying refractive index will not influence the scheme near the interface in a substantial way. Thus, for the sake of simplicity, we present two cases which will exhibit all the features of the scheme that are different from a standard first-order upwind scheme. The first is the Bucket of Water problem, introduced in Sect. 4, which can be solved analytically fairly easily. The second test case is the two-dimensional Compound Parabolic Concentrator (CPC).
In the second test case, the CPC, we shall compare our method with Monte Carlo ray tracing, see for instance [43]. This procedure involves choosing a number of rays with randomized initial conditions, tracing their trajectory through the optic and subsequently performing a box count. The box count is often necessary because a proper interpolation is too compli-cated. Thus, to make a comparison between ray tracing and our method, we should pick a representative number of rays per grid point, the control volume of each grid point forming a box in phase space. The average number of rays per box can be quite high, depending on the desired accuracy. In practice, 100 rays per box is a typical number.

Bucket of Water
We take a simple jump in refractive index as in (29) and we pick n 1 = 1.4 and n 2 = 1. The refractive index field is given by where we pick k such that q k+ 1 2 = 1 2 + O(Δq). This case corresponds roughly to a waterair transition, and we shall at times refer to this problem as the Bucket of Water problem. One important thing to note is that the optical axis is parallel to the interface, resulting in a refractive index field which does not depend on z, see Fig. 6.
The reason we choose this particular problem is because it exhibits both refraction and TIR. At the same time, the problem is solvable by the method of characteristics, see Appendix 2 for a complete description. The exact solution is given by where p = −S(− p; n 2 , n 1 , −ν) with S defined in (22a). In (37), the first statement is free propagation, the second is refraction, the third comes from total internal reflection and the fourth statement comes from the compact support of ρ 0 . We introduce the critical momentum p c = n 2 1 − n 2 2 = 0.9798, where with p < p c total internal reflection occurs. This condition  is perhaps counter-intuitive, however this is due to the choice of optical axis. We furthermore define the following quantity, A ray which passes through the point (q, p) at z will hit the interface at z − Δz(q, p). The variable a is simply the propagation speed on the screen. We see that, although in principle not too difficult to solve, the expressions become large and unwieldy. Going to a slightly more complicated geometry already precludes the solvability by hand. We take as initial condition the distribution see Fig. 7. These initial conditions contain the critical momentum p c , such that both refraction and total internal reflection will occur. The exact solution to the Bucket of Water problem is plotted in Fig. 8. We use an integration length of z = 0.4 with 1000 steps on an 800 × 800 grid. For the discretization of phase space we use 800 grid points for both position and momentum, see Fig. 9. We see from Fig. 8 and 9 that the numerical solution is very close to the exact solution. The only difference that is noticeable by eye is the numerical diffusion at the edges. However, note that the edge at the interface is sharp. The numerical diffusion is of course an effect of the first order accurate upwind scheme. The sharp edges are an effect of the discrete symplectic transformation which occurs at the interface. In fact, the scheme is entirely symplectic in principle, the diffusion is generated by the linear interpolation that is inherent in the upwind scheme.
In this case, the Bucket of Water problem, we can implement the scheme in a matrixvector formulation, where the evolution matrix is fixed. Hence, we only have to compute this evolution matrix once, and we can then use it for any initial condition. Every step in z then becomes a matrix-vector product, where the matrix is very sparse. This will not be the case in our next example, where the evolution matrix has to be recomputed every time step. However, we shall first investigate numerically the convergence properties of the scheme.

Convergence Results
The convergence of our scheme is tested numerically by performing several runs of the Bucket of Water problem with different grid sizes. The refractive index is once again given by (36). However, as our analysis of the error is only valid for functions that are smooth away from the interface, we will use a smooth initial condition, ρ 0 ∈ C ∞ 0 (P). This initial condition is constructed by using the "bump" function [44], given by One can check that the bump function is infinitely smooth and has compact support. The initial condition is constructed as follows, where which results in the support of ρ 0 being 1 10 ≤ q ≤ 2 5 , 1 10 ≤ p ≤ 11 10 . The integration time is chosen as z = 2 5 , which results in a large part of the initial condition to be refracted and reflected by the interface.
We compute the error by using the exact solution and taking the L 1 -norm of the difference, i.e., where z T = 2 5 . For simplicity, we shall choose Δq = Δp, the results are displayed in Table 1. As the analytical solution for this problem is known, we can use it as a benchmark test and determine the error of our solver exactly. The results clearly show that our scheme is first order accurate in the number of grid points in the q-direction. This is as expected since an ordinary upwind solver has an error of O(Δq) + O(Δp). The adjustments to the upwind scheme have virtually no impact on the error behaviour.

Compound Parabolic Concentrator
The compound parabolic concentrator (CPC) is a pair of mirrors that have a parabolic shape. The CPC represents a worst-case scenario for Liouville's equation, as many different parts of phase space are interacting. Furthermore, the curved mirrors force us to recompute the evolution matrix at each time step. There are rays present that travel perpendicular to the optical axis that have a significant influence on the output, which causes a severe time-step restriction. These issues are not present when one uses a ray tracing method. Moreover, the CPC allows analytical computation of the intersection point of a ray with the mirrors. Thus, the CPC presents many issues for our scheme while presenting many advantages to ray tracing methods.
We assume that the CPC is embedded in a medium of unit refractive index, say air. In two dimensions, the CPC has ideal transmission characteristics, meaning all incoming light within an acceptance angle θ at the entrance aperture is captured and concentrated onto the exit aperture [45]. The exit aperture is represented by the subset (−a, a) × (−1, 1) ⊂ P. The CPC can be constructed by tilting two parabolas, one over angle θ and one over angle −θ , and shifting them such that their focal points are at (−a, 0) and (a, 0), respectively. Finally, one has to pick the focal point distance such that the parabolae go through the points (a, 0) and (−a, 0), respectively, see Fig. 10.
After rotating and shifting a standard parabola and setting the focal point correctly, one wall of the CPC can be characterized by the equation, We can solve (44) for q and select the positive root, which we define as Q r : R + → R. For positive q, we use q = Q r (z) as the curve defining the wall of the CPC, for negative q we use −q = Q r (z). The normal can be found easily, since Q r is differentiable. The shape of a CPC with exit aperture 2a and acceptance angle θ is completely fixed. The optic has a length Z , given by We can apply Liouville's Theorem to find how area in phase space is transformed when traversing the CPC. The CPC accepts momenta in the range (− sin θ, sin θ), while the spatial aperture is given by (46), thus the total area in phase space representing the entrance aperture is 4a. The exit aperture has width 2a, so the range of momenta at the exit aperture must be (−1, 1), which corresponds to an angular range of (− π 2 , π 2 ). Hence, the CPC provides maximal concentration in the spatial coordinate, while the price to pay is maximum dilution in the angular coordinate.
This also provides us with a special case that allows us to construct an exact solution to Liouville's equation. If the solution is given by the characteristic function of (−a , a ) × (− sin θ, sin θ) at the entrance aperture, the solution at the exit aperture must be the characteristic function of (−a, a) × (−1, 1). Conversely, using as an initial condition the solution at z = Z is then given by The results are plotted in Figs. 11 and 12. We compute the numerical solution using a 300×300 grid and 500 time steps. We have used θ = π 6 and a = 1 2 . We have integrated Liouville's equation on the CPC to z = Z given by (45).
The figure shows a great resemblance between the numerical and exact solutions. Besides from some numerical diffusion at the top and bottom edges, the numerical solution is equal to the exact solution. An intermediate result at z = Z /3 is shown in Fig. 13. It shows that the numerical diffusion comes from the upwind scheme, whereas the sharp edges come from the symplectic transformation that is the reflection.
We shall now compare our scheme to Monte Carlo ray tracing, which is used in conjunction with a scaled histogram to represent the solution. The histogram is scaled with the average number of rays per grid point such that the average value of the scaled histogram is equal to unity. By the Central Limit Theorem, we see that this procedure converges with the square root of the average number of rays per box, i.e., e RT ∼ M N RT , where M is the number of boxes. In our comparison, we use the control volume of a grid point as the box, with N p = N q leading to N 2 q boxes and thus e RT ∼ N q √ N RT . However, it is important to note that even when It can be shown that the minimum error is made for N RT ∼ N 4 q . The Liouville solver is again implemented using a matrix-vector formulation. We first construct a basic evolution matrix, where every time step needs only 2N q elements adjusted from the basic evolution matrix. The exact solution is only known at z = Z , hence that is where we compare both methods against the exact solution. We choose N p = N q and we fix the CFL number for the q-direction to be 0.96. Note that for a given initial position of a ray one can analytically determine the intersection points with the CPC. Thus, this test case is certainly favouring the ray tracing method, as more general optics would require a root finding method to find all intersection points, adding to the computation time. Ray tracers may also generate errors by not finding the right intersection point.
Next, we will vary the number of grid points, and therefore the boxes, but we keep the average number of rays per box constant. We point out that typical choices for the grid size in applications are 300 × 300, hence we have chosen our comparison to include this particular value. We have fixed the average number of rays per grid point at 100, thus choosing N RT = 100N 2 q , the results are displayed in Table 2. For all grid sizes, our method is faster, in terms of computation time, than Monte Carlo ray tracing. One does notice that the scaling behaviour is different, ray tracing having a linear time scaling behaviour in the number of rays, t RT ∼ N RT and our method having quadratic time scaling in the total number of grid points t L ∼ (N q N p ) 2 , which in our case reduces to t L ∼ N 4 q . However, this is the same scaling as we have found earlier to produce a minimal error for the ray tracer.
Let us fix N q and assume e RT ∼ 1 √ N RT , thus ignoring the error due to the box count. Again, the ray tracer takes a computation time of t RT ∼ N RT , allowing an estimate of the computation time to achieve the desired error tolerance. Suppose we want to reduce the error by a factor of c, then we have to increase the number of rays by a factor of 1 c 2 . Consider the case N RT = 10 6 from Table 2, and let us estimate the necessary number of rays to obtain the same global error as by solving Liouville's equation with N q = 100. We find the scaling constant as c = 8.31·10 −2 0.1394 = 0.5961 and 1 c = 1.6775 leading to 2.8 · 10 6 rays needed. The time would then also multiply by a factor of 1 c 2 , leading to a computation time of 1 min 50 s. Rescaling the other computation times to obtain the same error as the corresponding Liouville solver cases leads to estimated times of 11 min 47 s, 33 min 9 s and 1 h 9 min 55 s respectively. These estimated times are also presented in Table 2 ast RT . Note, however, that the assumption of a negligible error due to the box count results in these estimates being lower bounds. The fact that these estimates are lower bounds can be demonstrated by numerically investigating the convergence behaviour of the Monte Carlo ray tracing method. We shall fix N q = 100 and run the ray tracer for 1,4,9 and 16 million rays, the results are displayed in Table 3. The table clearly shows that only for 16 million rays is the ray tracing method really approaching the same global error as our scheme. Since the computation time for ray tracing only depends on the number of rays, we may read the computation times from Table 2. Thus, in reality, using ray tracing to obtain a solution that has the same global error as our scheme will take a much longer time than the estimates from Table 2.
Hence, solving Liouville's equation appears to be much more efficient in terms of computing time when a specific error tolerance has to be met. We conclude that when one is interested in obtaining the full phase space description, solving Liouville's equation is certainly a more efficient approach compared to Monte Carlo ray tracing. There are increasingly many applications where full information on phase space is desirable, e.g., the optical mixing of light and the studying of aberrations in free-form optics [9,10].

Conclusions
We have constructed a numerical method which is able to obtain the correct physical solution to Liouville's equation when the Hamiltonian is discontinuous. We have done this by building the correct physics into the numerical scheme. When there are no interfaces, i.e., discontinuities in the Hamiltonian, the scheme is simply the upwind scheme. When encountering an interface, the upwind grid point is selected using Snell's law.
Furthermore, we have derived a general jump condition assuming the initial conditions are smooth. This jump condition is what allowed us to derive a consistent scheme. In the case of a simple jump such as the Bucket of Water, the resulting advection speed is the harmonic average of the speeds on both sides of the interface.
For very simple geometries, we can apply the method of characteristics (MOC) to find a global solution analytically. We find excellent agreement between the exact and numerical solutions. The only difference being the numerical diffusion that is inherent in first order methods for hyperbolic PDEs. We were also able to apply our scheme to the compound parabolic concentrator (CPC) with great success. Liouville's equation has no global analytical solution for the CPC. However, in an important special case, we can find the solution at the entrance and exit apertures, based on Liouville's Theorem. Apart from some numerical diffusion around p = ± sin θ , where θ is the acceptance angle, the numerical solution is equal to the exact solution. We have also shown by these two examples that our solver for Liouville's equation gives an approach that is likely to be more efficient, in terms of computation time, than Monte Carlo ray tracing. Especially when a certain error tolerance is to be met, our approach is certainly faster.
We intend to extend our scheme to a three-dimensional optical setting. A naive approach of a uniform grid quickly grows to be computationally infeasible. Whereas phase space in a two-dimensional optical setting is two-dimensional, for a three-dimensional setting phase space becomes four-dimensional. Hence, a uniform grid will have N 4 grid points as opposed to the N 2 in the two-dimensional case. Our research effort will therefore be focused on nonuniform grids and hopefully we may drastically reduce the number of grid points needed.
In future work, we might also extend our scheme to include the Fresnel coefficients at sharp interfaces. When Fresnel reflection is taken into account, we must adjust (5) accordingly. Fresnel's equations specify exactly how much light is reflected and transmitted and thus, we will be able to find a numerical solution to Liouville's equation. We will also try to obtain higher-order accuracy by using more sophisticated methods such as high resolution schemes [46], (W)ENO reconstruction [47,48] and specialized time integrators [49][50][51]. We now observe that, from (54), we have d dz ρ(z, q(z), p(z)) = 0, whenever h is smooth. Integrating this relation, we find ρ (z, q(z), p(z)) = ρ 0 (q(0), p(0)) .
However, (58) gives us the physically correct solution, also when the Hamiltonian h is not smooth, provided we can find the starting point (q(0), p(0)) of any given (q(z), p(z)) = (q, p). We assume the initial condition, the function ρ 0 , has compact support in {q < 0, p ≥ 0}, which gives rise to four regions in phase space. First, we have the source region {q < 0, p ≥ 0}, the only region where there is light at z = 0. Second, we have the refracted region {q ≥ 0, p ≥ 0}, which is to say the region of phase space where light from the source region can only come by refraction. Third, there is the reflected region {q < 0, p < 0}, where light can only come due to reflection. Last, there is the empty region {q > 0, p < 0}, as there is no way for light to come from the source to this region. We can treat each region separately.

Source Region
The source region is the only region where ρ 0 is non-zero at z = 0. It is the collection of rays that have not encountered the interface yet. The rays inside this region are given by q(z) = q(0) + z p(0) We can thus find the solution in the source region by solving for (q(0), p(0)) and applying (58).

Refracted region
We will now treat the refracted region, where we must immediately introduce the critical momentum p c . The critical momentum satisfies δ = 0 in (7b). Solving for ψ, we obtain ψ c = − n 2 1 − n 2 2 , where ν = −1 leading us to p c = n 2 1 − n 2 2 .
Due to the choice of the optical axis, we have that momenta smaller than the critical momentum suffer total internal reflection. Hence, when the initial momentum exceeds the critical momentum, rays refract into the region {q ≥ 0, p ≥ 0}, see Fig. 14. (62b) We shall first find Δz, which is the travel time from the interface to the point (q, p) in the refracted region. Thus, a ray from the refracted region satisfies q(z − Δz) = 0. Therefore, Δz is given as Next, we integrate forward in time a characteristic that starts at (q(0), p ), where p is given by Snell's law in reverse, Corollary 1. However, this being the same ray, we have that it must satisfy q(z − Δz) = 0, giving the following equation for q(0), Furthermore, we have of course that p(0) = p . Thus, this completely expresses q(0) and p(0) in terms of z, q and p.

Reflected region
We now turn to the reflected region, which is the set of phase space coordinates that can only be reached by reflection, see Fig. 15. Since n 2 > 0, the only way a characteristic can reach the reflected region from the source region is by TIR. Therefore only the characteristics which have a momentum smaller than the critical momentum p c from (61) will reflect. The strategy is again to integrate the Hamilton equation backward in z, with the condition that where now (q, p) is in the reflected region. Again, the travel time Δz from (q, p) to the interface is given by (63), and q(z − Δz) = 0. Now we can apply (64), however, the reflected momentum is now explicitly given by p(0) = −p.
Thus, we find that (64) becomes the equation