Method for Computationally Efficient Design of Dielectric Laser Accelerators

Dielectric microstructures have generated much interest in recent years as a means of accelerating charged particles when powered by solid state lasers. The acceleration gradient (or particle energy gain per unit length) is an important figure of merit. To design structures with high acceleration gradients, we explore the adjoint variable method, a highly efficient technique used to compute the sensitivity of an objective with respect to a large number of parameters. With this formalism, the sensitivity of the acceleration gradient of a dielectric structure with respect to its entire spatial permittivity distribution is calculated by the use of only two full-field electromagnetic simulations, the original and adjoint. The adjoint simulation corresponds physically to the reciprocal situation of a point charge moving through the accelerator gap and radiating. Using this formalism, we perform numerical optimizations aimed at maximizing acceleration gradients, which generate fabricable structures of greatly improved performance in comparison to previously examined geometries.


Introduction
Dielectric laser accelerators (DLAs) are periodic dielectric structures that, when illuminated by laser light, create a near-field that may accelerate electrically charged particles such as electrons [1]. A principal figure of merit for these DLA structures is the acceleration gradient, which signifies the amount of energy gain per unit length achieved by a particle that is phased correctly with the driving field. DLAs may sustain acceleration gradients on the order of ∼GV m −1 when operating using the high peak electric fields supplied by ultrafast (femtosecond) lasers. These acceleration gradients are several orders of magnitude higher than conventional particle accelerators. As a result, the development of DLA can lead to compact particle accelerators that enable new applications.
In previous works, candidate DLA geometries were optimized for maximum acceleration gradient by scanning through parameters of a specified structure geometry [2][3][4][5][6][7][8][9]. However, this strategy has limited potential to produce higher acceleration gradient structures because it only searches a small portion of the total design space.
In this paper, we derive an analytical form for the sensitivity of the acceleration gradient of a DLA structure with respect to its permittivity distribution using the adjoint-variable method (AVM). We may calculate this by use of only two full-field simulations. The first corresponds to the typical accelerator setup, where the structure is illuminated with externally incident laser light. The second corresponds to the inverse process, where the same physical structure is simulated but now with a charged particle traversing the structure as the source. Thus, this formalism explicitly makes use of the reciprocal relationship between accelerators and radiators [10,11]. We use this sensitivity information to perform optimizations, which generate DLA structures of much higher gradients than previously explored geometries.
This work is the first application of the AVM technique to the design of DLA structures and gives examples of fabricable structures that may improve the energy gain achievable with current DLA technology. In addition, the optimized structures give insight into general design principles for DLAs, meaning that one may use the principle findings of this paper to design DLAs without having to run optimizations directly. As an example, it was found that high gradient structures often include dielectric mirrors surrounding the particle gap, leading to higher field enhancement.
This paper is organized as follows: We first outline the status of DLAs and basic design requirements in section 2. We introduce AVM in section 3, where we derive the sensitivity of the acceleration gradient of a DLA with respect to its permittivity distribution. In section 4, we show that the 'adjoint' solution corresponds to that of a radiating charge. In section 5, we describe and demonstrate algorithms for using the sensitivity information to design DLA structures numerically.

A Brief Review of Dielectric Laser Accelerators
DLAs take advantage of the fact that dielectric materials have high damage thresholds at short pulse durations and infrared wavelengths [1,4,12] when compared to metal surfaces at microwave frequencies. This allows DLAs to sustain peak electromagnetic fields, and therefore acceleration gradients, that are 1 to 2 orders of magnitude higher than those found in conventional radio frequency (RF) accelerators. Experimental demonstrations of these acceleration gradients have been made practical in recent years by the availability of robust nanofabrication techniques combined with modern solid state laser systems [13]. By providing the potential for generating relativistic electron beams in relatively short length scales, DLA technology is projected to have numerous applications where tabletop accelerators may be useful, including medical imaging, radiation therapy, and X-ray generation [1,14]. To achieve high energy gain in a compact size, it is of principle interest to design structures that may produce the largest acceleration gradients possible without exceeding their respective damage thresholds.
Several recently demonstrated candidate DLA structures consist of a planar dielectric structure that is periodic along the particle axis with either an semi-open geometry or a narrow (micron to sub-micron) vacuum gap in which the particles travel [2][3][4][5][6][7][8][9]. These structures are then sideilluminated by laser pulses. Fig. 1 shows a schematic of the setup, with a laser pulse incident from the bottom.
The laser field may also be treated with a pulse front tilt [15,16] to enable group velocity matching over a distance greater than the laser's pulse length. For acceleration to occur, the dielectric structure must be designed such that the particle feels an electric field that is largely parallel to its trajectory over many optical periods. In the following calculations, the geometry of the dielectric structure is represented by a spatially varying dielectric constant (x, y). We assume invariance in one coordinate (ẑ) in keeping with the planar symmetry of most current designs. However the methodology we present can be extended to include a third dimension. In eε(x,y) x y βc 0 βλ λ Fig. 1. Diagram outlining the system setup for side-coupled DLA with an arbitrary dielectric structure (x, y) (green). A charged particle moves through the vacuum gap with speed βc 0 . The periodicity is set at βλ where λ is the central wavelength of the laser pulse. addition, our work approximates the incident laser pulse as a monochromatic plane wave at the central frequency, which is a valid approximation as long as the pulse duration is large compared to the optical period.

Adjoint Variable Method
In a general DLA system, we may define the acceleration gradient 'G' over a time period 'T' mathematically as follows: where ì r(t) is the position of the electron and E | | signifies the (real) electric field component parallel to the electron trajectory at a given time.
To maximize this quantity, we employ AVM [17,18], which is a technique common to a wide range of fields, including seismology [19], aircraft design [20], and, recently, photonic device design [21][22][23]. Many engineering systems can be described by a linear system of equations A(γ)z = b, where γ is a set of parameters describing the system. For a given set of parameters γ, solving this equation results in the solution 'z', from which an objective J = J(z), which is a function of the solution, can be constructed. The optimization of the engineering system corresponds to maximizing or minimizing J with respect to the parameters γ. For this purpose, AVM allows one to calculate the gradient of the objective function ∇ γ J for an arbitrary number of parameters γ i with the only added computational cost of solving one additional linear system A Tz = − dJ dz T , which is often called the 'adjoint' problem. For a more comprehensive overview of the method, we refer the reader to [17].
Here we provide the derivation of AVM specifically for the optimization of the accelerator structures. Since the structure is invariant in theẑ direction, we work in two dimensions, examining only the H z , E x and E y field components. For an approximately monochromatic input laser source with angular frequency ω, the electric fields are, in general, of the form where now ì E is complex. Let us assume the particle we wish to accelerate is moving on the line y = 0 with velocity ì v = βc 0x , where c 0 is the speed of light in vacuum and β ≤ 1. The x position of the particle as a function of time is given by x(t) = x 0 + βc 0 t, where x 0 represents an arbitrary choice of initial starting position. For normal incidence of the laser (laser propagating in the +ŷ direction), phase velocity matching between the particle and the electromagnetic fields is established by introducing a spatial periodicity in our structure of period βλ alongx , where λ is the laser wavelength. In the limit of an infinitely long structure (or equivalently, T → ∞) we may rewrite our expression for the gradient in Eq. (1) as an integral over one spatial period, given by Here the quantity φ 0 = 2π x 0 βλ is representative of the phase of the particle as it enters the spatial period. In further calculations, we set φ 0 = 0, only examining the acceleration gradients experienced by particles entering the accelerator with this specific phase. Since we have arbitrarily control over our input laser phase, this does not impose any constraint on the acceleration gradient attainable.
To simplify the following derivations, we define the following inner product operation involving the integral over two vector quantities ì a and ì b over a single period volume V With this definition, we then have the gradient where Now, we wish to examine the sensitivity of G with respect to an arbitrary parameter, γ, which may represent a shifting of material boundary, changing of dielectric constant at a point, or any other change to the system. Differentiating Eq. (5) gives Here we have made use of the fact that ì η does not depend on γ. From Maxwell's equations in the frequency domain, we may express our electromagnetic problem in terms of a linear operatorÂ as Here, k 0 = ω/c 0 , r is the relative permittivity, ì J represents a current density source, and a non-magnetic material is assumed (µ = µ 0 ). Differentiating Eq. (8) with respect to γ, and assuming that the current source ( ì J) does not depend on γ, we see that A is self-adjoint under our inner product, Â ì a, ì b = ì a,Â ì b , and the same is true forÂ −1 and dÂ dγ . Using these facts and combining Eq. (7) with Eq. (9), we find that Thus, if we define a second simulation with a source of − ì η and fields ì E a j , then the field solution, ì E a j = −Â −1 ì η, can be easily identified in Eq. (10). The sensitivity of the acceleration gradient can finally be expressed as The only quantity in this expression that depends on the parameter γ is dÂ dγ . As we will soon discuss, this quantity will generally be trivial to compute. On the other hand, the full field calculations of ì E and ì E a j are computationally expensive, but may be computed once and used for an arbitrarily large set of parameters γ i . This gives the AVM approach a significant scaling advantage with respect to traditional direct sensitivity methods, which require a separate full-field calculation for each parameter being investigated. It is this fact that we leverage with AVM to do efficient optimizations over a large design space.
To confirm that this derivation matches the results obtained by direct sensitivity analysis, we examine a simple accelerator geometry composed of two opposing dielectric squares each of relative permittivity . We take a single γ parameter to be the relative permittivity of the entire square region. Because we only change the region inside the dielectric square, we may identify the dÂ dγ operator for this parameter by examining Eq. (8), giving Thus, given the form of the acceleration gradient sensitivity given in Eq. (12), combined with Eq. (13), the change in acceleration gradient with respect to changing the entire square permittivity is simply given by the integral of the two field solutions over the square region, labeled 'sq' In Fig. 2 we compare this result with the direct sensitivity calculation where the system is manually changed and simulated again. The two methods agree with excellent precision. Extending this example to the general case of perturbing the permittivity at an arbitrary position ì r , we see that

Reciprocity
With the AVM form derived, we now wish to re-examine the adjoint source term from Eq. (11) in another interpretation. Let us now consider the fields radiated by a point particle of charge q flowing through our domain at y = 0 with velocity ì v = βc 0x . In the time domain, we can represent the current density of this particle as We may take the Fourier transform of ì J r ad with respect to time to examine the current density in the frequency domain, giving Comparing with the source of our adjoint problem, ì J a j = −i ωµ 0 ì η, we can see that This finding shows that the adjoint field solution ( ì E a j ) corresponds (up to a complex constant) to the field radiating from a test particle flowing through the accelerator structure. To put this another way, in order to calculate the acceleration gradient sensitivity with AVM, we must simulate the same structure operating both as an accelerator (Â ì E = −iωµ 0 ì J acc ) and as a radiator . It is understood that one way to create acceleration is to run a radiative process in reverse. Indeed, this is the working principle behind accelerator schemes such as inverse free electron lasers [24,25], inverse Cherenkov accelerators [26,27], and inverse Smith-Purcell accelerators [28,29]. Here, we see that this relationship can be expressed in an elegant fashion using AVM.

Finite-Difference Frequency-Domain Modeling
Now that we have shown how to use AVM to compute the sensitivity of the acceleration gradient with respect to the permittivity distribution, we will show practical applications of these results. First, for computational modeling, the problem must be transitioned from a continuous space to a discrete space. Here we make the transition using a finite-difference frequency-domain (FDFD) formalism [30,31]. The electromagnetic fields now exist on a Yee lattice and the linear operatorÂ becomes a sparse, complex symmetric matrix, A, relating the vector of electric field components, e, to the input current source components b as To solve for the field components, this system must be solved numerically for e. In twodimensions, this is usually done directly by use of "lower-upper" (LU) decomposition methods for sparse matrices. Only the right hand side of Eq. (22) is different between the original and adjoint simulations. Therefore after the A matrix is factored to solve the original simulation, its factored form may be saved and reused for the adjoint calculation, which cuts the total computational running time roughly in half.
Written in terms of this discrete system, the acceleration gradient is where η is now a discretized version of ì η. Similarly, the sensitivity of the gradient with respect to changing the permittivity at pixel 'i' is given by where, as before,ē is the solution of the adjoint problem For all simulations, we use an FDFD program developed specifically for this work, although a commercial package would also be sufficient. We have chosen a grid spacing that corresponds to 200 grid points per free space wavelength in each dimension. Perfectly matched layers are implemented as absorbing regions on the edges parallel to the electron trajectory, with periodic boundary conditions employed on boundaries perpendicular to the electron trajectory. A total-field scattered-field [31] formalism is used to create a perfect plane wave input for the acceleration mode.

Gradient Maximization
Since we now have a highly efficient method to calculate dG d i , we proceed to use this information to maximize the acceleration gradient with respect to the permittivity distribution. We use an iterative algorithm based on batch gradient ascent [32]. During each iteration, we first calculate dG d i for all pixels 'i' within some specified design region. Then, we update each i grid as follows Here, α is a step parameter that we can tune. We need α to be small enough to find local maxima, but large enough to have the optimization run in reasonable amount of time. This process is then iterated until convergence on G. During the course of optimization, the permittivity distribution is considered as a continuous variable, which is not realistic in physical devices. To address this issue, we employ a permittivity capping scheme during optimization. We define a maximum permittivity ' m ' corresponding to a material of interest. During the iterative process, if the relative permittivity of any cell becomes either less than 1 (vacuum) or greater than m , that cell is pushed back into the acceptable range. It was found that with this capping scheme, the structures converged to binary (each pixel being either vacuum or material with a permittivity of m ) after a number of iterations without specifying this choice of binary materials as a requirement of the optimization. Therefore, only minimal post-processing of the structures was required. The results of this optimization scheme are shown in Fig. 3(b-d) for three different m values corresponding to commonly explored DLA materials. The design region was taken to be a rectangle fully surrounding but not including the particle gap. The design region was made smaller for higher index materials, since making it too large led to divergence during the iteration. We found that a totally vacuum initial structure worked well for these optimizations. However, initially random values between 1 and m for each pixel within the design region also gave reasonable results.
This optimization scheme seems to favor geometries consisting of a staggered array of fieldreversing pillars surrounding the vacuum gap, which is already a popular geometry for DLA. However, these optimal designs also include reflective mirrors on either side of the pillar array, which suggests that for strictly higher acceleration gradients, it is useful to use dielectric mirrors to resonantly enhance the fields in the gap. It was observed that for random initial starting permittivity distributions, the same structures as shown in Fig. 3 are generated every time. Furthermore, these geometries are remarkably similar to those recently proposed [33], although these designs do not include the reflective front mirror. These findings together suggest that the proposed structures may be close to the globally optimal structure for maximizing G.
It was further found that convergence could be achieved faster by a factor of about ten by including a 'momentum' term in the update equation. This term corresponds to the sensitivity calculated at the last iteration multiplied by a constant, α < 1. Explicitly, for iteration number ' j' and pixel 'i'

Acceleration Factor Maximization
DLAs are often driven with the highest input field possible before damage occurs. Therefore, another highly relevant quantity to maximize is the 'acceleration factor' given by the acceleration gradient divided by the maximum electric field amplitude in the structure. This quantity will ultimately limit the amount of acceleration gradient we can achieve when running at damage threshold. Explicitly, the acceleration factor is given by Here, | ì E | is a vector of electric field magnitudes in our structure, and the max{} function is designed to pick out the highest value of this vector in either our optimization or material region, depending on the context. We would like to use the same basic formalism to maximize f A . However, since the max{} function is not differentiable, this is not possible directly. Instead we may use a 'smooth-max' function to approximate max{} as a weighted sum of vector components Here, the parameter a ≥ 0 controls the relative strength of the exponential sum terms, for a = 0, this function simply gives the average value of the field amplitudes. By sweeping a and examining the acceleration factors of the resulting optimized structures, we determined that a = 3 gave the best improvement in f A . If a is too large, the calculation may induce floating point overflow or rounding error issues. Using this smooth-max function, one may calculate d f A d i analytically and perform structure optimizations in the same way that was discussed previously. The derivation of the adjoint source term is especially complicated and omitted for brevity, although the end result is expressed solely in terms of the original fields, the adjoint fields, and the dÂ dγ operator, as before. Two structures with identical parameters but optimized, respectively, for maximum G and f A are shown in Fig.  4. On the left, we see that the G maximized structure shows the characteristic dielectric mirrors, giving resonant field enhancement. On the right is the structure optimized for f A , which has eliminated most of its dielectric mirrors and also introduces interesting pillar shapes. In Table 1 the main DLA performance quantities of interest are compared between these two structures. Whereas the acceleration gradient is greatly reduced when maximizing for f A , the f A value itself is improved by about 25% or 23% depending on whether one measures the maximum field in the design region or the material-only region, respectively. These findings suggest that the AVM strategy is effective in designing not only resonant, high acceleration gradient structures, but also non-resonant structures that are more damage resistant. In the future, when more components of DLA are moved on-chip (such as the optical power delivery), it will be important to have control over the resonance characteristics of the DLA structures to prevent damage breakdown at the input facet. Our technique may be invaluable in designing structures with tailor-made quality factors for this application.

Discussion
We found that AVM is a reliable method for optimizing DLA structures for both maximum acceleration gradient and also acceleration factor. The optimization algorithm discussed shows good convergence and rarely requires further post-processing of structures to create binary permittivity distributions. Therefore, it is a simple and effective method for designing DLAs. Whereas most structure optimization in this field uses parameter sweeps to search the design space, the efficiency of our method allows us to more intelligently find optimal geometries without   [5] show gradients far below those presented here. We had limited success designing DLA structures in the relativistic (β ≈ 1) regime, especially for higher index materials, such as Si. We believe this is largely due to the stronger coupling between electron beam and incident plane wave at this energy. The characteristics of the adjoint source change dramatically at the β = 1 point. Whereas in the sub-relativistic regime, the adjoint source generates an evanescent near-field extending from the gap particle position, at β ≥ 1, the adjoint fields become propagating by process of Cherenkov radiation. Upon using the above described algorithm, the gradients diverge before returning to low values, no matter the step size α. The only way to mitigate this problem is to choose prohibitively small design regions or low index materials, such as SiO 2 .
The AVM formalism presented in this work may also be extended to calculate higher order derivatives of G. For each higher order, the form of the derivative of G can be derived in a fashion very similar to the one outlined for first order. Given the full Hessian H i, j = d 2 G d i d j , as calculated by AVM, one could use Newton's method to do optimization. However, to perform exactly, this calculation would require as many additional simulations as there are grid points within the design region. Therefore, these higher order methods are inconvenient for our purposes where there are generally tens of thousands of design pixels. This limitation may be averted by using approximate methods for finding the inverse Hessian [34], which may provide substantial improvement to optimization results and convergence speeds in certain cases. However, in our case there was no need to explore beyond first order due to the relative success and speed of the algorithm presented.
As future works, our goal is to fabricate and test these structures experimentally, as well as include further metrics into the optimization if necessary, such as favoring larger feature sizes and incorporating focusing effects. Furthermore, this method is of great interest in designing waveguide-coupled accelerator structures, where typical designs optimized for plane wave input are likely suboptimal. This will be of critical importance when moving the optical power delivery source on-chip.
In addition to the side-incident geometry explored, this technique is applicable to designing other dielectric-based accelerator structures. This includes particle-laser co-propagating schemes [35] and perhaps dielectric wakefield acceleration [36], among others. Therefore, we expect that our results may find use in the larger advanced accelerator community.

Conclusion
We have introduced the adjoint variable method as a powerful tool for designing dielectric laser accelerators for high gradient acceleration and high acceleration factor. We have further shown that the adjoint simulation is sourced by a point charge flowing through the accelerator, which quantifies the reciprocal relationship between an accelerator and a radiator.
Optimization algorithms built on this approach allow us to search a substantially larger design space and generate structures that give gradients far above those normally used for each material. Furthermore, the structures designed by AVM are fundamentally not constrained by shape parameterization, allowing never-before-seen geometries to be generated and tested.

Funding
This work was supported by the Gordon and Betty Moore Foundation under grant GBMF4744 (Accelerator on a Chip), and by the U.S. Department of Energy under Contract DE-AC02-76SF00515.