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

Topology optimization of on-chip integrated laser-driven particle accelerator 1 2 corresp first-author y-f-he@foxmail.com y-f-he@foxmail.com y-f-he@foxmail.com https://orcid.org/0000-0001-7172-9403 1 2 3 https://orcid.org/0000-0001-9242-3454 4 1 2 1 2 1 2 1 2 1 corresp zhaozongqing99@caep.ac.cn zhaozongqing99@caep.ac.cn zhaozongqing99@caep.ac.cn Laser Fusion Research Center, CAEP, P. O. Box 919 986, Mianyang 621900, China Laser Fusion Research Center, CAEP, P. O. Box 919 986, Mianyang 621900, China The Sciences and Technology on Plasma Physics Laboratory, CAEP, Mianyang 621900, China The Sciences and Technology on Plasma Physics Laboratory, CAEP, Mianyang 621900, China Department of Plasma Physics and Fusion Engineering, CAS Key Laboratory of Geospace Environment, University of Science and Technology of China, Hefei 230026, China Department of Plasma Physics and Fusion Engineering, CAS Key Laboratory of Geospace Environment, University of Science and Technology of China, Hefei 230026, China The Institute of Optics, University of Rochester, Rochester, NY 14627, USA The Institute of Optics, University of Rochester, Rochester, NY 14627, USA 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. Laser-driven particle acceleration Dielectric grating accelerator Inverse Smith-Purcell effect Topology optimization Introduction 1 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 1 - 5 5 1-5 ]. Radio frequency (RF)-powered devices are the traditional choice for acceleration elements [ 5 5 - 8 8 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 9 - 17 17 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 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 19 , 20 20 ]. Theoretically, the DLA can provide acceleration gradients in excess of GV/m [ 21 21 - 24 24 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 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 26 , 27 27 ]. To make DLA competitive with larger devices, the Accelerator on a Chip International Program (ACHIP) [ 28 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 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 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 29 ]. In general, past research has focused on shape and size optimization, with examples shown in Fig. 1a and b 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 Fig. 1c )) of DLA structures was performed using adjoint methods by Hughes et al. [ 30 30 ]. When using the adjoint method for topology optimization, the adjoint source must be given [ 31 31 , 32 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. Theoretical descriptions 1 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 : Ω → R 3 such that the Maxwell equations <math display="block" id="M1"><mrow><mn>2</mn><mo>∇</mo><mo>⋅</mo><mi>B</mi><mo>=</mo><mn>0,</mn></mrow></math> <math display="block" id="M2"><mrow><mo>∇</mo><mo>⋅</mo><mi>D</mi><mo>=</mo><mi>ρ</mi><mn>,</mn></mrow></math> hold for a given 𝑱 : Ω → R 3 and ρ : Ω → R. We call 𝑬 the electric field, 𝑫 electric displacement field (or electric flux density), 𝑯 magnetic field intensity, 𝑩 magnetic flux density, and 𝑱 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 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. <math display="block" id="M3"><mrow><mi>E</mi><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mi>E</mi><mstyle mathvariant="bold" mathsize="normal"><mo stretchy="false">(</mo></mstyle><mi>r</mi><mstyle mathvariant="bold" mathsize="normal"><mo stretchy="false">)</mo></mstyle><msup><mi>e</mi><mrow><mo>−</mo><mi>i</mi><mi>ω</mi><mi>t</mi></mrow></msup><mtext>  and  </mtext><mi>H</mi><mo stretchy="false">(</mo><mi>r</mi><mn>,</mn><mi>t</mi><mo stretchy="false">)</mo><mo>=</mo><mi>E</mi><mstyle mathvariant="bold" mathsize="normal"><mo stretchy="false">(</mo></mstyle><mi>r</mi><mstyle mathvariant="bold" mathsize="normal"><mo stretchy="false">)</mo></mstyle><msup><mi>e</mi><mrow><mo>−</mo><mi>i</mi><mi>ω</mi><mi>t</mi></mrow></msup></mrow></math> We assume that all fields have the e -iωt dependence, with ω being the angular frequency. Through the curl operation, Maxwell’s equations can be expressed as <math display="block" id="M4"><mrow><mfrac><mn>1</mn><mi>μ</mi></mfrac><mo>∇</mo><mo>×</mo><mo>∇</mo><mo>×</mo><mi>E</mi><mstyle mathvariant="bold" mathsize="normal"><mo stretchy="false">(</mo></mstyle><mi>r</mi><mstyle mathvariant="bold" mathsize="normal"><mo stretchy="false">)</mo></mstyle><mo>−</mo><msup><mi>ω</mi><mn>2</mn></msup><mi>ϵ</mi><mi>E</mi><mstyle mathvariant="bold" mathsize="normal"><mo stretchy="false">(</mo></mstyle><mi>r</mi><mstyle mathvariant="bold" mathsize="normal"><mo stretchy="false">)</mo></mstyle><mo>=</mo><mi>i</mi><mi>ω</mi><mi>J</mi><mstyle mathvariant="bold" mathsize="normal"><mo stretchy="false">(</mo></mstyle><mi>r</mi><mstyle mathvariant="bold" mathsize="normal"><mo stretchy="false">)</mo></mstyle><mn>,</mn></mrow></math> 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 μ r and permittivity ϵ = ϵ 0 ϵ r is considered, where μ 0 is the vacuum permeability, and ϵ 0 is the vacuum permittivity. The working principle of FDFD [ 34 34 , 35 35 ] is to represent the electric field, dielectric constant distribution, and source on the Yee grid [ 33 33 ], and then represent the <math id="M5"><mrow><mo>∇</mo><mo>×</mo><mo>∇</mo><mo>×</mo></mrow></math> operator as a finite-difference differential matrix <math id="M6"><mrow><mover accent="true"><mtext>A</mtext><mo stretchy="true">˜</mo></mover></mrow></math> . In the discrete case, Eq. (1) can be written in matrix form as <math display="block" id="M7"><mrow><mo stretchy="false">[</mo><mover accent="true"><mi>A</mi><mo stretchy="true">˜</mo></mover><mo>−</mo><msup><mi>ω</mi><mn>2</mn></msup><msub><mi>ϵ</mi><mn>0</mn></msub><mover accent="true"><mrow><msub><mi>ϵ</mi><mtext>r</mtext></msub></mrow><mo stretchy="true">˜</mo></mover><mo stretchy="false">]</mo><mi>E</mi><mstyle mathvariant="bold" mathsize="normal"><mo stretchy="false">(</mo></mstyle><mi>r</mi><mstyle mathvariant="bold" mathsize="normal"><mo stretchy="false">)</mo></mstyle><mo>=</mo><mover accent="true"><mi>S</mi><mo stretchy="true">˜</mo></mover><mn>,</mn></mrow></math> where <math id="M8"><mrow><mover accent="true"><mi>A</mi><mo stretchy="true">˜</mo></mover></mrow></math> is a matrix that performs the finite difference derivative <math id="M9"><mrow><mfrac><mn>1</mn><mi>μ</mi></mfrac><mo>∇</mo><mo>×</mo><mo>∇</mo><mo>×</mo></mrow></math> , and does not depend on our permittivity distribution. <math id="M10"><mrow><mover accent="true"><mrow><msub><mi>ϵ</mi><mi>r</mi></msub></mrow><mo stretchy="true">˜</mo></mover></mrow></math> is a diagonal complex-valued matrix and each element along the diagonal corresponds to the relative permittivity of a cell in our domain. <math id="M11"><mrow><mover accent="true"><mi>S</mi><mo stretchy="true">˜</mo></mover></mrow></math> is a complex vector that represents the distribution of the driving current source. In practice, when applying the FDFD method, we first construct the matrix <math id="M12"><mrow><mover accent="true"><mi>ϕ</mi><mo stretchy="true">˜</mo></mover><mo>=</mo><mfrac><mn>1</mn><mi>μ</mi></mfrac><mo>∇</mo><mo>×</mo><mo>∇</mo><mo>×</mo><mo>−</mo><msup><mi>ω</mi><mn>2</mn></msup><msub><mo>∈</mo><mn>0</mn></msub><msub><mrow><mover accent="true"><mo>∈</mo><mo stretchy="true">˜</mo></mover></mrow><mi>r</mi></msub></mrow></math> and the source vector <math id="M13"><mrow><mover accent="true"><mi>S</mi><mo stretchy="true">˜</mo></mover></mrow></math> and then solve for the electric field <math id="M14"><mrow><mi>E</mi><mo>=</mo><mover accent="true"><mrow><msup><mover accent="true"><mi>ϕ</mi><mo>˜</mo></mover><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow><mo stretchy="true">˜</mo></mover><mover accent="true"><mi>S</mi><mo stretchy="true">˜</mo></mover></mrow></math> using an iterative or direct method. Because the fundamental features of DLA appear in the simplified 2D analysis, the full-field variation with z is not considered for the time being. This means that <math display="block" id="M15"><mrow><mo>∂</mo><mi>E</mi><mo>/</mo><mo>∂</mo><mi>z</mi><mo>=</mo><mn>0</mn><mtext>  and  </mtext><mo>∂</mo><mi>H</mi><mo>/</mo><mo>∂</mo><mi>z</mi><mo>=</mo><mn>0</mn></mrow></math> 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 H z and the TE mode by E z , as shown in Fig. 2 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 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 Fig. 3a . To obtain the device structure under the optimal objective function ξ ( x i ) 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 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 be run 2 10 times. 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 ξ ( x i ). To make the computation of the derivative of the objective function ξ ( x i ) 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. The objective function ξ( x i ) 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 Fig. 4 . The purpose is to find the derivative <math id="M16"><mrow><mo>∇</mo><mi>ξ</mi><mo stretchy="false">(</mo><msub><mi>x</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>=</mo><mo stretchy="false">[</mo><mo>∂</mo><mi>ξ</mi><mo>/</mo><mo>∂</mo><msub><mi>x</mi><mn>1</mn></msub><mtext> </mtext><mn>,</mn><mtext> </mtext><mo>∂</mo><mi>ξ</mi><mo>/</mo><mo>∂</mo><msub><mi>x</mi><mn>2</mn></msub><mo stretchy="false">]</mo></mrow></math> , so we apply the chain rule to the first term, <math id="M17"><mrow><mo>∂</mo><mi>ξ</mi><mo>/</mo><mo>∂</mo><msub><mi>x</mi><mn>1</mn></msub></mrow></math> , which yields <math display="block" id="M18"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mrow><mo>∂</mo><mi>ξ</mi></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>1</mn></msub></mrow></mfrac><mo>=</mo><mfrac><mrow><mo>∂</mo><mi>ξ</mi></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>6</mn></msub></mrow></mfrac><mo stretchy="false">(</mo><mfrac><mrow><mo>∂</mo><msub><mi>x</mi><mn>6</mn></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>5</mn></msub></mrow></mfrac><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><msub><mi>x</mi><mn>5</mn></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>3</mn></msub></mrow></mfrac><mfrac><mrow><mo>∂</mo><msub><mi>x</mi><mn>3</mn></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>1</mn></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>+</mo><mfrac><mrow><mo>∂</mo><msub><mi>x</mi><mn>6</mn></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>4</mn></msub></mrow></mfrac><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><msub><mi>x</mi><mn>4</mn></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>3</mn></msub></mrow></mfrac><mfrac><mrow><mo>∂</mo><msub><mi>x</mi><mn>3</mn></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>1</mn></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo stretchy="false">)</mo><mn>,</mn></mrow></mtd></mtr></mtable></mrow></math> or <math display="block" id="M19"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mfrac><mrow><mo>∂</mo><mi>ξ</mi></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>1</mn></msub></mrow></mfrac><mo>=</mo><mo stretchy="false">(</mo><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><mi>ξ</mi></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>6</mn></msub></mrow></mfrac><mfrac><mrow><mo>∂</mo><msub><mi>x</mi><mn>6</mn></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>5</mn></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mfrac><mrow><mo>∂</mo><msub><mi>x</mi><mn>5</mn></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>3</mn></msub></mrow></mfrac><mo>+</mo><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><mi>ξ</mi></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>6</mn></msub></mrow></mfrac><mfrac><mrow><mo>∂</mo><msub><mi>x</mi><mn>6</mn></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>4</mn></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mfrac><mrow><mo>∂</mo><msub><mi>x</mi><mn>4</mn></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>3</mn></msub></mrow></mfrac><mo stretchy="false">)</mo><mfrac><mrow><mo>∂</mo><msub><mi>x</mi><mn>3</mn></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mn>1</mn></msub></mrow></mfrac><mn>.</mn></mrow></mtd></mtr></mtable></mrow></math> We require only one forward sweep to obtain the result. The calculation graph from left to right corresponds to Eq. (3) from inside to outside. If we use the chain rule in another way, such as using Eq. (4), 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 d x i /d x j for any i, j using the rules of differentiation and some knowledge of the partial derivatives of each operation (see [ 37 37 , 38 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 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: <math display="block" id="M20"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mi>f</mi><mo stretchy="false">(</mo><msub><mi>μ</mi><mtext>r</mtext></msub><mo stretchy="false">)</mo><mo>=</mo><mi>Θ</mi><mi>H</mi><mo stretchy="false">(</mo><msub><mi>μ</mi><mn>1</mn></msub><mn>,</mn><msub><mi>μ</mi><mn>2</mn></msub><mo stretchy="false">)</mo><mn>,</mn></mrow></mtd></mtr></mtable></mrow></math> <math display="block" id="M21"><mrow><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mi>H</mi><mo stretchy="false">(</mo><msub><mi>μ</mi><mn>1</mn></msub><mn>,</mn><msub><mi>μ</mi><mn>2</mn></msub><mo stretchy="false">)</mo><mo>=</mo><msub><mi>μ</mi><mn>1</mn></msub><msup><mi>Θ</mi><mrow><mo>−</mo><mn>1</mn></mrow></msup><mo>+</mo><mrow><mo>(</mo><mrow><msub><mi>μ</mi><mn>2</mn></msub><mo>−</mo><msub><mi>μ</mi><mn>1</mn></msub></mrow><mo>)</mo></mrow><mn>,</mn><mi>Θ</mi><mo>=</mo><msub><mi>Θ</mi><mn>1</mn></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mn>1</mn><mo>+</mo><msup><mi>e</mi><mrow><mo>−</mo><mi>a</mi><mi>x</mi></mrow></msup></mrow></mfrac><mn>,</mn></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><mi>H</mi><mo stretchy="false">(</mo><msub><mi>μ</mi><mn>1</mn></msub><mn>,</mn><msub><mi>μ</mi><mn>2</mn></msub><mo stretchy="false">)</mo><mo>=</mo><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>μ</mi><mn>1</mn></msub><mo>+</mo><msub><mi>μ</mi><mn>2</mn></msub></mrow><mn>2</mn></mfrac></mrow><mo>)</mo></mrow><msup><mi>Θ</mi><mrow><mo>−</mo><mn>1</mn></mrow></msup><mo>+</mo><mrow><mo>(</mo><mrow><mfrac><mrow><msub><mi>μ</mi><mn>2</mn></msub><mo>−</mo><msub><mi>μ</mi><mn>1</mn></msub></mrow><mn>2</mn></mfrac></mrow><mo>)</mo></mrow><mn>,</mn><mi>Θ</mi><mo>=</mo><msub><mi>Θ</mi><mn>2</mn></msub><mo>=</mo><mfrac><mrow><mn>1</mn><mo>−</mo><msup><mi>e</mi><mrow><mo>−</mo><mi>a</mi><mi>x</mi></mrow></msup></mrow><mrow><mn>1</mn><mo>+</mo><msup><mi>e</mi><mrow><mo>−</mo><mi>a</mi><mi>x</mi></mrow></msup></mrow></mfrac><mn>.</mn></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> where Θ is a common s-shaped function, as shown in Fig. 4 Fig. 4 . We demonstrate two specific modification methods involving Θ 1 and Θ 2 , respectively, where Θ 1 is a “Sigmoid-like function”, and Θ 2 is a “Tanh-like function”, and a represents the hyperparameters. S Θ 1 and S Θ 2 are 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. <math display="block" id="M22"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mi>κ</mi><mo>=</mo><mfrac><mrow><mi>tanh</mi><mo stretchy="false">(</mo><mi>τ</mi><mi>σ</mi><mo stretchy="false">)</mo><mo>+</mo><mi>tanh</mi><mo stretchy="false">(</mo><mi>τ</mi><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><mo>−</mo><mi>σ</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo></mrow><mrow><mi>tanh</mi><mo stretchy="false">(</mo><mi>τ</mi><mi>σ</mi><mo stretchy="false">)</mo><mo>+</mo><mi>tanh</mi><mo stretchy="false">(</mo><mi>τ</mi><mo stretchy="false">(</mo><mn>1</mn><mo>−</mo><mi>σ</mi><mo stretchy="false">)</mo><mo stretchy="false">)</mo></mrow></mfrac></mrow></mtd></mtr></mtable></mrow></math> 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 32 ]. Results and discussion 1 The two important parameters ( ϵ and μ in Eq. (1)) that affect the electromagnetic distribution can be written in the following form: <math display="block" id="M23"><mrow><msub><mi>ϵ</mi><mtext>r</mtext></msub><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>ε</mi><mrow><mn>11</mn></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>ε</mi><mrow><mn>12</mn></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>ε</mi><mrow><mn>13</mn></mrow></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>ε</mi><mrow><mn>21</mn></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>ε</mi><mrow><mn>22</mn></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>ε</mi><mrow><mn>23</mn></mrow></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>ε</mi><mrow><mn>31</mn></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>ε</mi><mrow><mn>32</mn></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>ε</mi><mrow><mn>33</mn></mrow></msub></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mn>,</mn><msub><mi>μ</mi><mtext>r</mtext></msub><mo>=</mo><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>μ</mi><mrow><mn>11</mn></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>μ</mi><mrow><mn>12</mn></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>μ</mi><mrow><mn>13</mn></mrow></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>μ</mi><mrow><mn>21</mn></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>μ</mi><mrow><mn>22</mn></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>μ</mi><mrow><mn>23</mn></mrow></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>μ</mi><mrow><mn>31</mn></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>μ</mi><mrow><mn>32</mn></mrow></msub></mrow></mtd><mtd><mrow><msub><mi>μ</mi><mrow><mn>33</mn></mrow></msub></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mn>.</mn></mrow></math> Now, consider a material that consists of infinite columns in the x 3 -direction modeled by a permittivity such that ( x 1 , x 2 ) ↦ ϵ( x 1 , x 2 ). We assume an isotropic linear dielectric medium with permeability μ = μ 0 μ r and 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 <math display="block" id="M24"><mrow><msub><mi>ε</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>=</mo><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>ϵ</mi><mi>r</mi></msub><mn>,</mn></mrow></mtd><mtd columnalign="left"><mrow><mi>i</mi><mo>=</mo><mi>j</mi></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><mn>0,</mn></mrow></mtd><mtd columnalign="left"><mrow><mi>i</mi><mo>≠</mo><mi>j</mi></mrow></mtd></mtr></mtable></mrow></mrow><mtext>   and   </mtext><msub><mi>μ</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>=</mo><mrow><mo>{</mo><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mn>1</mn><mo>,</mo></mrow></mtd><mtd columnalign="left"><mrow><mi>i</mi><mo>=</mo><mi>j</mi></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><mn>0</mn><mo>,</mo></mrow></mtd><mtd columnalign="left"><mrow><mi>i</mi><mo>≠</mo><mi>j</mi><mn>.</mn></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> 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). <math display="block" id="M25"><mtable columnalign="left"><mtr><mtd><mrow><mo>(</mo><mrow><msubsup><mo>∇</mo><mi>x</mi><mtext>E</mtext></msubsup><msubsup><mi>ε</mi><mrow><mn>22</mn></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msubsup><msubsup><mo>∇</mo><mi>x</mi><mtext>H</mtext></msubsup><mo>+</mo><msubsup><mo>∇</mo><mi>y</mi><mtext>E</mtext></msubsup><msubsup><mi>ε</mi><mrow><mn>11</mn></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msubsup><msubsup><mo>∇</mo><mi>y</mi><mtext>H</mtext></msubsup><mo>+</mo><msup><mi>ω</mi><mn>2</mn></msup><mo>/</mo><msup><mi>c</mi><mn>2</mn></msup></mrow><mo>)</mo></mrow><msub><mi>H</mi><mi>z</mi></msub><mo>=</mo><mi>i</mi><mi>ω</mi><mi>J</mi><mo stretchy="false">(</mo><mi>z</mi><mo stretchy="false">)</mo></mtd></mtr><mtr><mtd><msubsup><mo>∇</mo><mi>y</mi><mtext>H</mtext></msubsup><msub><mi>H</mi><mi>z</mi></msub><mo>=</mo><mi>i</mi><mi>ω</mi><msub><mi>ϵ</mi><mn>0</mn></msub><msub><mi>ε</mi><mrow><mn>11</mn></mrow></msub><msub><mi>E</mi><mi>x</mi></msub><mn>,</mn><msubsup><mo>∇</mo><mi>x</mi><mtext>H</mtext></msubsup><msub><mi>H</mi><mi>z</mi></msub><mo>=</mo><mi>i</mi><mi>ω</mi><msub><mi>ϵ</mi><mn>0</mn></msub><msub><mi>ε</mi><mrow><mn>22</mn></mrow></msub><msub><mi>E</mi><mi>y</mi></msub></mtd></mtr></mtable></math> The terms <math id="M26"><mrow><msubsup><mo>∇</mo><mi>x</mi><mtext>E</mtext></msubsup></mrow></math> and <math id="M27"><mrow><msubsup><mo>∇</mo><mi>y</mi><mtext>E</mtext></msubsup></mrow></math> are 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 <math id="M28"><mrow><msubsup><mo>∇</mo><mi>x</mi><mtext>E</mtext></msubsup></mrow></math> and <math id="M29"><mrow><msubsup><mo>∇</mo><mi>y</mi><mtext>E</mtext></msubsup></mrow></math> for this simple case are <math display="block" id="M30"><mrow><mn>2</mn><msubsup><mo>∇</mo><mi>x</mi><mtext>E</mtext></msubsup><mo>=</mo><mfrac><mn>1</mn><mrow><mi>Δ</mi><mi>y</mi></mrow></mfrac><mrow><mo>[</mo><mrow><mtable><mtr><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mn>1</mn></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mn>1</mn></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow><msub><mi>U</mi><mi>x</mi></msub></mrow></mtd><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mn>1</mn></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mn>1</mn></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><msub><mi>U</mi><mi>x</mi></msub></mrow></mtd><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mn>1</mn></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mn>1</mn></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><msub><mi>U</mi><mi>x</mi></msub></mrow></mtd><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd></mtr></mtable></mrow><mo>]</mo></mrow></mrow></math> <math display="block" id="M31"><mrow><mn>2</mn><msubsup><mo>∇</mo><mi>y</mi><mtext>E</mtext></msubsup><mo>=</mo><mfrac><mn>1</mn><mrow><mi>Δ</mi><mi>x</mi></mrow></mfrac><mrow><mo>[</mo><mrow><mtable><mtr><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mn>1</mn></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mn>1</mn></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mn>1</mn></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mn>1</mn></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mn>1</mn></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mn>1</mn></mtd></mtr><mtr><mtd><mrow><msub><mi>V</mi><mi>y</mi></msub></mrow></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow><msub><mi>V</mi><mi>y</mi></msub></mrow></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd><mtd><mrow/></mtd></mtr><mtr><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><msub><mi>V</mi><mi>y</mi></msub></mrow></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow/></mtd><mtd><mrow><mo>−</mo><mn>1</mn></mrow></mtd></mtr></mtable></mrow><mo>]</mo></mrow></mrow></math> where <math id="M32"><mrow><mi>Δ</mi><mi>x</mi><mo stretchy="false">(</mo><mi>y</mi><mo stretchy="false">)</mo></mrow></math> denotes the FDFD Yee grid size along the x ( y ) direction, <math id="M33"><mrow><msub><mi>U</mi><mi>x</mi></msub><mo>=</mo><mi>exp</mi><mo stretchy="false">(</mo><mi>i</mi><msubsup><mi>k</mi><mi>x</mi><mrow><mtext>bloch</mtext></mrow></msubsup><mo stretchy="false">)</mo></mrow></math> , <math id="M34"><mrow><msub><mi>V</mi><mi>y</mi></msub><mo>=</mo><mi>exp</mi><mo stretchy="false">(</mo><mi>i</mi><msubsup><mi>k</mi><mi>y</mi><mrow><mtext>bloch</mtext></mrow></msubsup><mo stretchy="false">)</mo></mrow></math> , and <math id="M35"><mrow><msubsup><mi>k</mi><mrow><mi>x</mi><mo stretchy="false">(</mo><mi>y</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>bloch</mtext></mrow></msubsup></mrow></math> denotes the phase across the periodic boundary in the x ( y ) direction. Similarly, the terms <math id="M36"><mrow><msubsup><mo>∇</mo><mi>x</mi><mtext>H</mtext></msubsup></mrow></math> and <math id="M37"><mrow><msubsup><mo>∇</mo><mi>y</mi><mtext>H</mtext></msubsup></mrow></math> calculate first order spatial derivatives of the magnetic fields across the grid. The derivative operators for the electric and magnetic fields are different because of the staggered nature of the Yee grid ( Fig. 2 Fig. 2 ). However, they are related through <math display="block" id="M38"><mrow><mn>2</mn><msubsup><mo>∇</mo><mi>x</mi><mtext>H</mtext></msubsup><mo>=</mo><mo>−</mo><msup><mrow><mo stretchy="false">(</mo><msubsup><mo>∇</mo><mi>x</mi><mtext>E</mtext></msubsup><mo stretchy="false">)</mo></mrow><mrow><mtext>HT</mtext></mrow></msup><mn>,</mn><mtext> </mtext><mtext> </mtext><msubsup><mo>∇</mo><mi>y</mi><mtext>H</mtext></msubsup><mo>=</mo><mo>−</mo><msup><mrow><mo stretchy="false">(</mo><msubsup><mo>∇</mo><mi>y</mi><mtext>E</mtext></msubsup><mo stretchy="false">)</mo></mrow><mrow><mtext>HT</mtext></mrow></msup><mn>.</mn></mrow></math> The “HT” superscript indicates a Hermitian transpose operation. For the TE mode, averaging is needed for ϵ 11 and ϵ 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 <math id="M39"><mrow><msub><mi>Ξ</mi><mrow><mtext>gain</mtext></mrow></msub><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo></mrow></math> can be calculated using a periodic structure as follows. <math display="block" id="M40"><mrow><msub><mi>Ξ</mi><mrow><mtext>gain</mtext></mrow></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mi>y</mi></msub><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>1</mn><mi>λ</mi></mfrac><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mi>λ</mi><mo>/</mo><mn>2</mn></mrow><mrow><mo>−</mo><mi>λ</mi><mo>/</mo><mn>2</mn></mrow></msubsup></mrow></mstyle><msub><mi>E</mi><mi>y</mi></msub><mrow><mo>(</mo><mrow><mi>x</mi><mn>,</mn><mi>y</mi><mo>;</mo><mi>t</mi><mo>=</mo><mfrac><mi>y</mi><mrow><mi>β</mi><mi>c</mi></mrow></mfrac></mrow><mo>)</mo></mrow><mtext>d</mtext><mi>y</mi></mrow></math> Here, we focus only on the relationship between the permittivity distribution ( ϵ r ) and the energy gain ( <math id="M41"><mrow><msub><mi>Ξ</mi><mrow><mtext>gain</mtext></mrow></msub><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo></mrow></math> ). In other words, this is exactly the same case for nonmagnetic materials mentioned earlier. As described in Sect. 2, we must compute <math id="M42"><mrow><msub><mi>Ξ</mi><mrow><mtext>gain</mtext></mrow></msub><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo></mrow></math> for ϵ r as a derivative of <math display="block" id="M43"><mrow><mfrac><mrow><mo>∂</mo><msub><mi>Ξ</mi><mrow><mtext>gain</mtext></mrow></msub><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><msub><mi>ϵ</mi><mtext>r</mtext></msub></mrow></mfrac><mo>=</mo><mi>R</mi><mi>e</mi><mi>a</mi><mi>l</mi><mrow><mo>{</mo><mrow><mn>2</mn><mfrac><mrow><mo>∂</mo><msub><mi>Ξ</mi><mrow><mtext>gain</mtext></mrow></msub><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>E</mi></mrow></mfrac><mfrac><mrow><mo>∂</mo><mi>E</mi></mrow><mrow><mo>∂</mo><msub><mi>ϵ</mi><mtext>r</mtext></msub></mrow></mfrac></mrow><mo>}</mo></mrow><mn>,</mn></mrow></math> where <math id="M44"><mrow><msub><mi>Ξ</mi><mrow><mtext>gain</mtext></mrow></msub><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo></mrow></math> is equivalent to the objective function ( ξ ( x i )) mentioned previously. <math id="M45"><mrow><mo>∂</mo><msub><mi>Ξ</mi><mrow><mtext>gain</mtext></mrow></msub><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo><mo>/</mo><mo>∂</mo><msub><mi>ϵ</mi><mtext>r</mtext></msub></mrow></math> is an important parameter that indicates how to change ϵ r to increase <math id="M46"><mrow><msub><mi>Ξ</mi><mrow><mtext>gain</mtext></mrow></msub><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo></mrow></math> . Two optimization schemes for maximizing the acceleration gradient function <math id="M47"><mrow><msub><mi>Ξ</mi><mrow><mtext>gain</mtext></mrow></msub><mo stretchy="false">(</mo><msub><mi>E</mi><mi>y</mi></msub><mo stretchy="false">)</mo></mrow></math> 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 <math display="block" id="M48"><mrow><mtext>maximum</mtext><mtext> </mtext><mi>Γ</mi><mo>=</mo><mfrac><mrow><msub><mi>Ξ</mi><mrow><mtext>gain</mtext></mrow></msub><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mi>E</mi><mn>0</mn></msub></mrow></mfrac><mn>,</mn></mrow></math> where E 0 is the initial injection electric field, which should be less than the material damage threshold ( E material 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 <math display="block" id="M49"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><mtext>maximum</mtext></mrow></mtd><mtd columnalign="left"><mrow><mi>χ</mi><mo>=</mo><mfrac><mrow><msub><mi>Ξ</mi><mrow><mtext>gain</mtext></mrow></msub><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mi>E</mi><mrow><mtext>max</mtext></mrow></msub></mrow></mfrac><mo>,</mo></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow><mtext>subject to</mtext></mrow></mtd><mtd columnalign="left"><mrow><msub><mi>E</mi><mrow><mtext>max</mtext></mrow></msub><mo><</mo><msub><mi>E</mi><mrow><mtext>material</mtext><mtext> </mtext><mtext>damage</mtext></mrow></msub><mo>.</mo></mrow></mtd></mtr></mtable></mrow></math> Here, E max is the maximum value of the electric field magnitude vector in the optimization design region (ODR) and E material 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. The results of “ Case 1 ” are shown in Fig. 6 Fig. 6 and Fig. 7 Fig. 7 . The parameter settings in Fig 6 Fig 6 are the same as those in [ 30 30 ], in which the purpose is to test the solver and optimization scheme involved in this study. To accomplish this, a plane wave ( E 0 is the initial injected electric field) is introduced into the DLA chip from the left side using the periodic condition. The laser wavelength ( λ ) used in the simulations is 2000 nm. The normalized velocity ( β ) of the injected electrons is 0.5, and the square grid size (Δ x ( y )) used in the simulation is 20 nm. The material used in the optimization process is fused silica (SiO 2 ), the index is n = 1.45, that is, ε 11 = ε 22 = 2.1). The cross-sectional area utilized for the optimization design region is 5 β λ × λ , and the particle injection channel region is 5 β λ × <math id="M50"><mrow><mfrac><mn>1</mn><mn>5</mn></mfrac><mi>λ</mi></mrow></math> . An acceleration gradient function ( <math id="M51"><mrow><msub><mi>Ξ</mi><mrow><mtext>gain</mtext></mrow></msub><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo></mrow></math> ) of 0.32E 0 ( Γ = 0.32) can be achieved by optimizing the maximum conversion factor ( Γ ). The results of the gradient calculation and structure optimization are consistent with those in Ref. [ 30 30 ] (the achieved gradient is 0.31 E 0 ). As calculated based on [ 30 30 ], Fig. 6 Fig. 6 corresponds to an acceleration gradient of 1.37 GV/m. It can be observed from Fig. 6 Fig. 6 that the structural design given by the optimization results tends to include mirrors on both sides of the structure. By observing the normalized electric field diagram in the y direction ( E y / <math id="M52"><mrow><mrow><mo>|</mo><mrow><msub><mi>E</mi><mn>0</mn></msub></mrow><mo>|</mo></mrow></mrow></math> ), it can be observed that for higher acceleration gradients, it is useful to use dielectric mirror resonance to enhance the field in the gap. The parameters in Fig. 7 Fig. 7 are the same as those used in Ref. [ 29 29 ], and the purpose is to highlight the advantage of using the topology optimization method. A plane wave ( E 0 is the initial injected electric field) is introduced into the DLA chip from the left side using the periodic condition. The material used in the optimization process is crystal quartz, the index is n = 1.55, i.e., ε 11 = ε 22 = 2.4. The cross-sectional area utilized for the optimization design region is 4 β λ × λ , and the particle injection channel region is 4 β λ × <math id="M53"><mrow><mfrac><mn>1</mn><mn>4</mn></mfrac><mi>λ</mi></mrow></math> . The acceleration gradient function ( <math id="M54"><mrow><msub><mi>Ξ</mi><mrow><mtext>gain</mtext></mrow></msub><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo></mrow></math> ) of 3.1E 0 ( Γ = 3.1) is achieved by optimizing the maximum conversion factor ( Γ ). In Ref. [ 29 29 ], the gradient achieved was 0.49 E 0 ( Γ = 0.49). Therefore, the topology-optimized structures achieve an approximately six-fold improvement in the acceleration gradient. It is worth noting from Fig. 7a) Fig. 7a) that the optimized structure at this time deviates significantly from the intuitive physical structure. Although the structure differs greatly from the intuitive structure, it can be observed from the normalized electric field ( E y / <math id="M55"><mrow><mrow><mo>|</mo><mrow><msub><mi>E</mi><mn>0</mn></msub></mrow><mo>|</mo></mrow></mrow></math> ) in the y -direction, as shown in Fig. 7b Fig. 7b ), that the topologically optimized DLA structure has an obvious standing wave field distribution in the particle injection channel region (PICR), which benefits the acceleration of charged particles significantly. The results of “ Case 2 ” are shown in Fig. 8 Fig. 8 . The parameters used are the same as those used in Ref. [ 29 29 ]. The purpose is to point out that the method using automatic derivation can optimize the objective function more flexibly than the adjoint method, because in the “ Case 2 ” scheme, the max function ( E max ) is not differentiable, which is a case when it is impossible to directly provide the accurate adjoint source. The cross-sectional area utilized for the optimization design region is 4 β λ × λ , and the particle injection channel region is 4 β λ × <math id="M56"><mrow><mfrac><mn>1</mn><mn>4</mn></mfrac><mi>λ</mi></mrow></math> . An acceleration gradient function ( <math id="M57"><mrow><msub><mi>Ξ</mi><mrow><mtext>gain</mtext></mrow></msub><mo stretchy="false">(</mo><mi>E</mi><mo stretchy="false">)</mo></mrow></math> ) of 1.4 E 0 ( Γ = 1.4) is achieved by optimizing the maximum conversion factor ( Γ ). In Ref. [ 29 29 ], the gradient achieved is 0.49 E 0 ( Γ = 0.49). Therefore, the topology-optimized structures show an approximately three-fold improvement in the acceleration gradient compared to previously established structures. As in Fig. 8a Fig. 8a ), it can also be clearly observed that the optimized structure is very different from the intuitive physical structure. According to the normalized electric field ( E y / <math id="M58"><mrow><mrow><mo>|</mo><mrow><msub><mi>E</mi><mn>0</mn></msub></mrow><mo>|</mo></mrow></mrow></math> ), Fig. 8b Fig. 8b ) indicates that the electric field intensity of the topologically optimized DLA chip is weak. For comparison, we calculated the normalized maximum electric field ( E max / <math id="M59"><mrow><mrow><mo>|</mo><mrow><msub><mi>E</mi><mn>0</mn></msub></mrow><mo>|</mo></mrow></mrow></math> ) shown in Fig. 7a) Fig. 7a) and Fig. 8a) Fig. 8a) , and the results indicate that under the same parameters, the maximum electric field strength in the “ Case 2 ” optimized structure is only one-fifth of the maximum electric field strength in the “ Case 1 ” optimized structure. Therefore, it can be shown that when the “ Case 2 ” scheme is implemented, the electric field intensity in the structure is reduced, which means that the structure is less likely to suffer damage. 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 Fig. 7 is three times that in Fig. 8 Fig. 8 . If Γ is used as the evaluation index, the result for the acceleration gradient function in Fig. 8 Fig. 8 is twice that of Fig. 7 Fig. 7 . The acceleration gradients of the two schemes are different, but they are significantly improved when compared to the grating structure [ 29 29 ]. A phenomenological explanation is that in the y-direction of the topologically optimized DLA structures ( Fig. 7b Fig. 7b ) and Fig. 8b Fig. 8b )), the electric field ( E y / <math id="M60"><mrow><mrow><mo>|</mo><mrow><msub><mi>E</mi><mn>0</mn></msub></mrow><mo>|</mo></mrow></mrow></math> ) 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 29 ]. Summary 1 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 E x that the particle feels. However, the effect of E x will 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 E x in 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. References B.W.J. McNeil , N.R. Thompson , X-ray free-electron lasers . Nature Photon 4 , 814 - 821 ( 2010 ). doi: 10.1038/nphoton.2010.239 http://doi.org/10.1038/nphoton.2010.239 X-ray free-electron lasers Y. Zhang , W.C. Fang , X.X. Huang et al. , Design, fabrication, and cold test of an S-band high-gradient accelerating structure for compact proton therapy facility . Nucl. Sci. Tech. 32 , 38 ( 2021 ). doi: 10.1007/s41365-021-00869-z http://doi.org/10.1007/s41365-021-00869-z Design, fabrication, and cold test of an S-band high-gradient accelerating structure for compact proton therapy facility S. V. Bulanov and V. S. Khoroshkov , Feasibility of using laser ion accelerators in proton therapy . Plasma Phys. Rep. 28 , 453 - 456 ( 2002 ). doi: 10.1134/1.1478534 http://doi.org/10.1134/1.1478534 Feasibility of using laser ion accelerators in proton therapy W. Nan , B. Guo , C.J. Lin et al. , First proof-of-principle experiment with the post-accelerated isotope separator on-line beam at BRIF: measurement of the angular distribution of 23 Na + 40 Ca elastic scattering . Nucl. Sci. Tech. 32 , 53 ( 2021 ). doi: 10.1007/s41365-021-00889-9 http://doi.org/10.1007/s41365-021-00889-9 First proof-of-principle experiment with the post-accelerated isotope separator on-line beam at BRIF: measurement of the angular distribution of 23Na + 40Ca elastic scattering E. L. Ginzton , W. W. Hansen , and W. R. Kennedy , A Linear Electron Accelerator . Rev. Sci. Instrum. 19 , 89 - 108 ( 1948 ) doi: 10.1063/1.1741225 http://doi.org/10.1063/1.1741225 A Linear Electron Accelerator Z.-Y. Ma , S.-J. Zhao , X.-M. Liu et al. , High RF power tests of the first 1.3 GHz fundamental power coupler prototypes for the SHINE project . Nucl. Sci. Tech. 33 , 10 ( 2022 ). doi: 10.1007/s41365-022-00984-5 http://doi.org/10.1007/s41365-022-00984-5 High RF power tests of the first 1.3 GHz fundamental power coupler prototypes for the SHINE project C. Wang , J.H. Tan , X.X. Huang et al. , Design optimization and cold RF test of a 2.6-cell cryogenic RF gun . Nucl. Sci. Tech. 32 , 97 ( 2021 ). doi: 10.1007/s41365-021-00925-8 http://doi.org/10.1007/s41365-021-00925-8 Design optimization and cold RF test of a 2.6-cell cryogenic RF gun G.-R. Li , S.-X. Zheng , H.-J. Zeng et al. , Design and test of an RF acceleration system loaded with magnetic alloy for the proton synchrotron of the Xi’an Proton Application Facility . Nucl. Sci. Tech. 29 , 94 ( 2018 ). doi: 10.1007/s41365-018-0434-9 http://doi.org/10.1007/s41365-018-0434-9 Design and test of an RF acceleration system loaded with magnetic alloy for the proton synchrotron of the Xi’an Proton Application Facility R.J. England , R.J. Noble , K. Bane et al. , Dielectric laser accelerators . Rev. Mod. Phys. 86 , 1337 ( 2014 ). doi: 10.1103/RevModPhys.86.1337 http://doi.org/10.1103/RevModPhys.86.1337 Dielectric laser accelerators J. Breuer and P. Hommelhoff , Laser-based acceleration of nonrelativistic electrons at a dielectric structure . Phys. Rev. Lett. 111 , 134803 ( 2013 ). doi: 10.1103/PhysRevLett.111.134803 http://doi.org/10.1103/PhysRevLett.111.134803 Laser-based acceleration of nonrelativistic electrons at a dielectric structure E.A. Peralta , K. Soong , R.J. England et al. , Demonstration of electron acceleration in a laser-driven dielectric microstructure . Nature 503 , 91 - 94 ( 2013 ). doi: 10.1038/nature12664 http://doi.org/10.1038/nature12664 Demonstration of electron acceleration in a laser-driven dielectric microstructure K.J. Leedle , R.F. Pease , R.L. Byer et al. , Laser acceleration and deflection of 96.3 keV electrons with a silicon dielectric structure . Optica 2 , 158 - 161 ( 2015 ) doi: 10.1364/OPTICA.2.000158 http://doi.org/10.1364/OPTICA.2.000158 Laser acceleration and deflection of 96.3 keV electrons with a silicon dielectric structure K.J. Leedle , R.F. Pease , R.L. Byer et al. , Dielectric laser acceleration of sub-100 keV electrons with silicon dual-pillar grating structures . Optics Letters 40 , 4344 - 4347 ( 2015 ) doi: 10.1364/OL.40.004344 http://doi.org/10.1364/OL.40.004344 Dielectric laser acceleration of sub-100 keV electrons with silicon dual-pillar grating structures U. Niedermayer , T. Egenolf , O. Boine-Frankenheim et al. , Alternating phase focusing for dielectric laser acceleration . Phys. Rev. Lett. 121 , 214801 ( 2018 ) doi: 10.1103/PhysRevLett.121.214801 http://doi.org/10.1103/PhysRevLett.121.214801 Alternating phase focusing for dielectric laser acceleration D.S. Black , K.J. Leedle , Y. Miao et al. , Laser-driven electron lensing in silicon microstructures . Phys. Rev. Lett. 122 , 104801 ( 2019 ) doi: 10.1103/PhysRevLett.122.104801 http://doi.org/10.1103/PhysRevLett.122.104801 Laser-driven electron lensing in silicon microstructures N. Schönenberger , A. Mittelbach , P. Yousefi et al. , Generation and characterization of attosecond microbunched electron pulse trains via dielectric laser acceleration . Phys. Rev. Lett. 123 , 264803 ( 2019 ) doi: 10.1103/PhysRevLett.123.264803 http://doi.org/10.1103/PhysRevLett.123.264803 Generation and characterization of attosecond microbunched electron pulse trains via dielectric laser acceleration D.S. Black , U. Niedermayer , Y. Miao et al. , Net acceleration and direct measurement of attosecond electron pulses in a silicon dielectric laser accelerator . Phys. Rev. Lett. 123 , 264802 ( 2019 ) doi: 10.1103/PhysRevLett.123.264802 http://doi.org/10.1103/PhysRevLett.123.264802 Net acceleration and direct measurement of attosecond electron pulses in a silicon dielectric laser accelerator K.P. Wootton , J. McNeur and K. J. Leedle , Dielectric laser accelerators: Designs, experiments, and applications . Rev. Accelerator Sci. Technol. 09 , 105 - 126 ( 2016 ) doi: 10.1142/S179362681630005X http://doi.org/10.1142/S179362681630005X Dielectric laser accelerators: Designs, experiments, and applications K. Soong , R.L. Byer , E.R. Colby et al. , Laser damage threshold measurements of optical materials for direct laser accelerators . AIP Conference Proceedings 1507 , 511 ( 2012 ) doi: 10.1063/1.4773749 http://doi.org/10.1063/1.4773749 Laser damage threshold measurements of optical materials for direct laser accelerators R. Agustsson , E. Arab , A. Murokh et al. , Measuring single-shot, picosecond optical damage threshold in Ge, Si, and sapphire with a 5.1-µm laser . Optical Materials Express 5 , 2835 - 2842 ( 2015 ) doi: 10.1364/OME.5.002835 http://doi.org/10.1364/OME.5.002835 Measuring single-shot, picosecond optical damage threshold in Ge, Si, and sapphire with a 5.1-µm laser D. Cesar , J. Maxson , P. Musumeci et al. , Optical design for increased interaction length in a high gradient dielectric laser accelerator . Nucl. Instrum. Meth. A 909 , 252 - 256 ( 2018 ) doi: 10.1016/j.nima.2018.01.012 http://doi.org/10.1016/j.nima.2018.01.012 Optical design for increased interaction length in a high gradient dielectric laser accelerator J. Breuer , J. McNeur , P. Hommelhoff , Dielectric laser acceleration of electrons in the vicinity of single and double grating structures—theory and simulations . J. Phys. B At. Mol. Opt. Phys. 47 , 234004 ( 2014 ) doi: 10.1088/0953-4075/47/23/234004 http://doi.org/10.1088/0953-4075/47/23/234004 Dielectric laser acceleration of electrons in the vicinity of single and double grating structures—theory and simulations T. Plettner , R.L. Byer , B. Montazeri , Electromagnetic forces in the vacuum region of laser-driven layered grating structures . J. Mod. Opt. 58 , 1518 ( 2011 ) doi: 10.1080/09500340.2011.611914 http://doi.org/10.1080/09500340.2011.611914 Electromagnetic forces in the vacuum region of laser-driven layered grating structures D.S. Black , Z.X. Zhao , K.J. Leedle et al. , Operating modes of dual-grating dielectric laser accelerators . Phys. Rev. Accel. Beams 23 , 114001 ( 2020 ) doi: 10.1103/PhysRevAccelBeams.23.114001 http://doi.org/10.1103/PhysRevAccelBeams.23.114001 Operating modes of dual-grating dielectric laser accelerators D. Cesar , S. Custodio , J. Maxson et al. , High-field nonlinear optical response and phase control in a dielectric laser accelerator . Commun. Phys. 1 , 46 ( 2018 ). doi: 10.1038/s42005-018-0047-y http://doi.org/10.1038/s42005-018-0047-y High-field nonlinear optical response and phase control in a dielectric laser accelerator T.W. Hughes , S. Tan , Z.X. Zhao et al. , On-Chip Laser-Power Delivery System for Dielectric Laser Accelerators . Phys. Rev. Appl. 9 , 054017 ( 2018 ). doi: 10.1103/PhysRevApplied.9.054017 http://doi.org/10.1103/PhysRevApplied.9.054017 On-Chip Laser-Power Delivery System for Dielectric Laser Accelerators N.V Sapra , K.Y. Yang , D. Vercruysse et al. , On-chip integrated laser-driven particle accelerator . Science 367 , 79 - 83 ( 2020 ). doi: 10.1126/science.aay5734 http://doi.org/10.1126/science.aay5734 On-chip integrated laser-driven particle accelerator https://www.laserphysics.nat.fau.eu/research/achip/ https://www.laserphysics.nat.fau.eu/research/achip/ T. Plettner , P.P. Lu , R.L. Byer , Proposed few-optical cycle laser-driven particle accelerator structure . Phys. Rev. ST Accel. Beams 9 , 111301 ( 2006 ) doi: 10.1103/PhysRevSTAB.9.111301 http://doi.org/10.1103/PhysRevSTAB.9.111301 Proposed few-optical cycle laser-driven particle accelerator structure T. Hughes , G. Veronis , K.P. Wootton et al. , Method for computationally efficient design of dielectric laser accelerator structures . Opt. Express 25 , 15414 - 15427 ( 2017 ) doi: 10.1364/OE.25.015414 http://doi.org/10.1364/OE.25.015414 Method for computationally efficient design of dielectric laser accelerator structures C.M Lalau-Keraly , S. Bhargava , O.D Miller et al. , Adjoint shape optimization applied to electromagnetic design . Opt. Express 21 , 21693 - 21701 ( 2013 ). doi: 10.1364/OE.21.021693 http://doi.org/10.1364/OE.21.021693 Adjoint shape optimization applied to electromagnetic design A. Michaels , E. Yablonovitch , Leveraging continuous material averaging for inverse electromagnetic design . Opt. Express 26 , 31717 - 31737 ( 2018 ). doi: 10.1364/OE.26.031717 http://doi.org/10.1364/OE.26.031717 Leveraging continuous material averaging for inverse electromagnetic design K. Yee , Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media . IEEE T. Antennas Propagation 14 , 302 - 307 ( 1966 ). doi: 10.1109/TAP.1966.1138693 http://doi.org/10.1109/TAP.1966.1138693 Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media R.C. Rumpf , Simple implementation of arbitrarily shaped total-field/scattered-field regions in finite-difference frequency-domain . Prog. Electromagnetics Research B , 36 , 221 - 248 , ( 2012 ) doi: 10.2528/PIERB11092006 http://doi.org/10.2528/PIERB11092006 Simple implementation of arbitrarily shaped total-field/scattered-field regions in finite-difference frequency-domain T.W. Hughes , I.A.D. Williamson , M. Minkov et al. , Forward-mode differentiation of Maxwell’s equations . ACS Publications 6 , 3010 - 3016 ( 2019 ). doi: 10.1021/acsphotonics.9b01238 http://doi.org/10.1021/acsphotonics.9b01238 Forward-mode differentiation of Maxwell’s equations S.G. Johnson , Notes on perfectly matched layers (PMLs) . arxiv: 2108.05348 Notes on perfectly matched layers (PMLs) M. Minkov , I.A.D. Williamson , L.C. Andreani et al. , Inverse design of photonic crystals through automatic differentiation . ACS Photonics 7 , 1729 - 1741 ( 2020 ). doi: 10.1021/acsphotonics.0c00327 http://doi.org/10.1021/acsphotonics.0c00327 Inverse design of photonic crystals through automatic differentiation Autograd: Efficiently computes derivatives of numpy code . https://github.com/HIPS/autograd https://github.com/HIPS/autograd Autograd: Efficiently computes derivatives of numpy code D.P. Kingma , J. Ba , Adam: A method for stochastic optimization . arXiv:1412.6980 Adam: A method for stochastic optimization All 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. The 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. 10.1007/s41365-022-01101-2 This work was supported by the National Natural Science Foundation of China (Nos. 12004353, 11975214, 11991071, 11905202, 12174350), Key Laboratory Foundation of The Sciences and Technology on Plasma Physics Laboratory (No.6142A04200103), and Independent scientific research (No. JCKYS2021212011). 2022-04-08 2022-07-26 2022-07-31 2022-09-22 2022-10-24 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2022 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-022-01101-2 http://dx.doi.org/10.1007/s41365-022-01101-2 . Springer Nature or its licensor holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law. 2022 Nuclear Science and Techniques 09 33