Topology optimization of on-chip integrated laser-driven particle accelerator

Yang-Fan He• Bin Sun,3• Ming-Jiang Ma • Wei Li • Qiang-You He •Zhi-Hao Cui • Shao-Yi Wang • Zong-Qing Zhao

Abstract Particle accelerators are indispensable tools in both science and industry. However, the size and cost of conventional RF accelerators limits the utility and scope of this technology.Recent research has shown that a dielectric laser accelerator (DLA) made of dielectric structures and driven at optical frequencies can generate particle beams with energies ranging from MeV to GeV at the tabletop level. To design DLA structures with a high acceleration gradient,we demonstrate topology optimization,which is a method used to optimize the material distribution in a specific area based on given load conditions, constraints,and performance indicators. To demonstrate the effectiveness of this approach,we propose two schemes and design several acceleration structures based on them. The optimization results demonstrate that the proposed method can be applied to structure optimization for on-chip integrated laser accelerators, producing manufacturable structures with significantly improved performance compared with previous size or shape optimization methods.These results provide new physical approaches to explore ultrafast dynamics in matter, with important implications for future laser particle accelerators based on photonic chips.

Keywords Laser-driven particle acceleration ∙Dielectric grating accelerator ∙Inverse Smith–Purcell effect ∙Topology optimization

1 Introduction

A particle accelerator is a device that accelerates charged particles(such as protons and electrons).The input energy that induces acceleration is achieved by applying an electric field. The field provides kinetic energy to the particles such that they can reach speeds at a significant fraction of the speed of light. Originally invented for scientific research, particle accelerators now play an important role in improving health and prosperity. They are relevant to our daily lives in a variety of applications,including cancer treatment, material analysis, and removal of harmful microorganisms from food and water [1–5].Radio frequency (RF)-powered devices are the traditional choice for acceleration elements [5–8]. However, their large size, high input power requirements, and expensive infrastructure limit their utility and accessibility to the wider scientific community. This has promoted the exploration of more compact and economical alternative technologies. In recent years, dielectric laser accelerators have attracted significant interest as a method for accelerating charged particles driven by solid-state lasers [9–17]. The use of infrared lasers to power particle accelerators fabricated using optical-scale lithography is a developing area of research,in which evanescent acceleration modes within dielectric structures can synchronize charged particles.This enables a more compact accelerator design and a much higher acceleration gradient compared to RF accelerators, which provides the opportunity to achieve higher acceleration gradients and reduce cost and size simultaneously (see [18] and the references therein).

A dielectric laser accelerator (DLA) aims to reduce the dimensions associated with RF accelerators by using a near-infrared (NIR) ultrafast femtosecond laser to drive dielectric structures and to provide a higher acceleration gradient. Additionally, dielectric materials such as silicon and silicon dioxide can survive at energy levels one to two orders of magnitude higher than the metals used in RF accelerators [19, 20]. Theoretically, the DLA can provide acceleration gradients in excess of GV/m [21–24] compared to the 30-100 MV/m achieved by RF accelerators,which are limited because of the need to avoid material damage. In various experiments, an acceleration gradient of nearly 1 GV/m has been achieved (see [25] and the references therein). Compared with other advanced acceleration techniques, DLA currently provides the highest gradient among non-plasma accelerators. Because dielectric laser accelerators consist of dielectric structure devices and are driven at optical frequencies, cascade acceleration can be easily integrated on-chip [26, 27]. To make DLA competitive with larger devices, the Accelerator on a Chip International Program (ACHIP) [28], which is working toward building miniature particle accelerators on a chip using advanced laser and nanofabrication technology, has divided research topics in this field into several branches for detailed study; however, as stated in [29], the nanostructures of DLA are particularly important, and different nanostructures can be expected to have very different effects on particles driven by the same field distribution.

Most optimization studies and experiments in the field of DLA are based on grating structures [29]. Additionally,related theories have been gradually developed to explain the acceleration phenomenon. However, these relevant theories either do not consider how the structural parameters affect the particle dynamics or are limited to the grating structure and how its local parameters affect the particle dynamics [29]. In general, past research has focused on shape and size optimization, with examples shown in Fig. 1a and b. Obviously, the results obtained in this manner are limited by the initial architecture, and the results obtained are a subset of various possibilities depending on the structure used. A pioneering study on topology optimization ( Fig. 1c) of DLA structures was performed using adjoint methods by Hughes et al. [30].When using the adjoint method for topology optimization,the adjoint source must be given [31, 32]. In principle, the adjoint source is the first derivative of the objective function with respect to the design parameters.This means that there are situations in which this derivative cannot be solved for, so it is not possible to provide an accurate accompanying source for simulation purposes. In other words, it may be difficult to use adjoint methods to implement the inverse design of arbitrary objective functions. Additionally, to produce structures that can be practically fabricated, a mathematical method that can generate additional binary structure distributions is needed.Therefore, an effective method for designing an acceleration structure that can be manufactured is highly desirable.

This study was inspired by the utilization of a self-guided derivation technique used in artificial intelligence to compute gradients for an arbitrary objective function, and applies a hyperparameter control function distribution to guide the generation of a more binary structure distribution to achieve more general topology optimization. In this paper, we outline in detail how topology optimization can be used to achieve the inverse design of DLA structures.The remainder of this paper is organized as follows. In Sect. 2,we first outline methods for simulating the physical phenomena associated with DLA,and then demonstrate the optimization principle for realizing a DLA structure using automatic derivation techniques. Finally, we present a control method that can realize a manufacturable structure.In Sect. 3, the validity of the proposed method is verified first, two optimization schemes that can be used in DLA are discussed, and design demonstrations of these two optimization schemes are provided.

2 Theoretical descriptions

All macroscopic electromagnetism, including the propagation of light in a photonic crystal, is governed by Maxwell’s four macroscopic equations.In SI units,we seek vector fields D,E,B,H : Ω →R3such that the Maxwell equations

hold for a given J :Ω →R3and ρ:Ω →R.We call E the electric field, D electric displacement field (or electric flux density), H magnetic field intensity, B magnetic flux density, and J and ρ electric current intensity and electric charge density, respectively. Although the results obtained by analytically solving Maxwell’s equations have theoretical significance, when combined with the research object of this study, it is very difficult to obtain their analytical solutions. Therefore, only numerical solutions of Maxwell’s equations are considered in this study,which can be divided into two categories: time domain (TD) methods and frequency domain (FD) methods. The most famous time domain method is the finite difference time domain(FDTD) [33], which is usually suitable for solving for the transient change process under external excitation. The frequency domain method is based on time harmonic differential and integral equations to study the steady-state field distribution after infinite time under a time-harmonic excitation. Considering the small memory and low time consumption requirements, this study uses the finite-difference frequency-domain (FDFD) approach to solve for the field distribution of the structure at a given frequency(ω)for a DLA coupled with non-magnetic materials under normal conditions. For mathematical convenience, we employ the standard trick of using a complex-valued field and taking the real part to obtain the physical fields. This allows us to write a harmonic mode as a spatial pattern(or‘‘mode profile’’) multiplied by a complex exponential as follows.

where E(r)is a complex vector representing the(unknown)electric field distribution in the domain. Generally, ignoring higher-order relationships, a dielectric medium with permeability μ=μ0μrand permittivity ∊=∈0∊ris considered,where μ0is the vacuum permeability,and ∈0is the vacuum permittivity. The working principle of FDFD[34,35]is to represent the electric field,dielectric constant distribution, and source on the Yee grid [33], and then represent the ∇×∇× operator as a finite-difference differential matrix ～A. In the discrete case, Eq. (1) can be written in matrix form as

In 2D problems, the electric and magnetic fields can be separated into two orthogonal modes:TM and TE ,where the TM mode can be fully described by Hzand the TE mode by Ez,as shown in Fig. 2.It is worth mentioning that the simulation area of the DLA cannot be infinite, that is,the simulation boundary must be truncated, and the commonly used condition is the perfect matching layer and periodic boundary (see [36] and the references therein).The process will be discussed in more detail in Sect. 3 for specific examples.

To implement topology optimization,the DLA structure is divided into three parts: laser injection region (LIR),optimization design region (ODR), and particle injection channel region (PICR), as shown in Fig. 3a. To obtain the device structure under the optimal objective function ζ(xi)value, a straightforward idea is to mesh the design area such that the permittivity value in each grid is selected at random as 0 or 1, where 0 represents vacuum and 1 represents the permittivity constant of the material, and then change the permittivity value state of one grid, run the simulation with the permittivity value of the other grid unchanged, record the objective function ζ(x) value and repeat the calculations until all the combined results are traversed, and finally select an optimal structure from all simulation results, as shown in Fig. 3b.

However,this enumeration method will become difficult to implement in practice because even when enumerating only 10 grids in the optimization design region, it must berun 210times.This also means that with further refinement of the grid in the optimization design region,the number of calculations will increase exponentially.In other words,the idea is straightforward,but not practical.Therefore,we use the gradient of the objective function to update and optimize the parameters to minimize the objective function ζ(xi). To make the computation of the derivative of the objective function ζ(xi)with respect to the design variables more flexible and extensive, we employ an automatic derivation technique.Automatic derivation is a method that lies between symbolic derivation and numerical derivation(it is a numerical calculation method that calculates an approximation of the derivative). In general, automatic derivation is the application of a symbolic derivation to the most basic operations, such as constant, power function,exponential function, logarithmic functions, trigonometric function, and other basic function. The value of the independent variable is substituted to obtain its derivative value,which is retained as an intermediate result.Then,the derivative value of the entire function is calculated according to the derivation results of these basic operation units. The difference between automatic derivation and symbolic derivation is that instead of computing the analytical solution, the chain rule is used to decompose the complex function into individual units, differentiate these small units,and finally combine them to obtain the solution to the full derivative.

Fig. 2 (Color online) Fundamental features of DLA used in the simplified 2D analysis, in which the full-field variation with z is not considered. In 2D problems, the electric and magnetic fields can be separated into orthogonal TM and TE modes, where the TM mode can be fully described by Ez,and the TE mode can be fully described by Hz.The field distribution of the TM(TE)mode is discretized onto a square Yee grid, as shown in a and b

Fig.3 (Color online)In our topology optimization method,the DLA structure is divided into three parts: laser injection region (LIR),optimization design region (ODR) and particle injection channel region(PICR),as shown in a.The permittivity function of the design area is randomly distributed as 0 (white block) or 1 (dark block),where 0 represents vacuum,and 1 represents the permittivity function of the material, as shown in b

Fig. 4 (Color online) An objective function is divided into a calculation graph consisting of multiple calculation units obtained according to the calculation order. This directed acyclic graph charts the flow of data

The objective function ζ(xi) is divided into multiple computing units obtained according to the order of computation, which we call a computation graph, which is a directed acyclic graph that charts the flow of data,as shown in Fig. 4. The purpose is to find the derivative∇ζ(xi)=[∂ζ/∂x1, ∂ζ/∂x2], so we apply the chain rule to the first term, ∂ζ/∂x1, which yields

Fig. 5 (Color online) We show two specific modification methods using Θ1 and Θ2, where Θ1 is a ‘‘Sigmoid function,’’ and Θ2 is a‘Tanh function.’’ SΘ1 and SΘ2 are the regions corresponding to the non-extreme points of the function,and the size of this region can be controlled by introducing one or more hyperparameters. That is, the function can be made more binary by introducing certain hyperparameters a

We require only one forward sweep to obtain the result.The calculation graph from left to right corresponds to Eq. (4) from inside to outside. If we use the chain rule in another way, such as using Eq. (5), the method needs to first calculate the rightmost derivative and then calculate the derivative to its left in turn. Therefore, not only a forward propagation, but a reverse sweep is required. In automatic derivation,the first method is called the forward mode and the second is called the reverse mode. An automatic differentiation library allows the user to compute the exact Jacobian dxi/dxjfor any i, j using the rules of differentiation and some knowledge of the partial derivatives of each operation (see [37, 38] and the references therein). When the objective function is convex, the solution of the gradient descent method is globally optimal.However, in general, the objective function is often not a convex function; therefore, it cannot be guaranteed that its solution is the global optimal solution. To prevent the solution from falling into a local optimum, adaptive moment estimation [39] is used in our study.

To obtain a device that can actually be manufactured or meaningful, we need to modify the equipment parameterization to encourage the use of more binary features in the optimization,because if the optimized equipment is simply represented as a set of data, the eigenvalues during optimization will change constantly, and the final output may not produce a device that can be manufactured.In general,we can modify the parameters into the following form:

where Θ is a common s-shaped function, as shown in Fig. 4.We demonstrate two specific modification methods involving Θ1and Θ2, respectively, where Θ1is a ‘‘Sigmoid-like function’’,and Θ2is a‘‘Tanh-like function’’,and a represents the hyperparameters. SΘ1and SΘ2are the regions corresponding to the non-extremum points of the function,and the size of either region can be controlled by introducing hyperparameters, that is, the function can be made more binary by introducing certain hyperparameters a.Alternatively,Heaviside’s step function(κ,Eq. (8))with projection significance can be used.

Here, σ is the projection threshold (when the value is less than σ,and the Heaviside function(κ)projects the value of the unit in the 0 direction. In contrast, when the density of the unit is greater than σ, the Heaviside function projects the value of the unit in the direction of 1), where τ is a set of hyperparameters that controls the steepness around the threshold parameter σ. In the actual process, convolution filtering is required to smooth the edges of the design structure [32].

3 Results and discussion

The two important parameters (∊and μ in Eq. (1)) that affect the electromagnetic distribution can be written in the following form:

Now, consider a material that consists of infinite columns in the x3-direction modeled by a permittivity such that (x1,x2)∊(x1, x2). We assume an isotropic linear dielectric medium with permeability μ=μ0μrand permittivity∊=∈0∊r. However, only isotropic non-magnetic materials are considered in this study, which are the preferred materials for chip accelerators. Their properties can be described by

The coordinate system used for DLA optimization is shown in the figure; therefore, the TE mode should be used to realize particle acceleration. At this time, the electromagnetic distribution of the system can be calculated using Eq. (11).

The terms ∇Exand ∇Eyare banded matrices that calculate first order spatial derivatives of the electric fields across the grid. As a quick example of these matrices, they were computed for a small two-dimensional grid composed of 3×3 cells.Using periodic boundary conditions,the matrix derivative operators ∇Exand ∇Eyfor this simple case are

The ‘‘HT’’ superscript indicates a Hermitian transpose operation. For the TE mode, averaging is needed for ∈11and ∈22.All involved matrices are sparse;therefore,we can apply sparse matrix techniques to save computation time and memory.

The acceleration gradient (or the particle energy gain per unit length) is an important figure of merit. For a laser operating in normal incidence mode (laser propagating in the y direction), phase velocity matching between the particle and the electromagnetic fields is established by introducing a spatial periodicity in our structure having a period of βλ along the x direction, where λ is the laser wavelength. The acceleration gradient function Ξgain(E)can be calculated using a periodic structure as follows.

Here, we focus only on the relationship between the permittivity distribution (∊r) and the energy gain (Ξgain(E)).In other words, this is exactly the same case for nonmagnetic materials mentioned earlier. As described in Sect. 2,we must compute Ξgain(E) for ∊ras a derivative of

where Ξgain(E) is equivalent to the objective function(ζ(xi)) mentioned previously. ∂Ξgain(E)/∂∊ris an important parameter that indicates how to change ∊rto increase Ξgain(E). Two optimization schemes for maximizing the acceleration gradient function Ξgain(Ey) are presented below.

Case 1: Considering the energy conversion rate, the acceleration gradient should be higher than the peak value of the incident electric field as much as possible. This means that the ‘‘conversion factor’’ (Γ) given by dividing the acceleration gradient by the peak amplitude of the incident electric field needs to be optimized. This quantity will reveal the maximum conversion rate that the structure can achieve.Using mathematical language,we express this principle as

where E0is the initial injection electric field,which should be less than the material damage threshold(Ematerial damage).

Case 2: Specifically considering the protection or reuse of the structure,the optimal energy gain should be selected so that the acceleration gradient is as high as possible without causing damage to the structure. This implies that the value of the maximum electric field in the structure must be determined so that it cannot exceed the local peak laser fluence at the damage threshold. Therefore, another relevant quantity to maximize is the ‘‘acceleration factor’’(χ) given by the acceleration gradient divided by the maximum electric field magnitude in the system. This amount ultimately limits the amount of acceleration gradient that can be achieved so that the damage threshold is not exceeded. Using mathematical language, we express this principle as Here, Emax is the maximum value of the electric field magnitude vector in the optimization design region (ODR)and Ematerial damage is the material damage threshold expressed in terms of the electric field. The damage threshold of the DLA structure can be improved by the choice of material.

Fig. 6 (Color online) The parameter settings are the same as those used in Ref. [30], The material used in the optimization process is fused silica (SiO2). A plane wave (E0 is the initial injected electric field) is introduced into the chip from the left side using the periodic condition. The laser wavelength (λ) is 2000 nm, The normalized velocity (β) of the injected electrons is 0.5, ε11 =ε22 =2.1, the square grid size (Δx(y)) used in the simulation is 20 nm, and the achieved gradient is 0.32E0 (Γ = 0.32)

Fig. 7 (Color online). Topology optimization results of DLA structure using ‘Case 1’’ scheme. The parameters used are the same as those used in[29].A plane wave(E0 is the initial injected electric field) is introduced into the chip from the left side using the periodic condition. The laser wavelength (λ) is 800 nm, the normalized velocity (β) of the injected electrons is 0.99, and the square grid size(Δx(y)) used in the simulation is 8 nm

Fig.8 (Color online)Topology optimization results of DLA structure using‘Case 2’scheme.The parameter settings used are the same as those used in [29]

The optimization results show that the optimized structures obtained by the two schemes (‘‘Case 1’’, and‘‘Case 2’’) are significantly different from the intuitive structures, and it is difficult to explain these results using physical intuition.It is worth noting that in the comparison between ‘‘Case 1’’ and ‘‘Case 2,’’ the structures obtained by the two optimization schemes are completely different.If χ is used as the evaluation index, the result for the acceleration gradient function in Fig. 7 is three times that in Fig. 8.If Γ is used as the evaluation index,the result for the acceleration gradient function in Fig. 8 is twice that of Fig. 7. The acceleration gradients of the two schemes are different, but they are significantly improved when compared to the grating structure [29]. A phenomenological explanation is that in the y-direction of the topologically optimized DLA structures (Figs. 7b and Fig. 8b), the electric field (Ey/|E0|) exhibits a distinct standing wave field pattern in the particle injection channel region(PICR),which is very favorable for the acceleration of charged particles. This also indicates that the DLA optimized by either of these two schemes can pass more particle beams,which undoubtedly can increase particle excitation (by appropriately increasing the width of the particle injection channel region). Although increasing the width of the particle injection channel region (PCIR) reduces the acceleration gradient of the DLA, the optimized structures obtained by these two schemes provide an intrinsic acceleration gradient that is several times greater than that achieved by the grating-structured accelerator [29].

4 Summary

The topology optimization algorithm presented in this study can realize the optimization of any objective function, and the obtained structure also has a high-contrast binary structure. The advantage of this method is that the structure can be designed according to the expected goal,so that the optimal geometry can be found more intelligently without tedious adjustment of the associated parameters. By changing the material parameterization function, we can transition to shape optimization or size optimization. We study the acceleration structure of an accelerator based on the high-gradient and compact characteristics of the DLA. The two design schemes presented ignore the influence of Exthat the particle feels. However,the effect of Exwill lead to the deflection of accelerated particles over long distances; therefore, it may require that the design include other structures embedded in the junction of the accelerating structure to focus the beam in order to prevent the effect of Exin the material.

In future work, our goal is to design a complete dielectric laser accelerator that supports larger sizes and increased focus effects. To achieve simple manufacturing,we may need to add more constraints to the current parameterization.

This technique is also applicable to the design of other dielectric-based accelerator structures.This is based on the inverse Cherenkov Effect of acceleration. As such, our work provides an opportunity to significantly increase the energy gain of DLAs, which is crucial for the practical application of these exciting micro accelerators.

AcknowledgementsThe authors thank Dr.Fu-Long Liu of the China Institute of Atomic Energy, Dr.Wen-Bo Mo of Tsinghua University,Dr. Dan Wang of Sichuan University, and Dr. Han-Sheng Ye of Tsinghua University for the illuminating discussions. We thank the reviewers for their constructive comments.

Author ContributionsAll authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Yang-Fan He and Bin Sun. The first draft of the manuscript was written by Yang-Fan He and all authors commented on previous versions of the manuscript.All authors read and approved the final manuscript.