Jocelyn Meyron

Quentin Mérigot

Boris Thibert

Light in Power: A General and Parameter-free Algorithm for Caustic Design

We present in this paper a generic and parameter-free algorithm to efficiently build a wide variety of optical components, such as mirrors or lenses, that satisfy some light energy constraints. In all of our problems, one is given a collimated or point light source and a desired illumination after reflection or refraction and the goal is to design the geometry of a mirror or lens which transports exactly the light emitted by the source onto the target. We first propose a general framework and show that eight different optical component design problems amount to solving a light energy conservation equation that involves the computation of visibility diagrams. We then show that these diagrams all have the same structure and can be obtained by intersecting a 3D Power diagram with a planar or spherical domain. This allows us to propose an efficient and fully generic algorithm capable to solve these eight optical component design problems. The support of the prescribed target illumination can be a set of directions or a set of points located at a finite distance. Our solutions satisfy design constraints such as convexity or concavity. We show the effectiveness of our algorithm on simulated and fabricated examples.

: Our algorithm can be used to design mirrors and lenses that reflect or refract collimated or point light sources onto a prescribed distribution of light. From left to right: Three lenses that refract the three channels of a color image; Mirror that reflects a point light source (located inside the mirror); Fabricated mirror that reflects a collimated light source; Fabricated lens that refracts a collimated light source.

Introduction

The field of non-imaging optics deals with the design of optical components whose goal is to transfer the radiation emitted by a light source onto a prescribed target. This question is at the heart of many applications where one wants to optimize the use of light energy by decreasing light loss or light pollution. Such problems appear in the design of car beams, public lighting, solar ovens and hydroponic agriculture. This problem has also been considered under the name of caustic design, with applications in architecture and interior decoration [START_REF] Finckh | Geometry construction from caustic images[END_REF].

In this paper, we consider the problem of designing a wide variety of mirrors and lenses that satisfy different kinds of light energy constraints. To be a little bit more specific, in each problem that we consider, one is given a light source and a desired illumination after reflection or refraction which is called the target. The goal is to design the geometry of a mirror or lens which transports exactly the light emitted by the source onto the target. The design of such optical components can be thought of as an inverse problem, where the forward problem would be the simulation of the target illumination from the description of the light source and the geometry of the mirror or lens.

In practice, the mirror or lens needs to satisfy aesthetic and pragmatic design constraints. In many situations, such as for the construction of car lights, physical molds are built by milling and the mirror or lens is built on this mold. Sometimes the optical component itself is directly milled. This imposes some constraints that can be achieved by imposing convexity or smoothness conditions. The convexity constraint is classical since it allows in particular to mill the component with a tool of arbitrary large radius. Conversely, concavity allows to mill its mold. Also, convex mirrors are easier to chrome-plate, because convex surfaces have no bumps in which the chrome would spuriously concentrate [START_REF] Cork | Method of and mixture for aluminizing a metal surface[END_REF].

In this paper, we propose a generic algorithm capable of solving eight different caustic design problems, see Figures 1 and2. Our approach relies on the relation between these problems and optimal transport. The algorithm is fully generic in the sense that it can deal with any of the eight caustic-design problems just by changing a formula, and can handle virtually any ideal light source and target. Our contributions are the following:

• We propose a general framework for eight different optical component design problems (i.e. four non-imaging problems, for which we can produce either concave or convex solutions). These problems amount to solving the same light energy conservation equation (see Sec. 3), which involves prescribing the amount of light reflected or refracted in a finite number of directions. • We propose a single algorithm with no parameter capable to solve this equation for the eight different problems.

We will see that, in order to solve this equation, we need to compute integrals over visibility cells, which can be obtained in all cases by intersecting a 3D Power diagram with a planar of spherical domain. The equation is then solved using a damped Newton algorithm. • In all of the four non-imaging problems, we can construct either a concave or convex optical component, easing their fabrication. Several components can then be combined to produce a single caustic, providing resilience to small obstacles and providing degrees of freedom to to control the shape of the optical system. • We show that we can solve near-field problems (when the target is at a finite distance) by iteratively solving farfield problems (when the target is a set of directions) for the eight optical component design problems.

Related work

The field of non-imaging optics has been extensively studied in the last thirty years. We give below an overview of the main In each case, the goal is to build the surface R of a mirror or a lens. For each problem, we provide two solutions (for instance, convex or concave surfaces when the light source is collimated). Top/Bottom: Collimated light sources/Point light sources. Left/Right: Mirror/Lens design. approaches to tackle several of the problems of this field. A survey on inverse surface design from light transport behaviour is provided by [START_REF] Patow | A survey of inverse surface design from light transport behavior specification[END_REF] Energy minimization methods. Many different methods to solve inverse problems arising in non-imaging optics rely on variational approaches. When the energies to be minimized are not convex, they can be handled by different kind of iterative methods. One class of methods deals with stochastic optimization. [START_REF] Finckh | Geometry construction from caustic images[END_REF] propose to represent the optical component (mirror or lens) as a C 2 B-spline triangle mesh and to use stochastic optimization to adjust the heights of the vertices so as to minimize a light energy constraint. Note that this approach is very costly, since a forward simulation needs to be done at every step and the number of steps is very high in practice. Furthermore, using this method, lots of artifacts in the final caustic images are present. Stochastic optimization has also been used by [START_REF] Papas | Goal-based caustics[END_REF] to design reflective or refractive caustics for collimated light sources. At the center of the method is the Expectation Minimization algorithm initialized with a Capacity Constrained Voronoi Tessellation (CCVT) using a variant of Lloyd's algorithm [Lloyd, 1982]. The source is a uniform directional light and is modeled using an array of curved microfacets. The target is represented by a mixture of Gaussian kernel functions. This method cannot accurately handle low intensity regions and artifacts due to the discretization are present. Microfacets were also used by [START_REF] Weyrich | Fabricating microgeometry for custom surface reflectance[END_REF] to represent the mirror. Due to the sampling procedure, this method cannot correctly handle smooth regions and does not scale well with the size of the target. More recently, [START_REF] Piovarči | Directional screens[END_REF] used microgeometry to design directional screens which provide increased gain and brightness. In their approach, the screen is decomposed into many small patches, each patch reflecting a set of rays toward a prescribed cone of directions. Their problem for each patch is similar to the one we solve for the special case of a directional source and a target at infinity and corresponds to one pillow, see Section 7. Similarly to us, their approach is based on convex optimization and produce convex patches. They prescribe the areas of the facets, whereas we prescribe their measures. Their numerical approach relies on gradient descent rather than on Newton's method, and only deals with collimated source and far-field target.

The approaches proposed by [START_REF] Kiser | Architectural caustic-sâĂŤcontrolling light with geometry[END_REF], Yue et al., 2014and Schwartzburg et al., 2014 have in common that they first compute some kind of relationship between the incident rays and their position on the target screen and then use an iterative method to compute the shape of the refractive surface. [START_REF] Yue | Poisson-based continuous surface generation for goal-based caustics[END_REF] use a continuous parametrization and thus cannot correctly handle totally black and high-contrast regions (boundaries between very dark and very bright areas). [START_REF] Yue | Pixel art with refracted light by rearrangeable sticks[END_REF] proposed to use sticks to represent the refractive surface. This allows to reduce production cost, to be more entertaining for the user since a single set of sticks can produce different caustic patterns. The main problem with this approach is the computational complexity since they need to solve a NP-hard assignment problem. The problem of designing lenses for collimated light sources has also been considered by [START_REF] Schwartzburg | High-contrast computational caustic design[END_REF] They propose a method to build lenses that can refract complicated and highly contrasted targets. They first use optimal transport on the target space to compute a mapping between the refracted rays of an initial lens and the desired normals, then perform a post-processing step to build a surface whose normals are close to the desired ones.

Monge-Ampère equations When the source and target lights are modeled by continuous functions, the problem amounts to solving a generalized Monge-Ampère equation, either in the plane for collimated light sources, or on the sphere for point light sources. Let us explain this link more precisely for a collimated light source, assuming that the source rays are collinear to the constant vector e z = (0, 0, 1) and emitted from a horizontal domain Ω ⊂ R 2 × {0}. We assume that the optical component is smooth and parameterized by a heightfield function ϕ : Ω → R and denote by µ and ν the source and target measures. At every point ϕ(x) of the optical component, the gradient ∇ϕ(x) encodes the normal, and we denote by F (∇ϕ(x)) ∈ S 2 the direction of the ray that is reflected at (x, ϕ(x)) using Snell's law. The conservation of light energy thus reads ν(A) = µ((F • ∇ϕ) -1 (A)) for every set A. This is equivalent to having ν(B) = µ((∇ϕ) -1 (B)) for every set B, where ν(B) = ν(F (B)). When F and ∇ϕ are one-toone (which is the case if the optical component is convex or concave) and µ and ν are modeled by continuous functions f and g, with the change of variable formula, the light energy conservation becomes equivalent to the following generalized Monge-Ampère equation

g(∇ϕ(x)) det(D 2 ϕ(x)) = f (x).
(1)

Similar equations are obtained for point light sources. The existence and regularity of their solutions, namely of the mirror or lens surfaces, have been extensively studied. When the light source is a point, this problem has been studied for mirrors [START_REF] Caffarelli | Weak solutions of one inverse problem in geometric optics[END_REF]Caffarelli et al., 2008] and lenses [START_REF] Gutiérrez | The refractor problem in reshaping light beams[END_REF] and when the light source is collimated one recovers (1) [START_REF] Gutiérrez | The parallel refractor[END_REF].

We refer to the book of Gutiérrez, 2016 for an introduction to Monge-Ampère equations.

Optimal transport based methods in non-imaging optics In fact, the Monge-Ampère equations corresponding to the non-imaging problems considered in this paper can be recast as optimal transport problems. This was first observed by Wang, 2004 and[START_REF] Glimm | [END_REF][START_REF] Glimm | [END_REF] for the mirror problem with a point light source. Many algorithms related to optimal transport have been developed to address non-imaging problems. For collimated sources, one could rely on widestencils finite difference schemes [START_REF] Prins | A numerical method for the design of free-form reflectors for lighting applications[END_REF], or on numerical solvers for quadratic semi-discrete optimal transport [Mérigot, 2011;de Goes et al., 2012]. For point sources, there exist variants of the Oliker-Prussner algorithm for the mirror problem [START_REF] Caffarelli | On the numerical solution of the problem of reflector design with given far-field scattering data[END_REF] or the lens problem [START_REF] Gutiérrez | The refractor problem in reshaping light beams[END_REF]. Both algorithms have a O(N 4 ) complexity, restricting their use to small discretizations. A quasi-Newton method is proposed by de Castro et al., 2016 for point-source reflector design, handling up to 10 4 Dirac masses.

Finally we note that the approach of [START_REF] Schwartzburg | High-contrast computational caustic design[END_REF] to build lenses also relies on optimal transport. However, the optimal transport step is used as a heuristic to estimate the normals of the surface, and not to directly construct a solution to the non-imaging problem. A post-processing step is then performed by minimizing a non-convex energy composed of five weighted terms. In contrast, all the results presented in this article use no post-processing.

Light energy conservation

We present in this section several mirror and lens design problems arising in non-imaging optics. In all the problems, one is given a light source (emitted by either a plane or a point) and a desired illumination "at infinity" after reflection or refraction, which is called the target, and the goal is to design the geometry of a mirror or lens which transports the energy emitted by the source onto the target. We do not take into account multiple reflections or refractions. We show that even though the problems we consider are quite different from one another, they share a common structure that corresponds to a so-called generalized Monge-Ampère equation, whose discrete version is given by Equation (2). This section gathers and reformulates in a unified setting results about mirror and lens design for collimated and point light sources. We refer to the work of [START_REF] Caffarelli | Weak solutions of one inverse problem in geometric optics[END_REF], Prins et al., 2013, de Castro et al., 2016 and references therein.

Collimated light source

Convex mirror design

In this first problem, the light source is collimated, meaning that it emits parallel rays. The source is encoded by a light intensity function ρ : Ω → R over a 2D domain Ω contained in the (xy) plane R 2 × {0} ⊂ R 3 , and that all the rays are parallel and directed towards e z = (0, 0, 1). For simplicity, we will conflate R 2 and R 2 × {0}. The desired target illumination is "at infinity", and is described by a set of intensity values σ = (σ i ) 1 i N supported on a finite set of directions Y = {y 1 , • • • , y N } included in the unit sphere S 2 . The problem is to find the surface R of a mirror that reflects the source intensity ρ to the target intensity σ, see Figures 2 (top left) and 3.

Since the number of reflected directions Y is finite, the mirror surface R is composed of a finite number of planar facets as illustrated in Figure 3. We will construct the mirror R as the graph of a piecewise-linear convex function of the form x ∈ R 2 → max i x|p i -ψ i , where •|• denotes the dot product. The vectors p 1 , . . . , p N ∈ R 2 and the elevations ψ = (ψ 1 , . . . , ψ N ) ∈ R N have to be determined. To choose p i , we require that the plane

P i = {(x, x|p i ) | x ∈ R 2 }
reflects rays with direction e z towards the direction y i ∈ S 2 . The downward pointing unit normal to P i is given by the formula ). The mirror surface R is the graph of a convex piecewise affine function. The support X ρ of ρ is decomposed into visibility cells (V i (ψ)) 1 i N . Every vertical ray above a point x ∈ X ρ belongs to a cell V i (ψ), touches a plane with slope p i and is reflected to the direction y i .

n i = (p i , -1)/( p i 2 + 1) ∈ R 3
, and Snell's law of reflection gives us y i = e z -2 n i |e z n i . Solving for p i and denoting p R 2 the orthogonal projection onto R 2 × {0}, we get

p i = -p R 2 (y i -e z )/ y i -e z |e z .

Given a vector of elevations

ψ := (ψ i ) 1 i N , we define the visibility cell V i (ψ) V i (ψ) = {x ∈ R 2 × {0} | ∀j, -x|p i + ψ i -x|p j + ψ j }.
By construction, for every i ∈ {1, . . . , N }, any vertical ray emanating from a point x ∈ V i (ψ) touches the mirror surface R at an altitude x|p i -ψ i , and is thus reflected towards direction y i . Consequently, the amount of light reflected towards direction y i equals the integral of ρ over V i (ψ). The Collimated Source Mirror problem (CS/Mirror) then amounts to finding elevations

ψ ∈ R N such that ∀i ∈ {1, • • • , n} Vi(ψ) ρ(x)x . = σ i .
(2)

By construction, a solution to Equation ( 2) provides a parameterization R ψ of a convex mirror that reflects the collimated light source ρ to the discrete target σ:

R ψ : x ∈ R 2 → (x, max i x|p i -ψ i ) ∈ R 3 ,
where R 2 × {0} and R 2 are identified. Notice that since the mirror is a graph over R 2 × {0}, the vectors y i cannot be upward vertical. In practice we assume that every direction y i belongs to the hemisphere S 2 -:= {x ∈ S 2 , x|e z 0}. Furthermore, we localize the position of the mirror by considering it only above the domain

X ρ := {x ∈ R 2 × {0}, ρ(x) = 0}.
Concave mirror. The same approach also allows the construction of concave mirrors, using a concave function of the form x → min i x|p i + ψ i . This amounts to replacing the visibility cells by

V i (ψ) = {x ∈ R 2 × {0} | x|p i + ψ i x|p j + ψ j ∀j}.
In that case, a solution to Equation ( 2) provides a parametrization of a concave mirror R ψ (x) = (x, min i x|p i + ψ i ) that sends the collimated light source ρ to the discrete target σ.

Convex lens design

In this second design problem, we are interested in designing lenses that refract a given collimated light source intensity ρ to a target intensity σ, see the top right diagram in Figure 2. We denote by n 1 the refractive index of the lens, by n 2 the ambient space refractive index and by κ = n 1 /n 2 the ratio of the two indices. We assume that the rays emitted by the source are vertical and that the bottom of the lens is flat and orthogonal to the vertical axis. There is no refraction angle when the rays enter the lens, and we only need to build the top part of the lens. By a simple change of variable, we show that this problem is equivalent to (CS/Mirror). More precisely, for every y i ∈ Y , we now define p i to be the slope of a plane that refracts the vertical ray e z to the direction y i . We define R as the graph of a convex function of the form x → max i x|p i -ψ i , where ψ = (ψ i ) 1 i N is the set of elevations. A calculation similar to the (CS/Mirror) case gives the following expression:

p i = -p R 2 (y i -κe z )/ y i -κe z |e z .
In that case, we define the visibility cell V i (ψ) to be the set of points x ∈ R 2 × {0} that are refracted to the direction y i :

V i (ψ) = {x ∈ R 2 × {0} | ∀j, -x|p i + ψ i -x|p j + ψ j }.
The Collimated Source Lens problem (CS/Lens) then amounts to finding weights (ψ i ) 1 i N that satisfy (2). In that case the lens surface is parameterized by

R ψ : x ∈ R 2 → (x, max i x|p i -ψ i ).
In practice, we choose the directions y i in S 2 + and the mirror to be parameterized over the support X ρ of ρ.

Concave lens. Note that we can also build concave lenses by considering parameterizations with concave functions of the form x → min i x|p i + ψ i . Figure 4 illustrates a concave and a convex solution to the same non-imaging optics problem. In this second mirror design problem, all the rays are now emitted from a single point in space, located at the origin, and the light source is described by an intensity function ρ on the unit sphere S 2 . As in the previous cases, the target is "at infinity" and is described by a set of intensity values

σ = (σ i ) 1 i N supported on the finite set of directions Y = {y 1 , • • • , y N } ⊂ S 2 .
The problem we consider is to find the surface R of a mirror that sends the light intensity ρ to the light intensity σ (Fig. 2, bottom left).

Following [START_REF] Caffarelli | Weak solutions of one inverse problem in geometric optics[END_REF] , we build a concave surface R that is composed of pieces of confocal paraboloids. More precisely, we denote by P (y i , ψ i ) the solid (i.e filled) paraboloid whose focal point is at the origin with focal distance ψ i and with direction y i . We define the surface R ψ as the boundary of the intersection of the solid paraboloids, namely R ψ = ∂ (∩ i P (y i , ψ i )). The visibility cell V i (ψ) is the set of ray directions x ∈ S 2 emanating from the light source that are reflected in the direction y i . Since each paraboloid ∂P (y i , ψ i ) is parameterized over the sphere by x → ψ i x/(1 -x|y i ), one has [START_REF] Caffarelli | Weak solutions of one inverse problem in geometric optics[END_REF]]

V i (ψ) = x ∈ S 2 | ∀j, ψ i 1 -x|y i ψ j 1 -x|y j .
The Point Source Mirror problem (PS/Mirror) then amounts to finding (ψ i ) that satisfy the light energy conservation equation (2). The mirror surface is then parameterized by

R ψ : x ∈ S 2 → min i ψ i 1 -x|y i x.
In practice, we assume that the target Y is included in S 2 -, that the support X ρ of ρ is included S 2 + := {x ∈ S 2 , x|e z 0}, and that the mirror is parameterized over X ρ .

One can also define the mirror surface as the boundary of the union (instead of the intersection) of a family of solid paraboloids. Then, the visibility cells become

V i (ψ) = x ∈ S 2 | ∀j, ψ i 1 -x|y i ψ j 1 -x|y j
and a solution to Equation (2) provides a parameterization R ψ (x) = x max i ψ i /(1 -x|y i ) of the mirror surface. Let us note that in this case the mirror is neither convex nor concave.

Convex lens design.

We now consider the lens design problem for a point light source. As in the collimated setting, we fix the bottom part of the lens. We choose a piece of sphere centered at the source, so that the rays are not deviated. Following [START_REF] Gutiérrez | The refractor problem in reshaping light beams[END_REF] , the lens is composed of pieces of ellipsoids of constant eccentricities κ > 1, where κ is the ratio of the indices of refraction. Each ellipsoid ∂E(y i , ψ i ) can be parameterized over the sphere by x → ψ i x/(1 -κ x|y i ). The visibility cell of y i is then

V i (ψ) = x ∈ S 2 | ∀j, ψ i 1 -κ x|y i ψ j 1 -κ x|y j .
The Point Source Lens problem (PS/Lens) then amounts to finding weights (ψ i ) 1 i N that satisfy (2). Note that the top surface of the lens is then parameterized by

R ψ : x ∈ S 2 → min i ψ i 1 -κ x|y i x.
In practice, we choose the set of directions y i to belong to S 2 + and the lens to be parameterized over the support

X ρ ⊂ S 2 + of ρ.
One can also choose to define the lens surface as the boundary of the union (instead of the intersection) of a family of solid ellipsoids. In that case, the visibility cells are given by

V i (ψ) = x ∈ S 2 | ∀j, ψ i 1 -κ x|y i ψ j 1 -κ x|y j
and a solution to Equation (2) provides a parameterization R ψ (x) = x max i ψ i /(1 -κ x|y i ) of the lens surface. Let us note that in this case the lens is neither convex nor concave.

General formulation

Let X be a domain of either the plane R 2 × {0} or the unit sphere S 2 , ρ : X → R a probability density and Y = {y 1 , • • • , y N } ⊂ S 2 be a set of N points. We define the function G : R N → R N by

G i (ψ) = Vi(ψ) ρ(x)x .
where G(ψ) = (G i (ψ)) 1 i N and V i (ψ) ⊂ X is the visibility cell of y i , whose definition depends on the non-imaging problem. Using this notation, Equation (2) can be rephrased as finding weights ψ = (ψ i )

1 i N such that ∀i ∈ {1, • • • , N }, G i (ψ) = σ i .
(3)

Many other problems arising in non-imaging optics amount to solving equations of this form. For example, the design of a lens that refracts a point light source to a desired nearfield target can also be modeled by a Monge-Ampère equation that has the same structure [START_REF] Gutiérrez | The refractor problem in reshaping light beams[END_REF]. In this case, the visibility diagram correspond to the radial projection onto the sphere of pieces of confocal ellipsoids with non constant eccentricities and is not associated to an optimal transport problem.

Visibility and Power cells

The main difficulty to evaluate the function G appearing in Equation ( 3) is to compute the visibility cells V i (ψ) associated to each optical modeling problem. We show in this section that the visibility cells have always the same structure, allowing us to build a generic algorithm in Section 5. We first need to introduce the notion of Power diagram.

Power diagrams. Let P ⊆ R 3 × R be a weighted point cloud, i.e. P = {(p i , ω i )} 1 i N with p i ∈ R 2 and ω i ∈ R. The Power cell of the ith point p i is given by

Pow i (P) := {x ∈ R 3 | ∀j, x -p i 2 + ω i x -p j 2 + ω j }.
Power cells partition R 3 into convex polyhedra up to a negligible set. Power diagrams are well-studied objects appearing in computational geometry [Aurenhammer, 1987], and can be computed efficiently in dimension 2 and 3. When all the weights are equal, one recover the usual Voronoi diagram.

Visibility diagram as a restricted Power diagram. We now show that in all the non-imaging problems of Section 3, the visibility cells are of the form

V i (ψ) = Pow i (P) ∩ X. (4)
For a collimated source, X denotes the plane R 2 × {0} and for a point source, X is the unit sphere S 2 . The expression of the weighted point cloud P = {(p i , ω i )} depends on the problem. We refer to

p i = - p R 2 (yi-ez) yi-ez|ez ω i = 2ψ i -p i 2 .
Ccv (CS/Mirror)

p i = p R 2 (yi-ez) yi-ez|ez ω i = 2ψ i -p i 2 .
Cvx (PS/Mirror) p i = -yi 2 ln(ψi)

ω i = -1 ln(ψi) - 1 4 ln(ψi) 2 Ccv (PS/Mirror) p i = y i /(2 ln(ψ i )) ω i = 1 ln(ψi)) - 1 4 ln(ψi) 2 Cvx (CS/Lens) p i = - p R 2 (yi-κez) yi-κez|ez ω i = 2ψ i -p i 2 Ccv (CS/Lens) p i = p R 2 (yi-κez) yi-κez|ez ω i = 2ψ i -p i 2 Cvx (PS/Lens) p i = -κ yi 2 ln(ψi) ω i = -1 ln(ψi) -κ 2 4 ln(ψi) 2
Ccv (PS/Lens)

p i = κy i /(2 ln(ψ i )) ω i = 1 ln(ψi) -κ 2 4 ln(ψi) 2
the formula in the (CS/Mirror) case, where the ith visibility cell is given by

V i (ψ) = {x ∈ R 2 × {0} | ∀j, -x|p i + ψ i -x|p j + ψ j } = {x ∈ R 2 × {0} | ∀j, x -p i 2 + ω i x -p j 2 + ω j },
where ω i = 2ψ i -p i 2 . We conclude that the visibility cells for a convex mirror of the (CS/Mirror) problem are indeed given by (4), where the weighted point cloud is given by the first line of Table 1.

A generic algorithm

For each optical design problem, given a light source intensity function, a target light intensity function and a tolerance, Algorithm 1 outputs a triangulation of a mirror or a lens that satisfies the light energy conservation equation (2).

The main problem is to find weights ψ such that G(ψ) = σ (see Equation ( 3)). This is done using a damped Newton algorithm similar to recent algorithms that have been shown to have a quadratic local convergence rate for optimal transport problems [START_REF] Kitagawa | A newton algorithm for semi-discrete optimal transport[END_REF] or for Monge-Ampère equations in the plane [Mirebeau, 2015]. A key point of this algorithm is to enforce the Jacobian matrix DG(ψ) to always be of rank N -1. To this purpose, we need to enforce all along the process that ∀i ∈ {1, . . . , N }, G i (ψ) > 0.

(5)

Indeed, first note that since G is invariant under the addition of a constant, the kernel of DG(ψ) always contains the constant vector (1, . . . , 1). Now note that if we have G i (ψ) = 0, then ALGORITHM 1: Mirror / lens construction Input A light source intensity function ρin.

A target light intensity function σin.

A tolerance η > 0.

Output A triangulation RT of a mirror or lens.

Step 1 Initialization (Section 5.1)

T, ρ ← DISCRETIZATION_SOURCE(ρin) Y, σ ← DISCRETIZATION_TARGET(σin) ψ 0 ← INITIAL_WEIGHTS(Y )
Step 2 Solve Equation (3):

G(ψ) = σ (Section 5.2) ψ ← DAMPED_NEWTON(T, ρ, Y, σ, ψ 0 , η)
Step 3 Construct a triangulation RT of R (Section 5.3) RT ← SURFACE_CONSTRUCTION(ψ, R ψ ) the corresponding visibility cell V i (ψ) is empty, which implies that ∇G i (ψ) = 0 (the gradient being taken with respect to ψ). This is because the gradient of G i involves integral on the boundary ∂V i (ψ), as shown for instance by [START_REF] Kitagawa | A newton algorithm for semi-discrete optimal transport[END_REF] in Theorem 1.3. Hence, if G i (ψ) = 0, then the rank of DG(ψ) is at most N -2 which prevents from using the Damped Newton method. Our method consists of three steps, described in Algorithm 1:

• Initialization (Sec. 5.1): We first discretize the source density into a piecewise affine density and the target one into a finitely supported measure. Then, we construct initial weights ψ 0 satisfying ∀i, G i (ψ (0) ) > 0. • Damped Newton (Sec 5.2): We construct a sequence

ψ k following Algorithm 2 until G(ψ k ) -σ ∞ η.
The main difficulty here is to evaluate G(ψ k ) and DG(ψ k ). • Surface construction (Sec 5.3): Finally, we convert the solution ψ k ∈ R N into a triangulation. Depending on the non-imaging problem, this amounts to approximating an intersection (or union) of half-spaces (or solid paraboloids, or ellipsoids) by a triangulation.

Initialization

Discretization of light intensity functions Our framework allows to handle any kind of collimated or point light source or target light intensity functions. It can be for example any positive function on the plane or the sphere (depending on the problem) or a greyscale image, which we see as piecewise affine function. We first approach the support of the source density ρ by a triangulation T and assume that the density ρ : T → R + is affine on each triangle. We then normalize ρ by dividing it by the total integral T ρ(x)x . . Input The source ρ : T → R + and target σ = i σiδy i ; an initial vector ψ 0 and a tolerance η > 0.

Step 1 Transformation to an Optimal Transport problem

If X = R 2 × {0}, then ψ 0 = ψ 0 (and G = G). If X = S 2 , then ψ 0 = (ln(ψ 0 i )) 1 i N (and G = (Gi • exp) 1 i N ).
Step 2 Solve the equation:

G( ψ) = σ Initialization: ε0 := min mini Gi(ψ 0 ), mini σi > 0, k := 0. While G( ψ k ) -σ ∞ η -Compute d k = -D G( ψ k ) + ( G( ψ k ) -σ) -Find the smallest ∈ N s.t. ψ k, := ψ k + 2 -d k satisfies      min i Gi( ψ k, )) ε0 G(ψ k, ) -σ ∞ (1 -2 -( +1) ) G( ψ k ) -σ ∞ -Set ψ k+1 = ψ k + 2 -d k and k ← k + 1. Return ψ := ( ψ k i ) 1 i N if X = R 2 × {0} or ψ := (exp( ψ k i )) 1 i N if X = S 2 .
Similarly, the target light intensity function can also be any discrete probability measure. If the user provides an image, one can transform it into a discrete measure of the form σ = i σ i δ yi using Lloyd's algorithm or more simply by taking one Dirac mass per pixel. We do the latter in all experiments. The target measure is also normalized by dividing with the discrete integral i σ i . We need min i σ i > 0 for the damped Newton algorithm, but this is not a restriction: if σ i = 0, we simply remove the corresponding Dirac mass δ yi , thus ensuring that no light is sent towards y i .

Choice of the initial family of weights ψ 0 . As mentioned at the beginning of this section, we need to ensure that at each iteration all the visibility cells have non-empty interiors.

In particular, we need to choose a set of initial weights ψ 0 = (ψ 0 i ) 1 i N such that the initial visibility cells are not empty.

• For the collimated light sources cases (CS/Mirror) and (CS/Lens), we see that if we choose ψ 0 i = p i 2 /2 then ω i = 0, where p i is obtained using the formulas of the Section 4. Then, the visibility diagram becomes a Voronoi diagram, hence p i ∈ V i (ψ 0 ).

• For the Point Source Mirror (PS/Mirror) case, an easy calculation shows that if we choose ψ 0 i = 1, then -y i ∈ V i (ψ 0 ).

• For the Point Source Lens (PS/Lens) case, we can show that if we also choose ψ 0 i = 1, then y i ∈ V i (ψ 0 ). Note that the previous expressions for ψ 0 ensure that G i (ψ 0 ) = ρ(V i (ψ 0 )) > 0 only when the support X ρ of the light source is large enough. As an example in the (PS/Mirror) case, if -y i / ∈ X ρ , then we may have G i (ψ 0 ) = 0. To handle this difficulty, we use a linear interpolation between ρ and a constant density supported on a set that contains the (-y i )'s. This strategy also works for the (CS/Mirror), (PS/Lens) and (CS/Lens) cases.

Damped Newton algorithm

When the light source is collimated (i.e. X = R 2 × {0}), the problem is known to be an optimal transport problem in the plane for the quadratic cost, the function G is the gradient of a concave function, its Jacobian matrix DG is symmetric and DG 0. Moreover, if G i (ψ) > 0 for all i and if X ρ is connected, then the kernel of DG is spanned by ψ = cst. This ensures the convergence of the damped Newton algorithm [START_REF] Kitagawa | A newton algorithm for semi-discrete optimal transport[END_REF] presented as Algorithm 2, where A + denotes the pseudo-inverse of the matrix A. Practically, taking the pseudo-inverse of D G( ψ k ) guarantees that the mean of the ψ k remains constant. In practice, we remove a line and a column of the matrix to make it full rank.

When the light source is a point source, we make the change of variable ψ = ln(ψ) and G = G • exp, so that G(ψ) = σ if and only if G( ψ) = σ. This change of variable turns the optical component design problem into an optimal transport problem, ensuring that G is the gradient of a concave function and that D G is symmetric negative [de Castro et al., 2016], thus easily invertible. In the (PS/Mirror) problem with convex mirrors, the damped Newton algorithm is also provably converging [START_REF] Kitagawa | A newton algorithm for semi-discrete optimal transport[END_REF].

Computation of G and DG By Section 4, the visibility cells V i (ψ) can be computed by intersecting a certain 3D Power diagram with a triangulation T of the support X ρ of ρ. Such intersection can be computed using the algorithm developed by Lévy, 2015 . Then,

G i (ψ) = Vi(ψ)
ρ(x)dx can be computed using first order quadrature formulas. The computation of DG is done using forward-mode automatic differentiation, where we store the gradient of G i (ψ) as a sparse vector. Note that this works quite efficiently since all numbers that occur in the computation of G i (ψ) depend only on ALGORITHM 3: Optical component design for a NF target.

Input The source ρ : T → R + and target σ = N i=1 σiδz i ; an initial vector ψ 0 and two tolerances η, ηNF > 0.

Initialization ∀i, c 0 i = O While c k+1 -c k 1 /N > ηNF -Compute v k i = R ψ k (c k i ) -Set y k i = (zi -v k i )/ zi -v k i -Solve ψ k+1 ← SOLVE_FF(T, ρ, Y k , σ, η) -Update c k+1 i to be the centroid of Vi(ψ k+1 )
the values ψ j where j is such that (i, j) are neighbors in the visibility diagram, i.e. V i (ψ) ∩ V j (ψ) = ∅.

Linear system Since D G is sparse and symmetric negative, we solve the linear systems using preconditioned conjugate gradient.

Surface construction

In the last step of Algorithm 1, we build a triangulation of the mirror or lens surface. The input is a family of weights ψ solving Equation (3) and the parameterization function R ψ whose formula is given in Section 3 and depends on the eight different cases. We triangulate each visibility cell by taking the convex hull of the vertices of its boundary. A vertex of the triangulation will belong to at least one visibility cell. For each vertex, we can compute exactly the normal to the (continuous surface) using Snell's law since we know the incident ray and the corresponding reflected/refracted direction y i .

Finite-distance caustics

In this section, we show that we solve the eight optical component design problems when the target is at a finite distance. This setting is called near-field (NF) in contrast with the previous far-field (FF) case that deals with targets at infinity. This setting is interesting since in most applications, one wants the focused image to be at a finite distance and not at infinity. To be more precise, given one of the eight optical component design problems mentioned above with a target illumination σ = N i=1 σ i δ zi supported on a set of points Z = {z 1 , . . . , z N } ⊂ R 3 , we propose an algorithm that iteratively solves a FF problem, namely Equation ( 2), and converges to a solution of the NF problem.

Algorithm

The procedure consists in solving a sequence of FF problems that quickly converges to a solution of the NF one. Details can be found in Algorithm 3. In this algorithm, SOLVE_FF(T, ρ, Y, σ, η) denotes an algorithm that solves the FF problem between a source ρ : T → R + and a target σ = i σ i δ yi supported on Y ⊂ S 2 for a numerical error η. It can for instance be Step 2 of Algorithm 1. This algorithm is used to produce all the lenses and mirrors of the article, except for the first image of Figure 15.

Convergence analysis

It is clear that when a fixed point is reached in Algorithm 3, the corresponding weight vector ψ is a good approximation of the NF problem. In practice, the algorithm converges very quickly: in all our examples, after 6 iterations, we get an error η N F of 10 -6 , see Figure 5.

We have no guarantee on the convergence of the discretization. However, note that the set of reflected (or refracted) rays emanating from a visibility cell has a diameter proportional to the diameter of the visibility cell. This generates a blur whose diameter is also proportional to the diameter of the cells. In practice, we observe that these cells have a small diameter (see Figures 7 and10) and indeed the blur is unnoticeable. The convergence is illustrated in Figure 6 where we show LuxRender renderings at iteration 1 (which corresponds to the FF setting), at iteration 2 and after 6 steps (when the error c k+1 -c k 1 /N becomes smaller than η N F = 10 -6 ). In all the images, we look at the projection of the refracted illumination onto the target screen. One can see that the first one is distorted since it corresponds to the projection of an image supported on the sphere onto a planar screen. Starting from the second iteration of the algorithm, the image is not distorted anymore.

The convergence of this algorithm is also illustrated in Fig- ure 15. On the first image of Figure 15, the optical component is composed of three lenses that solve the FF problem. We observe that the LuxRender rendering creates three shifted copies of the same image. To be more precise, the translation between the 3 images is the same as the one between the lenses i.e. if each lens has a width of 1 then each image is also shifted by 1. In the second image, the nine lenses solve the NF problem for an image at a finite distance. The fact that the images do superimpose show that the NF problem is accurately solved. The quality of the first image of Figure 1 also assesses that the NF problem is solved accurately.

Performance

We report performance of algorithms 1 and 3 to solve the (CS/Mirror) case for the Train target on a laptop with an i7 CPU. Table 2 shows the running times in the FF setting, We underline that the number of iterations in the Newton step remains low: in all our examples it varies from 10 to 20. This means that the computational cost of the method is concentrated in the computation of the functional G, its Jacobian matrix DG and the resolution of the linear system. We believe that there is much room for improvement in the first two steps, by optimizing the computation of visibility cells, and by using an explicit computation of DG instead of relying on automatic differentiation.

In the NF setting, with the same configuration, the results are summarized in Table 3. The total running time is the one after k = 6 iterations. Note that the running time of the second iteration is greater than the first one. Indeed, in Algorithm 3, we saw that we use the weights found at one step to initialize the next one. Since the target directions change greatly between the first and second steps, these weights do not provide a good initialization. To fix that, we automatically apply a simple perturbation to obtain good initial weights. Starting from the third step, the running times decrease as the target directions do not move a lot at each iteration.

Results and Discussion

In this section, we present several numerical examples for the different problems previously described as well as some other applications. In the experiments, we take κ = 1.5. Unless stated otherwise, the light source is chosen uniform and the discretization of the target (number of Diracs N ) is equal to the size of the image. The stopping criterion to Newton's algorithm (Algorithm 2) is set to η = 10 -8 . Since the convergence of Algorithm 3 is always fast, we do not use here the stopping criteria η N F and stop the algorithm after k = 6 iterations in all the examples presented hereafter.

Evaluation strategy

The output of our algorithm is a triangulation equipped with a normal at each vertex. In all the simulations, we use the LuxRender rendering engine, with Bidirectional Path Tracing combined with a Sobol sampler and the Fresnel coefficient is not taken into account.

Results

Genericity. Our algorithm is able to solve eight different optical component design problems. We present for instance in Figure 7 four examples for which we display the visibility diagram of X ρ as well as the optical component (lens or mirror) above it, a mesh of the optical component and a forward simulation using LuxRender.

Then, for the examples of Figures 8,9, 11, and 12, we display the target distribution as an image; a mean curvature plot (blue represents low mean curvature and red high mean curvature) of the constructed mesh R T and a forward simulation using LuxRender.

High-contrast and complex target lights. We can handle any kind of target distribution. Figures 8 and9 shows several examples of mirror design for respectively a collimated and a point light source. Note that we are able to construct mirrors for smooth images such as the Train image as well as images with totally black areas (third and fourth rows). We are also able to handle target supported on non-convex sets such as the Hikari and Siggraph images. One can notice that since the area of the visibility cells are equal to the greyscale values of the image then the triangles have roughly the same size, implying that one can recognize the target image in the mesh of the surface, see Figure 7 for zooms on different meshes. The mean curvature plot shows the discontinuities in the surface which come from the black areas in the image. Figures 11 and12 show the same kind of results for the lens design problems (CS/Lens) and (PS/Lens).

Non-uniform light sources. Our algorithms can be used with non-uniform light sources. Below, we compare the meshes that are generated in the (CS/Lens) case when the source is either uniform or a Gaussian (Fig. 10 left). Because of the higher concentration of light, the details of the triangulation are more concentrated in the middle in the Gaussian case (middle) than in the uniform case (right).

Convex / concave optical components. As shown in Section 3, for each problem, one can choose between two different parameterizations. For instance, for the (CS/Lens) problem, one can build a lens which is either concave or convex, see Figure 4 for an illustration of these differences. Note that in the (PS/Lens) setting (which corresponds to the last row of Figure 7 and Figure 12), the light source is not supported on the full hemisphere S 2 + but instead on a smaller part of it. Indeed, choosing a smaller support for µ enforces that R T is a graph above the plane instead of the hemisphere and thus avoids potential inter-refractions. Furthermore, since we parametrize the lens as a union of ellipsoids, it is neither convex nor concave. As for all figures, we have performed no post-processing on R T in order to emphasize the benefit of designing convex or concave optical components (convexity is a form of regularity). One also observes that when the lens is rotated with respect to the light source (Figure 14 and first row of Figure 19), or when the target screen is not at the right distance (Figure 18), the image is deformed in a monotonic and regular way. We believe this is due to the monotonicity properties of optimal transport and to the convexity/concavity properties of the optical components.

Comparison with previous work. Figure 13 compares the state of the art results obtained by Schwartzburg et al., 2014 (second column) and the LuxRender renderings obtained by our method (third column) on two target distributions for the (CS/Lens) case with a collimated uniform light source in the near-field setting. Although the results are comparable, one can notice, in the second column, the presence of small artifacts between the black and white regions, for instance around the rings (notably in the center). The contrast is more accurate with our convex lenses. Application to pillows This problem consists in decomposing the optical component (mirror or lens) into several smaller optical components that are called pillows, as illustrated in Figure 1, and are widely used in car headlight design. Each pillow independently satisfies a non-imaging problem with the same target light, but with a different source (since it receives only a portion of the light). Hence, the optical component made with all the pillows glued together is more reliable and allows for example to reduce the artifacts due to small occluders. Indeed if one object is in front of one or more lenses, the quality of the refracted image decreases but the image can still be recognized. Using pillows also gives some flexibility to the designer to improve the appearance and the volume occupied by the component. An example with 9 pillows can be found in Figure 15, and the effect of a small occluder. In practice, in order to avoid a shift between the nine simulated images, we solve the near-field problem (see Section 6).

Application to color images Using pillows, we can also target color images. Indeed, we can build one component for each of the Red, Green and Blue channels of an image. If we then place three lights (red, green and blue) in front of each component, using Algorithm 3 the 3 images will be perfectly align and thus produce the original color image, see the first image of Figure 1.

Physical prototypes

We built three lenses (see Figure 16) and two mirrors (see figures 17 and 1) corresponding to collimated light sources. The lenses are fabricated in PMMA (whose index of refraction is 1.49) and the mirrors in aluminum. All the lenses and mirrors have size 100mm × 100mm and were milled in one pass on 3-axis CNC machines after milling the blank. For the Hikari lens, we choose to focalize the target image on a wall at 2 meters and the target is a square of size 600mm × 600mm. For the Train and Einstein signature lenses and the two mirrors, we choose to focalize the target image on a wall at 1 meter and the target is a square of size 300mm × 300mm. The five components were milled in one pass on the DMG-DC100V machine with a 10mm radius end-mill.

The milling process is very sensitive. For the lenses, the endmill is a super finishing ball mill D10, 3 teeth and is following a concentric spiral trajectory. For the mirrors the end-mill is a PCD ball mill hooped D10, 2 teeth and is following a parallel scanning trajectory. We observe that the precision of the milling is not accurate enough: when the collimated light source is traversing a lens or reflected by a mirror with no sandpapering and polishing, the light is dispersed and we do not recognize the target (see Figure 19, second row). We had to sandpaper them by hand before polishing them with a polishing paste. This clearly damages the lens surface: there is a tradeoff between removing the artifacts due to the milling and smoothing too much the surface (see Figure 19 rows 3-5), thus damaging its refractive properties. Note that thanks to the convexity property (see Figure 19 row 4), the lens surface is quite regular and is more robust to sandpapering.

We can also observe some artifacts in the milling process. For instance, some corrugations are present in the lens (Figure 19 first row) and induce some artifacts in the projected image (Figure 16, first row). We observe that although the image are very contrasted, the projected image are very accurate. The boundary of the target is often slightly blurred and this is due to the boundary of the lens or mirror where the milling was less good. Our model do not take into account the different wavelengths of the white color and we observe on the boundary of the projected images a small chromatic aberration (the boundary is slightly blue).

Conclusion and perspectives

We presented a general framework for eight different optical component design problems satisfying light energy constraints. We proposed an efficient algorithm able to solve them whether the target is at infinity (far-field) or at a finite distance (nearfield). The main limitation of the approach is the fact that we only deal with ideal light sources (a light bulb is for instance neither collimated nor punctual). Another limitation is that we do not account for self shadowing and internal reflections (although, this is not a problem in the situations we have encountered). In the future, we also want to try fabricating physical prototypes when the source is punctual. This problem becomes more difficult for the lenses since we also have to mill the two sides of the lens, and in particular the inner sphere. We also believe that the robustness and versatility of the proposed approach can make it a useful component for the design of heuristics able to deal with extended light sources and in computer graphics for caustic design.

Figure 2 :

 2 Figure2: Four inverse problems arising in non-imaging optics. In each case, the goal is to build the surface R of a mirror or a lens. For each problem, we provide two solutions (for instance, convex or concave surfaces when the light source is collimated). Top/Bottom: Collimated light sources/Point light sources. Left/Right: Mirror/Lens design.

Figure 3 :

 3 Figure 3: Convex Mirror for a collimated light source (when N = 16). The mirror surface R is the graph of a convex piecewise affine function. The support X ρ of ρ is decomposed into visibility cells (V i (ψ)) 1 i N . Every vertical ray above a point x ∈ X ρ belongs to a cell V i (ψ), touches a plane with slope p i and is reflected to the direction y i .

Figure 4 :

 4 Figure 4: Concave and convex lenses for a uniform collimated light source and the same target.

  3.2. Point light source3.2.1 Concave mirror design.

ALGORITHM 2 :

 2 Damped Newton method for G(ψ) = σ

Figure 5 :

 5 Figure 5: Mean error ( c k+1 -c k 1 /N ) in Algorithm 3 for a uniform collimated light source. Red: (CS/Lens) with the Train target; Blue: (CS/Mirror) with the Train target. The y-axis is in logarithmic scale.

Figure 6 :

 6 Figure 6: Forward simulation of different iterations of Algorithm 3. From left to right: target image, 1 st iteration; 2 nd iteration; after 6 iterations. First row: Concave (CS/Lens) Train. Second row: Convex (CS/Mirror) Train.

Figure 7 :

 7 Figure 7: Four non-imaging problems solved with Algorithm 1. From left to right: visibility diagram on X ρ (wireframe) with the optical component R, Triangulation R T of R; forward simulation using LuxRender. From top to bottom: Convex Collimated Source Mirror; Concave Point Source Mirror; Concave Collimated Source Lens; Point Source Lens (union of ellipsoids).

Figure 8 :

 8 Figure 8: Convex Collimated Source Mirror problem with a uniform light source for different target distributions. From left to right: target distribution, mean curvature of the mirror, forward simulation using LuxRender. Dimensions of images (top to bottom): 256 2 , 300 2 , 400 2 .

Figure 9 :

 9 Figure 9: Concave Point Source Mirror problem for a uniform point light source with different target distributions. From left to right: target distribution, mean curvature of the mirror (top view), forward simulation using LuxRender. Dimensions of images (top to bottom): 256 2 , 300 2 , 400 2 .

Figure 10 :

 10 Figure 10: Triangulation R T for a non-uniform light source. From left to right: non uniform collimated light source; mesh of the lens for this non-uniform light; mesh of the lens for a uniform light source.

Figure 11 :

 11 Figure 11: Concave Collimated Source Lens with a uniform light source for different target distributions. From left to right: target distribution, mean curvature of the lens (top view), forward simulation using LuxRender. Dimensions of images (top to bottom): 256 2 , 300 2 , 400 2 .

Figure 12 :

 12 Figure 12: Point Source Lens with a uniform light source for different target distributions. The lens surface is the boundary of the union of filled ellipsoids, hence is not convex, nor concave. From left to right: target distribution, mean curvature of the lens (top view), forward simulation using LuxRender. Dimensions of images (top to bottom): 256 2 , 300 2 , 400 2 .

Figure 13 :

 13 Figure13: Comparison with[START_REF] Schwartzburg | High-contrast computational caustic design[END_REF] From left to right: target distribution; images obtained by[START_REF] Schwartzburg | High-contrast computational caustic design[END_REF] and taken from their article; our forward simulation using LuxRender. Last row: meshes of the two corresponding convex lenses: Rings (left) and Siggraph (right).

Figure 14 :

 14 Figure 14: Stability under rotation of the lens. LuxRender renderings in the (CS/Lens) setting for the Train target while rotating the lens with respect to the direction of the collimated light source (0 • / 5 • / 15 • ).

Figure 15 :

 15 Figure 15: Pillows and differences between FF and NF. Top: The lens is composed of three pillows that solve the FF problem. Note in the last image the shift between the three projected images. Middle: A lens composed of nine pillows (each of them solving the NF problem) that refracts a uniform collimated light source; Bottom: The same lens with an obstacle in red.

Figure 16 :

 16 Figure 16: Fabricated lenses for a collimated light source. From left to right: experimental setup, zoom on the target screen. From top to bottom: Train, Hikari, Einstein's signature targets. Images are focused on a screen at 2 meters for the first two rows and 1 meter for the last one.

Figure 17 :

 17 Figure 17: Fabricated mirror for a collimated light source. Hikari target.

Figure 18 :

 18 Figure 18: Stability with respect to the depth of the focus plane. The lens of Hikari is designed to focus at a distance of 2 meters. The target screen is at different depths, top: 1.5 meters; bottom: 1m (left), 50cm (right).

Table 1 :

 1 Formulas for the weighted points used to define the Power cells in Equation (4) for the various problems. In the lens design problem, κ > 0 is the ratio of the indices of refraction, κ > 1 in the (PS/Lens) setting. Ccv means concave and Cvx convex. Ccv means that the optical component converges to a concave one when the discretization tends to infinity.

	Type	Points	Weights
	Cvx (CS/Mirror)

Table 1 and the work of de Castro et al., 2016 for formulas in the (PS/Mirror) case, the other ones being obtained in a similar fashion. Let us show the derivation of

Table 2 :

 2 Running time and number of Newton steps in Algorithm 2 (FF target) for the (CS/Mirror) and the Train target in the far-field setting.

	size	time	# Newton steps
	128 2	9s	11 iterations
	256 2	38s	13 iterations
	512 2	245s	15 iterations
	1024 2 1598s	18 iterations
	2048 2 7538s	24 iterations

Table 3 :

 3 Running time in Algorithm 3 (NF target) for the (CS/Mirror) and the Train target in the near-field setting for different discretizations.

Acknowledgements We would like to thank Alain Di Donato from the GINOVA platform, S.MART DS department for building the physical prototypes, as well as the GMP department of the Institute of Technology of Aix-en-Provence in France, and André Lieutier and Jimmy Cresson from Dassault-Systèmes for helping with the fabrication process, and Philippe Halot for taking the pictures. This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) and by ANR project MAGA (ANR-16-CE40-0014 MAGA).


P L Krapivsky
email: pkrapivsky@gmail.com

J M Luck

K Mallick
email: kirone.mallick@ipht.fr

Quantum return probability of a system of N non-interacting lattice fermions

We consider N non-interacting fermions performing continuous-time quantum walks on a one-dimensional lattice. The system is launched from a most compact configuration where the fermions occupy neighboring sites. We calculate exactly the quantum return probability (sometimes referred to as the Loschmidt echo) of observing the very same compact state at a later time t. Remarkably, this probability depends on the parity of the fermion number -it decays as a power of time for even N , while for odd N it exhibits periodic oscillations modulated by a decaying power law. The exponent also slightly depends on the parity of N , and is roughly twice smaller than what it would be in the continuum limit. We also consider the same problem, and obtain similar results, in the presence of an impenetrable wall at the origin constraining the particles to remain on the positive half-line. We derive closed-form expressions for the amplitudes of the power-law decay of the return probability in all cases. The key point in the derivation is the use of Mehta integrals, which are limiting cases of the Selberg integral.

Introduction

The advent of techniques to manipulate cold atoms has led to the experimental realization of low-dimensional quantum gases [START_REF] Bloch | Quantum simulations with ultracold quantum gases[END_REF] that, for many decades, were mainly thought of as being toy models for theoreticians [START_REF] Paredes | Tonks-Girardeau gas of ultracold atoms in an optical lattice[END_REF][START_REF] Kinoshita | A quantum Newton's cradle[END_REF]. The continuous tuning of the pair interaction between atoms thanks to Feshbach resonances allows one to create many-body systems ranging from free to strongly interacting. This revolution triggered a profusion of studies of the non-equilibrium dynamics of quantum systems [START_REF] Cazalilla | Focus on dynamics and thermalization in isolated quantum many-body systems[END_REF][START_REF] Polkovnikov | Nonequilibrium dynamics of closed interacting quantum systems[END_REF][START_REF] Eisert | Quantum many-body systems out of equilibrium[END_REF], including quenches, entanglement and thermalization, that resulted in many new ways of exploring the quantum world [START_REF] Haroche | Exploring the Quantum: Atoms, Cavities, and Photons[END_REF].

Non-interacting tight-binding fermions on a lattice form the simplest of all quantum gases. They can be viewed as continuous-time quantum walkers. Quantum walks were initially introduced in the context of quantum information theory, as a general framework to state quantum algorithms [START_REF] Aharonov | Quantum random walks[END_REF][START_REF] Ambainis | Quantum walks and their algorithmic applications[END_REF][START_REF] Kempe | Quantum random walks: An introductory overview[END_REF]. Many fundamental problems that have been studied for classical random walks have natural counterparts in the realm of quantum walks, often leading to results that look unexpected and surprising from a classical point of view. First of all, quantum walkers behave ballistically rather than diffusively [START_REF] Farhi | Quantum computation and decision trees[END_REF][START_REF] De Toro Arias | Anomalous dynamical scaling and bifractality in the onedimensional Anderson model[END_REF]. The survival in the presence of traps [START_REF] Parris | One-dimensional trapping kinetics at zero temperature[END_REF][START_REF] Krapivsky | Survival of classical and quantum particles in the presence of traps[END_REF] and the ballistic spreading of bound states [START_REF] Krapivsky | Interacting quantum walkers: Two-body bosonic and fermionic bound states[END_REF][START_REF] Krapivsky | Quantum centipedes: Collective dynamics of interacting quantum walkers[END_REF] provide yet other examples of qualitatively different behavior in the classical and quantum cases.

A foremost question that has been investigated since the earlier days of quantum mechanics concerns the decay of individual quantum states [START_REF] Gamow | Zur Quantentheorie des Atomkernes[END_REF][START_REF] Gurney | Quantum mechanics and radioactive disintegration[END_REF][START_REF] Weisskopf | Berechnung der natürlichen Linienbreite auf Grund der Diracschen Lichttheorie[END_REF][START_REF] Weisskopf | Über die natürliche Linienbreite in der Strahlung des harmonischen Oszillators[END_REF][START_REF] Khalfin | Contribution to the decay theory of a quasi-stationary state[END_REF][START_REF] Winter | Evolution of a quasi-stationary state[END_REF][START_REF] Adler | Weisskopf-Wigner decay theory for the energy-driven stochastic Schrödinger equation[END_REF]. In many circumstances the return probability of a quantum system to its very initial state, also referred to as the Loschmidt echo [START_REF] Peres | Stability of quantum motion in chaotic and regular systems[END_REF][START_REF] Jalabert | Environment-independent decoherence rate in classically chaotic systems[END_REF][START_REF] Manfredi | Loschmidt echo in a system of interacting electrons[END_REF][START_REF] Gorin | Dynamics of Loschmidt echoes and fidelity decay[END_REF], falls off exponentially in time. Quantum revival therefore corresponds to a rare event, somehow analogously to those described by the laws of large deviations in classical statistical mechanics [START_REF] Derrida | Non-equilibrium steady states: Fluctuations and large deviations of the density and of the current[END_REF][START_REF] Touchette | The large deviation approach to statistical mechanics[END_REF]. The return probability has been recently studied for various systems of bosons or fermions confined to one dimension. There, a wide variety of temporal decays can be found, ranging from a power law to a superexponential decay. These recent works include two-particle systems [START_REF] Taniguchi | Escape behavior of quantum two-particle systems with Coulomb interactions[END_REF][START_REF] García-Calderón | Time evolution of decay of two identical quantum particles[END_REF][START_REF] Hunn | Tunneling decay of two interacting bosons in an asymmetric double-well potential: A spectral approach[END_REF], many-particle systems [START_REF] Del Campo | Long-time behavior of many-particle quantum decay[END_REF][START_REF] Longhi | Many-particle quantum decay and trapping: The role of statistics and Fano resonances[END_REF][START_REF] Marchewka | Role of quantum statistics in multi-particle decay dynamics[END_REF][START_REF] Del Campo | Exact quantum decay of an interacting many-particle system: The Calogero-Sutherland model[END_REF], and systems with infinitely many particles [START_REF] Stéphan | Local quantum quenches in critical one-dimensional systems: Entanglement, the Loschmidt echo, and light-cone effects[END_REF][START_REF] Viti | Inhomogeneous quenches in a free fermionic chain: Exact results[END_REF][START_REF] Stéphan | Return probability after a quench from a domain-wall initial state in the spin-1/2 XXZ chain[END_REF]. For instance, for free lattice fermions with domain-wall initial condition, the return probability was found to obey a pure Gaussian decay in time [START_REF] Stéphan | Local quantum quenches in critical one-dimensional systems: Entanglement, the Loschmidt echo, and light-cone effects[END_REF][START_REF] Viti | Inhomogeneous quenches in a free fermionic chain: Exact results[END_REF]. More recently, the XXZ quantum spin chain with domain-wall initial condition has been investigated in the whole massless phase, where the anisotropy parameter obeys |∆| < 1 [START_REF] Stéphan | Return probability after a quench from a domain-wall initial state in the spin-1/2 XXZ chain[END_REF]. There, the decay of the return probability is found to be either Gaussian or exponential, depending on whether the angle parameter γ (defined by ∆ = cos γ) is commensurate to π or not. This highly discontinuous asymptotic behavior is attributed to integrability.

The aim of the present work is to study the return probability for a system of N non-interacting fermions hopping on a one-dimensional lattice. The fermions are launched from the most compact state where they occupy N successive lattice sites. We consider two different settings: an infinite chain (section 2) and a semiinfinite chain bounded by an impenetrable wall (section 3). In the latter situation, the compact initial state lies near the wall. We first derive general expressions for the return probability at arbitrary finite times (sections 2.1 and 3.1). We then perform an asymptotic analysis of the regime of late times (sections 2.3 and 3.2). In both settings the return probability manifests a dependence on the parity of the fermion number, oscillating forever for odd N and decaying monotonically for even N . In all cases it falls off as a power of time, whose exponent also slightly depends on the parity of N , and is roughly twice smaller than what it would be in the continuum limit. We also derive closed-form expressions for the amplitudes of the power-law decay of the return probability in all these situations. This derivation relies on the usage of Mehta integrals [START_REF] Mehta | Random Matrices[END_REF]Ch. 17]. Our results suggest a non-trivial crossover behavior in the scaling regime where time t and the fermion number N are both large and proportional to each other (sections 2.4 and 3.3). The determination of the corresponding scaling functions F and F (w) remains a challenging open problem. Section 4 contains a brief summary and a discussion of our findings. Two appendices are respectively devoted to the Andréief identity (Appendix A) and to Barnes' G-function (Appendix B).

Free fermions on the infinite chain

We consider N non-interacting tight-binding fermions hopping on an infinite chain. At the initial time (t = 0), the fermions are launched from N consecutive sites, labeled n = 1, . . . , N (see figure 1).

N

Figure 1. The initial configuration of the N fermions on the chain.

We are interested in the return probability R N (t), i.e., the probability that the fermions occupy the same sites at a later time t, and especially in the asymptotic decay of this quantity in the regime of long times.

General expressions at finite times

General expressions for the return probability R N (t) at arbitrary finite times can be derived as follows. Within the formalism of the second quantization, the Hamiltonian of the system reads

H = n a † n a n+1 + a † n+1 a n , (2.1)
in dimensionless units, with the standard fermionic anti-commutation rules

{a m , a n } = {a † m , a † n } = 0, {a m , a † n } = δ m,n . (2.2)
The many-body wavefunction |Ψ(t) obeys the time-dependent Schrödinger equation

i d|Ψ dt = H|Ψ , (2.3)
and so

|Ψ(t) = e -iHt |Ψ(0) . (2.4)
The return probability therefore reads

R N (t) = |A N (t)| 2 , (2.5)
with

A N (t) = Ψ(0)|Ψ(t) = Ψ(0)|e -iHt |Ψ(0) . (2.6)
The above expression can be made explicit by bringing the Hamiltonian H to the diagonal form

H = 2π 0 dq 2π ε q a † q a q , (2.7)
with ε q = 2 cos q (0 ≤ q ≤ 2π).

(2.8)

The operators

a q = n e inq a n , a † q = n e -inq a † n (2.9)
obey the anti-commutation relation { a q , a † q ′ } = 2π δ(q -q ′ ).

(2.10)

We find after some elementary algebra

A N (t) = 2π 0 • • • 2π 0 N n=1 dq n 2π e -2it cos qn | Ψ(0)|q | 2 , (2.11)
where q = (q 1 , . . . , q N ), and the q n are the momenta of the N fermions. The manybody amplitude Ψ(0)|q is given by a normalized Slater determinant. In the present situation, where the fermions are launched from the sites m = 1, . . . , N , we have ‡

Ψ(0)|q = 1 √ N ! det e imqn , (2.12)
and therefore

A N (t) = 1 N ! 2π 0 • • • 2π 0 N n=1
dq n 2π e -2it cos qn det e imqn 2 .

(2.13)

The formula (2.12) can be recast as

Ψ(0)|q = 1 √ N ! N n=1 e iqn 1≤m<n≤N e iqn -e iqm , (2.14)
where the last product is a Vandermonde determinant. We are thus left with the following integral expression for the amplitude:

A N (t) = 1 N ! 2π 0 • • • 2π 0 N n=1 dq n 2π e -2it cos qn × 1≤m<n≤N e iqn -e iqm 2 .
(2.15)

Applying the Andréief identity (A.1) to (2.13), we get the alternative expression

A N (t) = det (ψ m,n (t)) = det (J m-n (2t)) , (2.16)
where is the one-body wavefunction at site n at time t for a particle launched from site m at time 0, and the J n are Bessel functions. The expressions (2.15) and (2.16) are complementary. The first one is more amenable to some analytical investigations, including the asymptotic analysis of the long-time regime performed in section 2.3. The second one could have been obtained by more elementary means, along the lines of [START_REF] Krapivsky | Survival of classical and quantum particles in the presence of traps[END_REF][START_REF] Krapivsky | Interacting quantum walkers: Two-body bosonic and fermionic bound states[END_REF][START_REF] Krapivsky | Quantum centipedes: Collective dynamics of interacting quantum walkers[END_REF]. It demonstrates explicitly that the amplitude is properly normalized, i.e., it obeys A N (0) = 1. It is also more suited for numerical evaluations and power-series expansions in t.

ψ m,n (t) = 2π 0 dq 2π e -2it cos q+i(n-m)q = i m-n J n-m (2t) (2.

The first few values of N

It is interesting to first look at the first few values of the fermion number N .

• N = 1. For a single particle, we have A 1 (t) = J 0 (2t), and so

R 1 (t) = J 2 0 (2t) ≈ cos 2 (2t -π/4) πt .
(2.18)

The return probability therefore oscillates forever, becoming exactly zero at infinitely many times t k ≈ (k -1/4)π/2. These oscillations can be averaged out by replacing the squared cosine by a factor 1/2. We thus obtain the power-law decay

R 1 (t) ≈ 1 2πt (2.19)
for the mean return probability.

• N = 2. The situation of two fermions is the simplest one where quantum statistics plays a role. We have A 2 (t) = J 2 0 (2t) + J 2 1 (2t), and so

R 2 (t) = (J 2 0 (2t) + J 2 1 (2t)) 2 ≈ 1 π 2 t 2 .
(2.20)

In this case the return probability falls off monotonically to zero.

• N = 3. We have

A 3 (t) = (J 0 (2t) + J 2 (2t))(J 2 0 (2t) + 2J 2 1 (2t) -J 0 (2t)J 2 (2t)) = 2J 1 (2t) t (J 2 0 (2t) + J 2 1 (2t)) - J 0 (2t)J 2 1 (2t) t 2 , (2.21)
where the first line gives the raw determinantal expression (2.16). The second one, obtained by means of the recursion

J n+1 (2t) + J n-1 (2t) = n t J n (2t), (2.22)
is more suitable to study the late-time behavior of the return probability. Keeping only the first group of terms in the second line of (2.21), we obtain

R 3 (t) ≈ 4 cos 2 (2t -3π/4) π 3 t 5 . ( 2

.23)

The return probability again oscillates forever, becoming exactly zero at infinitely many times t k ≈ (k + 1/4)π/2. The situation is therefore qualitatively similar to that of a single particle. Averaging again over the oscillations, we obtain

R 3 (t) ≈ 2 π 3 t 5 .
(2.24)

Asymptotic analysis in the long-time regime

In this section we determine the asymptotic decay of the return probability R N (t) in the long-time regime, for arbitrary values of the fermion number N . To do so, we evaluate the multiple integral entering (2.15) by means of the saddle-point approximation. Saddle points are defined by the condition that every momentum q n is either 0 or π. The vicinity of the most general saddle point can thus be parametrized by choosing M momenta near 0, of the form q n = x n , and the remaining N -M momenta around π, of the form q n = π + y n . For a fixed M , there are N M ways of choosing which momenta are near 0 and near π. Expanding in (2.15) the arguments of the exponentials and the products to quadratic order in the variables x n and y n , we obtain the estimate

A N (t) ≈ N M=0 I N,M (t) M !(N -M )! , (2.25)
with

I N,M (t) = 2 2M(N -M) e 2i(N -2M)t × ∞ -∞ • • • ∞ -∞ M n=1 dx n 2π e itx 2 n 1≤m<n≤M (x n -x m ) 2 (2.26) × ∞ -∞ • • • ∞ -∞ N -M n=1 dy n 2π e -ity 2 n 1≤m<n≤N -M (y n -y m ) 2 .
The above integrals can be performed exactly. They are indeed given be the analytical continuation to imaginary values of a of the Mehta integral [40, Eq. (17.6.7), γ = 1]

I N (a) = ∞ -∞ • • • ∞ -∞ N n=1 dx n e -ax 2 n 1≤m<n≤N (x n -x m ) 2 = (2π) N/2 (2a) -N 2 /2 N j=1 j! = (2π) N/2 (2a) -N 2 /2 G(N + 2), (2.27)
where G is Barnes' G-function (see Appendix B). We thus obtain

I N,M (t) = 2 2M(N -M) (2π) -N/2 e i(N -2M)(2t-N π/4) × G(M + 2)G(N -M + 2)(2t) M(N -M)-N 2 /2 .
(2.28)

All the quantities I N,M (t) which enter the estimate (2.25) thus fall off as power laws in time, albeit with an M -dependent exponent. The behavior of the return probability at long times is governed by the term I N,M (t) with the slowest decay. Even and odd values of the fermion number N yield different behaviors, and have to be dealt with separately. The emerging picture fully corroborates the observations made in section 2.2 for the first few values of N .

• If N = 2m is even, the slowest decay is reached for M = m, i.e., for equal numbers of momenta near 0 and π. We thus obtain the estimate

A 2m (t) ≈ 2 m(m-1) G(m + 1) 2 π m t m 2 , (2.29)
showing that the amplitude is positive and decays monotonically, at least for large times. The echo therefore falls off monotonically as a power law

R 2m (t) ≈ 2 2m(m-1) G(m + 1) 4 π 2m t 2m 2 , (2.30)
with exponent 2m 2 = 1 2 N 2 . • If N = 2m + 1 is odd, the slowest decay is reached both for M = m and M = m + 1. We thus obtain the estimate

A 2m+1 (t) ≈ 2 m 2 G(m + 1)G(m + 2) π m+1/2 t m 2 +m+1/2 cos(2t -(2m + 1)π/4), (2.31)
showing that the amplitude behaves for large times as an oscillatory function modulated by a decaying power law. Averaging over the oscillations, we obtain the following power-law decay for the mean echo:

R 2m+1 (t) ≈ 2 2m 2 -1 G(m + 1) 2 G(m + 2) 2 π 2m+1 t 2m 2 +2m+1 , (2.32)
with exponent 2m 2 + 2m + 1 = 1 2 (N 2 + 1). Figure 2 illustrates the above results, showing log-log plots of the return probability R N (t) against time t for N = 1, 2, 3 and 4. Data have been obtained using the expression (2.16). For odd N (left), the echo oscillates forever. For even N (right), it falls off monotonically but exhibits damped oscillations, which can be attributed to the first subleading contributions I N,M (t) for N = 2m and M = m ± 1. The blue straight lines -slightly translated for a better readability -have the predicted slopes 1, 2, 5 and 8. It is worth comparing the above results for lattice fermions to the corresponding predictions in the continuum limit. Within the present framework, taking the continuum limit just amounts to approximating the dispersion relation (2.8) by a quadratic law of the form

ε q ≈ 2 -q 2 , (2.33)
and therefore to restricting the saddle-point approximation to the sector where all momenta are near zero (i.e., M = N ). This yields

R (cont) N (t) ≈ G(N + 1) 2 2 N (N +1) π N t N 2 .
(2.34)

We still obtain a power-law decay of the return probability, albeit with no parity effect and with a different -roughly twice larger -decay exponent.

The above decay of the return probability in the continuum can be alternatively derived by the following heuristic reasoning. Suppose the fermions have mass m and are launched from the interval 0 < x < ℓ. In reduced units (h = 1), the typical momentum of each particle is p ∼ 1/ℓ, and so its wavefunction spreads ballistically over a region of size

L(t) ∼ t mℓ . (2.35)
The return probability of one single particle therefore scales as R 1 (t) ∼ ℓ/L(t) ∼ mℓ 2 /t. The exponent of this decay is in agreement with (2.32) and (2.34). In the same setting, the modulus of the wavefunction of N non-interacting fermions at late times can be estimated as

|ψ N (x 1 , . . . , x N ; t)| ∼ C N (t) 1≤m<n≤N |x m -x n | , (2.36)
as long as all the coordinates x n are less than L(t), whereas it falls off very fast at larger distances. Dimensional analysis implies that the normalization of the wavefunction scales as C N (t) ∼ L(t) -N 2 /2 , and that the return probability scales as

R N (t) ∼ ℓ L(t) N 2 ∼ mℓ 2 t N 2
.

(2.37)

The exponent of this decay is in agreement with (2.34). Furthermore, comparing the prefactors of the estimates (2.34) and (2.37) at large N yields the identification mℓ 2 ≈ N/(2e 3/2 ). The scaling ℓ ∼ √ N of the length of the confining interval is an artifact of the continuum framework.

Large-N asymptotics

In the N → ∞ limit, the amplitude A N (t) admits the following remarkably simple expression for all finite times:

A ∞ (t) = lim N →∞ det (J m-n (2t)) = e -t 2 .
(2.38)

This expression has been derived independently in two recent works, devoted to quantum quenches of fermionic chains [START_REF] Viti | Inhomogeneous quenches in a free fermionic chain: Exact results[END_REF], and to volumes of balls in unitary groups [START_REF] Wei | From random matrix theory to coding theory: Volume of a metric ball in unitary groups[END_REF]. The return probability therefore reads

R ∞ (t) = e -2t 2 .
(2.39)

It can be verified, by expanding the result (2.16) as a power series in t for the first values of N , that the expressions for R N (t) and R ∞ (t) start differing at order t 2N +2 . This phenomenon has already been noticed in [START_REF] Wei | From random matrix theory to coding theory: Volume of a metric ball in unitary groups[END_REF].

On the other hand, for large but finite N , the results of section 2.3 can be given more explicit forms. Using the asymptotic expansion (B.5), the expressions (2.30) for N even and (2.32) for N odd respectively yield

N even: ln R N (t) = - N 2 2 ln t N + 3 2 - 1 3 ln N + 1 3 ln 2 + 4ζ ′ (-1) + • • • , (2.40)
N odd: ln R N (t) = - N 2 + 1 2 ln t N + 3 2 - 1 3 ln N - 2 3 ln 2 + 3 4 + 4ζ ′ (-1) + • • • , (2.41)
where the remainders go to zero for large N . To leading order, both expressions read

ln R N (t) ≈ - N 2 2 ln t N + 3 2 . (2.42)
The results (2.39) and (2.42) suggest a scaling law of the form

ln R N (t) ≈ - N 2 2 F (x), x = t N , (2.43)
all over the regime where t and N are large and comparable, with

F (x) ≈ 4x 2 (x ≪ 1), ln x + 3 2 (x ≫ 1). (2.44)
The form of the scaling variable x reflects the ballistic nature of the dynamics of a free tight-binding particle.

Free fermions near a wall

In this section we consider the same problem on a semi-infinite chain ending with an impenetrable wall. At time t = 0, the fermions are launched from the N sites which are closest to the wall (see figure 3).

1 N We are interested in the temporal decay of the return probability R (w)

N (t)
, where the superscript (w) reminds of the presence of a wall.

General expressions at finite times

General expressions for the return probability R The presence of an impenetrable wall imposes Dirichlet boundary conditions at site n = 0. With these boundary conditions, the Hamiltonian H is brought to the diagonal form

H = π 0 dq π ε q a † q a q , (3.1)
with

ε q = 2 cos q (0 ≤ q ≤ π). (3.2)
The operators

a q = √ 2 n sin nq a n , a † q = √ 2 n sin nq a † n (3.3)
obey the anti-commutation relation

{ a q , a † q ′ } = π δ(q -q ′ ). (3.4)
The many-body amplitude Ψ(0)|q again given by a normalized Slater determinant. Since the fermions are launched from the sites m = 1, . . . , N , the latter reads

Ψ(0)|q = 2 N N ! det (sin mq n ) . (3.5)
The analogue of (2.11) therefore reads

A (w) N (t) = 2 N N ! π 0 • • • π 0 N n=1 dq n π e -2it cos qn (det (sin mq n )) 2 . (3.6)
The formula (3.5) can be brought to a form similar to (2.14). We recall that sin mq = sin q U m-1 (cos q), (

where the U m are the Tchebyshev polynomials of the second kind [42, Vol. II, Ch. X].

The mth polynomial has degree m and the leading term

U m (z) = (2z) m + • • • We have therefore det (sin mq n ) = N n=1 sin q n det (U m-1 (q n )) = 2 N (N -1)/2 N n=1 sin q n × 1≤m<n≤N (cos q n -cos q m ), (3.8)
where the last product is again a Vandermonde determinant. We are thus left with the following integral expression for the amplitude:

A (w) N (t) = 2 N 2 N ! π 0 • • • π 0 N n=1 dq n π sin 2 q n e -2it cos qn × 1≤m<n≤N (cos q n -cos q m ) 2 . (3.9)
Applying the Andréief identity (A.1) to (3.6), we get the alternative expression

A (w) N (t) = det ψ (w) m,n (t) , (3.10)
where

ψ (w) m,n (t) = ψ m,n (t) -ψ -m,n (t) = i -(n-m) J n-m (2t) -i -(n+m) J n+m (2t) (3.11)
is the one-body wavefunction at site n at time t for a particle located on the semiinfinite chain and launched from site at time 0. The expression (3.11) can be recovered by applying the reflection principle (i.e., the method of images) to (2.17). The complementary expressions (3.9) and (3.10) are the exact analogues on the semi-infinite chain of their counterparts (2.15) and (2.16) on the infinite chain.

Asymptotic analysis in the long-time regime

In this section we determine the asymptotic decay of the return probability R (w) N (t). We again apply the saddle-point approximation to the multiple integral in (3.9). Saddle points are still defined by the condition that every momentum q n is either 0 or π. The vicinity of the most general saddle point will be parametrized by choosing M momenta near 0, of the form q n = x n , and the remaining N -M momenta around π, of the form q n = π -y n . At variance with the previous situation (section 2.3), the variables x n and y n are now positive. We thus obtain the estimate

A (w) N (t) ≈ 2 N 2 N M=0 I (w) N,M (t) M !(N -M )! , (3.12)
with

I (w) N,M (t) = 2 N -(N -2M) 2 e 2i(N -2M)t × ∞ 0 • • • ∞ 0 M n=1 dx n π x 2 n e itx 2 n 1≤m<n≤M (x 2 n -x 2 m ) 2 (3.13) × ∞ 0 • • • ∞ 0 N -M n=1 dy n π y 2 n e -ity 2 n 1≤m<n≤N -M (y 2 n -y 2 m ) 2 .
The above integrals can still be performed exactly. They are indeed proportional to the analytical continuation to imaginary values of a of the Mehta integral [40, Eq. (17.6.6), α = 3/2, γ = 1]

I (w) N (a) = ∞ -∞ • • • ∞ -∞ N n=1 dx n x 2 n e -ax 2 n 1≤m<n≤N (x 2 n -x 2 m ) 2 = a -N (2N +1)/2 N j=1 j! Γ(j + 1/2) = π N/2 (2a) -N (2N +1)/2 N ! G(2N + 2), (3.14)
where G is again Barnes' G-function (see Appendix B). We thus obtain

I (w) N,M (t) = 2 -(N -2M) 2 π -N/2 e i(N -2M)(2t-(N +1/2)π/4) × M !(N -M )! G(2M + 2)G(2N -2M + 2) × (2t) 2M(N -M)-N (2N +1)/2 . (3.15)
The behavior of the return probability at long times is still governed by the term with the slowest decay. Even and odd values of the fermion number N again yield different behaviors.

• If N = 2m is even, the slowest decay is reached for M = m, i.e., for equal numbers of momenta near 0 and π. We thus obtain the estimate

A (w) 2m (t) ≈ 2 m(2m-1) G(2m + 2) π m m! t m(2m+1) , (3.16)
showing that the amplitude is positive and decays monotonically, at least for large times. The echo therefore falls off monotonically as a power law

R (w) 2m (t) ≈ 2 2m(2m-1) G(2m + 2) 2 π 2m m! 2 t 2m(2m+1) , (3.17)
with exponent 2m(2m + 1) = N (N + 1).

• If N = 2m + 1 is odd, the slowest decay is reached both for M = m and M = m + 1.

We thus obtain the estimate (3.18) showing that the amplitude behaves for large times as an oscillatory function modulated by a decaying power law. Averaging over the oscillations, we obtain the following power-law decay for the mean echo:

A (w) 2m+1 (t) ≈ 2 m(2m+1) G(2m + 3) π m+1/2 m! t 2m 2 +3m+3/2 cos(2t -(2m + 3/2)π/4),
R (w) 2m+1 (t) ≈ 2 2m(2m+1)-1 G(2m + 3) 2 π 2m+1 m! 2 t 4m 2 +6m+3 , (3.19)
with exponent 4m 2 + 6m + 3 = N 2 + N + 1. Figure 4 illustrates the above results, showing log-log plots of the return probability R (w) N (t) against time t for N = 1, 2, 3 and 4. Data have been obtained using the expression (3.10). For odd N (left), the echo oscillates forever. For even N (right), it falls off monotonically but exhibits damped oscillations, which can again be attributed to the first subleading contributions I (w) N,M (t) for N = 2m and M = m ± 1. The blue straight lines -slightly translated for a better readability -have the predicted slopes 3, 6, 13 and 20.

It is again worth comparing the above results for lattice fermions to the corresponding prediction in the continuum limit. Restricting the saddle-point approximation to the sector where all momenta are near zero (i.e., M = N ) yields

R (w,cont) N (t) ≈ G(2N + 2) 2 N (2N +1) π N N ! t N (2N +1) . (3.20)
We again obtain a power-law decay of the return probability, albeit with no parity effect and with a different, larger decay exponent. The above power-law decay can still be recovered by means of heuristic reasoning. Suppose the fermions are launched from the interval 0 < x < ℓ near the wall. The modulus of the many-body wavefunction at time t now reads approximately

ψ (w) N (x 1 , . . . , x N ; t) ∼ C (w) N (t) N n=1 x n 1≤m<n≤N x 2 m -x 2 n . (3.21)
Dimensional analysis determines the scaling of the normalization of the wavefunction, C +1) , and of the return probability, R A very similar result can be found in [START_REF] Del Campo | Long-time behavior of many-particle quantum decay[END_REF]. The decay exponent of the above estimate agrees with (3.20). Furthermore, comparing the prefactors of the estimates (3.20) and (3.22) at large N yields mℓ 2 ≈ N/e 3/2 . The latter estimate is twice larger than its counterpart on the infinite line, given below (2.37). An intuitive interpretation of this factor two will be given below (3.27).

N (t) ∼ L(t) -N (2N (w)
N (t) ∼ ℓ L(t) N (2N +1) ∼ mℓ 2 t N (2N +1) . ( (w)

Large-N asymptotics

In the N → ∞ limit, the return probability is expected to be equal to the square root of its counterpart (2.39) on the infinite chain, namely R (w) ∞ (t) = e -t 2 .

(3.23)

A heuristic way of showing this goes as follows. On the infinite chain, the compact fermionic state has two ends, and can therefore decay through either end, whereas it has only one right end if confined near a wall. Hence we can expect R ∞

(t) = (R (w)
∞ (t)) 2 . The latter result can be checked by expanding the expression (3.10) as a power series in t for the first values of N . Doing so confirms our expectation and shows that R ∞ (t) again start differing at order t 2N +2 . On the other hand, using the asymptotic formula (B.5), we can still derive more explicit forms of the above results for large N . The expressions (3.17) for N even and (3.19) for N odd respectively yield

N even: ln R (w) N (t) = -N (N + 1) ln t 2N + 1 + 3 2 - 1 6 ln N + ln 2 - 3 8 + 2ζ ′ (-1) + • • • , (3.24)
N odd: ln R (w)

N (t) = -(N 2 + N + 1) ln t 2N + 1 + 3 2 - 1 6 ln N -2 ln 2 + 9 8 + 2ζ ′ (-1) + • • • (3.25)
The argument t/(2N + 1) of the logarithms has been chosen in order to minimize the order of magnitude of the correction terms given in the second lines of the above expressions. To leading order, both results read ln R

N (t) ≈ -N 2 F (w) (x), x = t 2N . (3.27)
The effective fermion number 2N entering the scaling variable x can be interpreted as the total distance to be traveled by an excitation entering the compact fermionic state from the right, bouncing at the wall, and exiting from the right. Finally, the scaling function F (w) (x) obeys the very same asymptotics (2.44) as F (x) both for x ≪ 1 and for x ≫ 1. This observation might be more than a coincidence, in the sense that both scaling functions might be identical.

Discussion

We studied the quantum return probability for a system of N non-interacting lattice fermions launched from N consecutive sites, either on the infinite chain or near an impenetrable wall ending a semi-infinite chain.

In each case we derived exact expressions for the return probability valid for all fermion numbers N and time t. We thus obtained two complementary kinds of expressions, namely integral formulas (see (2.15) and (3.9)), which are the natural outcome of the second-quantized formalism, and determinantal formulas (see (2.16) and (3.10)), which could have been obtained by a more elementary, first-quantized approach as well. We deduced the asymptotic long-time behavior of the return probability by evaluating the integral formulas by means of the saddle-point method.

Even and odd values of the fermion number N yield different qualitative behaviors, as well as slightly different expressions for the decay exponents. For even N , the echo falls off monotonically as a power law. For odd N , it exhibits periodic oscillations modulated by a decaying power law. This qualitative dependence on the parity of the fermion number is a pure lattice effect, which is absent in the continuum limit. The return probability of N particles thus provides yet another example of a situation where quantum dynamics exhibits qualitatively different features on the lattice and in the continuum. The exponents characterizing the temporal decay of the return probability are gathered in table 1. infinite chain infinite chain near a wall near a wall (lattice) (continuum) (lattice) (continuum)

N even: 1 2 N 2 N odd: 1 2 (N 2 + 1) N 2 N even: N (N + 1) N odd: N 2 + N + 1 N (2N + 1)
Table 1. Exponents of the temporal decay of the return probability on the infinite chain and near a wall. Comparison between the values of the exponents for lattice fermions and the predictions of the continuum limit.

Our results also yield explicit expressions for the prefactors of the asymptotic power-law decay of the return probability, in both geometries and for all values of the fermion number N . The key point in the derivation of these results has been the use of Mehta integrals, which have been extensively used in random matrix theory [START_REF] Mehta | Random Matrices[END_REF]Ch. 17] and can be viewed as limiting cases of the Selberg integral. Reference [START_REF] Forrester | The importance of the Selberg integral Bull[END_REF] provides a historical overview of the Selberg and related integrals, whereas the recent work [START_REF] Dean | Non-interacting fermions at finite temperature in a d-dimensional trap: Universal correlations[END_REF] and the references therein mention yet other connections between random matrix theory and systems of free fermions. Our expressions for the prefactors of the power-law decays involve as an essential ingredient Barnes' G-function, which is also ubiquitous in random matrix theory. The behavior of our results at large fermion numbers led us to hypothesize that the return probability obeys the scaling laws (2.43) and (3.27), with arguments x = t/N and x = t/(2N ), in the regime where t and N are both large and comparable. The form of these scaling variables reflects the ballistic nature of the motion of a free tight-binding particle. It would clearly be desirable to investigate this ballistic scaling regime in a more thorough fashion, in order to eventually determine the analytic form of the scaling functions F and F (w) , and thus prove or disprove their conjectural equality suggested by the similarity of their behavior at small and large values of their argument. Such an investigation would certainly imply the usage of more advanced techniques. There is indeed no simple way to generalize the saddle-point approach used here in order to encompass the whole ballistic scaling regime.

The framework of the present study could also be extended in several other directions. One could investigate the return probability of N fermions launched from other types of localized states, either on the chain, such as compactly supported but not most compact states, or in the continuum, with arbitrary initial states whose many-body wavefunction is localized in the center-of-mass coordinate. A further direction of research consists in replacing free fermions by other integrable systems of interacting particles in one dimension. Finally, higher-dimensional settings are also worth being envisaged. We intend to return to some of these matters in the future.

17

 17

  ) ‡ Throughout this work, determinants are of size N × N , with indices in the range 1 ≤ m, n ≤ N .

Figure 2 .

 2 Figure 2. Log-log plots of the return probability RN (t) on the infinite chain against time t, for N = 1, 2, 3 and 4. Left: odd N . Right: even N . The blue straight lines have the predicted slopes 1, 2, 5 and 8.

Figure 3 .

 3 Figure 3. The initial configuration of the N fermions near a wall ending a semi-infinite chain.

N

  (t) at arbitrary finite times can be derived along the very lines of section 2.1.

Figure 4 .

 4 Figure 4. Log-log plots of the return probability R

N

  (t) near a wall ending a semi-infinite chain against time t, for N up to 4. Left: odd N . Right: even N . The blue straight lines have the predicted slopes 3, 6, 13 and 20.

2

  Table 2 gives our predictions in factorized form up to N = 10. 33 3 6 5 2 π 5 t 31 6 2 16 π 6 t 18 2 52 3 8 5 4 π 6 t 42 7 2 23 3 2 π 7 t 25 2 71 3 12 5 6 7 2 π 7 t 57 8 2 32 3 4 π 8 t 32 2 96 3 16 5 8 7 4 π 8 t 72 9 2 45 3 6 π 9 t 41 2 125 3 24 5 10 7 6 π 9 t 91 10 2 60 3 8 π 10 t 50 2 160 3 32 5 12 7 8 π 10 t 110

		infinite chain	near a wall		infinite chain	near a wall
	N odd	R N (t)	R	(w) N (t)	N even	R N (t)	R N (t) (w)
		(2.32)	(3.19)		(2.30)	(3.17)
	1	1 2πt	1 2πt 3	2	1 π 2 t 2	2 4 π 2 t 6
	3	2 π 3 t 5	2 9 3 2 π 3 t 13	4	2 4 π 4 t 8	2 20 3 4 π 4 t 20
	5	2 9 π 5 t 13

Table 2 .

 2 Asymptotic temporal decay of the return probability on the infinite chain and near a wall, for fermion numbers up to N = 10. The exact prefactors are given in factorized form.

Acknowledgments

It is a pleasure to thank Jean-Paul Blaizot, Adolfo del Campo, Bertrand Eynard and Jean-Marie Stéphan for fruitful discussions.

Appendix A. The Andréief identity

The Andréief identity

relates the integral of the product of the determinants built upon two families of functions to the determinant of their scalar products. This identity comes in many guises in various branches of the mathematical literature, concerning especially orthogonal polynomials and random matrix theory. Although it seems to appear for the first time in print in an article by Andréief in 1883 [START_REF] Andréief | Note sur une relation entre les intégrales définies des produits des fonctions Mém[END_REF], it is also associated with other names, including Cauchy-Binet, Gram and Heine. Here, it allows us to respectively derive (2.16) and (3.10) from (2.13) and (3.6).

Let us give an elementary proof of the above identity for the sake of completeness. Starting from the left-hand side, let us introduce the Leibniz expansions of the determinants:

where σ and τ are permutations acting on N symbols, and sgn σ = ±1 and sgn τ = ±1 are their signatures. We thus obtain

The integrand is now a product, and so

For fixed permutations σ and τ , let us change the index from n to m = τ n . We have then σ n = µ m , where µ = σ • τ -1 , and so sgn µ = sgn σ sgn τ . For fixed τ , the sum over σ can be replaced by a sum over µ. The sum over τ simply yields a factor N !.

We thus obtain

The sum over µ is nothing but the Leibniz expansion of the determinant given in the right-hand side of (A.1).

On the infinite chain, the expression (2.13) of the amplitude A N (t) is proportional to the left-hand side of (A.1), with

Applying the identity (A.1) yields (2.16).

On the semi-infinite chain, the expression (3.6) of the amplitude A (w) N (t) is proportional to the left-hand side of (A.1), with the same ρ(x n ) and

Applying the identity (A.1) yields (3.10).

Appendix B. Barnes' G-function

Barnes' G-function shares many common features with Euler's Γ-function. This appendix summarizes the main properties of both functions, which can be found in the Wikipedia article [START_REF] Wikipedia | The Free Encyclopedia[END_REF] or in the Digital Library of Mathematical Functions [47, Ch. 5.17].

Euler's Γ-function and Barnes' G-function are meromorphic functions in the complex plane obeying the recursion relations

with appropriate regularity conditions. When z is a positive integer, Euler's Γ-function becomes the usual factorial:

whereas Barnes' G-function becomes the 'superfactorial':