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.

Fig. 1Full-screen View

(Color online) Three means of optimization: a) size optimization, which takes the structural shape or the shape of the hole as the optimization object to adjust the size; b) shape optimization, looking for new concave-convex distributions on the existing structure; c) topology optimization, which can obtain the optimal distribution in the material design space.

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

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³ such that the Maxwell equations $2\nabla \cdot B=\mathrm{0,}$ $\nabla \cdot D=\rho \mathrm{,}$ hold for a given 𝑱: Ω→ R³ 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], 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. $E(r\mathrm{,}t)=E(r)e^{-i\omega t}\text{ and }H(r\mathrm{,}t)=E(r)e^{-i\omega t}$ 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 $\frac{1}{\mu }\nabla \times \nabla \times E(r)-{\omega }^2ϵE(r)=i\omega J(r)\mathrm{,}$ (1) 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 μ = μ₀μr and permittivity ϵ = ϵ0ϵr is considered, where μ₀ is the vacuum permeability, and ϵ₀ is 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 $\nabla \times \nabla \times$ operator as a finite-difference differential matrix $\tilde{\text{A}}$. In the discrete case, Eq. (1) can be written in matrix form as $[\tilde{A}-{\omega }^2{ϵ}_0\widetilde{{ϵ}_{\text{r}}}]E(r)=\tilde{S}\mathrm{,}$ (2) where $\tilde{A}$ is a matrix that performs the finite difference derivative $\frac{1}{\mu }\nabla \times \nabla \times$, and does not depend on our permittivity distribution. $\widetilde{{ϵ}_r}$ is a diagonal complex-valued matrix and each element along the diagonal corresponds to the relative permittivity of a cell in our domain. $\tilde{S}$ 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 $\tilde{\varphi }=\frac{1}{\mu }\nabla \times \nabla \times -{\omega }^2{\in }_0{\tilde{\in }}_r$ and the source vector $\tilde{S}$ and then solve for the electric field $E=\widetilde{{\tilde{\varphi }}^{-1}}\tilde{S}$ 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 $\partial E/\partial z=0\text{ and }\partial H/\partial z=0$ (3) 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 Hz and 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.

Fig. 2Full-screen View

(Color online) Fundamental features of DLA used in the simplified 2D analysis, in which the full-field variation with 𝒛 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 E_z, and the TE mode can be fully described by H_z. The field distribution of the TM (TE) mode is discretized onto a square Yee grid, as shown in a and b

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.

Fig. 3Full-screen View

(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.

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 ¹⁰ 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 ξ(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.

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 $\nabla \xi (x_i)=[\partial \xi /\partial x_1\mathrm{,}\partial \xi /\partial x_2]$, so we apply the chain rule to the first term, $\partial \xi /\partial x_1$, which yields $\begin{array}{l}\frac{\partial \xi }{\partial x_1}=\frac{\partial \xi }{\partial x_6}(\frac{\partial x_6}{\partial x_5}\left(\frac{\partial x_5}{\partial x_3}\frac{\partial x_3}{\partial x_1}\right)+\frac{\partial x_6}{\partial x_4}\left(\frac{\partial x_4}{\partial x_3}\frac{\partial x_3}{\partial x_1}\right))\mathrm{,}\hfill \end{array}$ (4) or $\begin{array}{l}\frac{\partial \xi }{\partial x_1}=(\left(\frac{\partial \xi }{\partial x_6}\frac{\partial x_6}{\partial x_5}\right)\frac{\partial x_5}{\partial x_3}+\left(\frac{\partial \xi }{\partial x_6}\frac{\partial x_6}{\partial x_4}\right)\frac{\partial x_4}{\partial x_3})\frac{\partial x_3}{\partial x_1}\mathrm{.}\hfill \end{array}$ (5) 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 dxi/dxj for 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.

Fig. 4Full-screen View

(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.

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: $\begin{array}{l}f({\mu }_{\text{r}})=\Theta H({\mu }_1\mathrm{,}{\mu }_2)\mathrm{,}\hfill \end{array}$ (6) $\{\begin{array}{l}H({\mu }_1\mathrm{,}{\mu }_2)={\mu }_1{\Theta }^{-1}+\left({\mu }_2-{\mu }_1\right)\mathrm{,}\Theta ={\Theta }_1=\frac{1}{1+e^{-ax}}\mathrm{,}\hfill \\ H({\mu }_1\mathrm{,}{\mu }_2)=\left(\frac{{\mu }_1+{\mu }_2}{2}\right){\Theta }^{-1}+\left(\frac{{\mu }_2-{\mu }_1}{2}\right)\mathrm{,}\Theta ={\Theta }_2=\frac{1-e^{-ax}}{1+e^{-ax}}\mathrm{.}\hfill \end{array}$ (7) where Θ is a common s-shaped function, as shown in Fig. 4. We demonstrate two specific modification methods involving Θ₁ and Θ₂, respectively, where Θ₁ is a “Sigmoid-like function”, and Θ₂ is a “Tanh-like function”, and a represents the hyperparameters. SΘ₁ and SΘ₂ 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. $\begin{array}{l}\kappa =\frac{\tanh (\tau \sigma )+\tanh (\tau (\mu (r)-\sigma ))}{\tanh (\tau \sigma )+\tanh (\tau (1-\sigma ))}\hfill \end{array}$ (8) 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].

Fig. 5Full-screen View

(Color online) We show two specific modification methods using Θ₁ and Θ₂, where Θ₁ is a “Sigmoid function,” and Θ₂ is a “Tanh function.” SΘ₁ and SΘ₂ 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.

Results and discussion

The two important parameters (ϵ and μ in Eq. (1)) that affect the electromagnetic distribution can be written in the following form: ${ϵ}_{\text{r}}=\left(\begin{array}{ccc}{\epsilon }_{11}& {\epsilon }_{12}& {\epsilon }_{13}\\ {\epsilon }_{21}& {\epsilon }_{22}& {\epsilon }_{23}\\ {\epsilon }_{31}& {\epsilon }_{32}& {\epsilon }_{33}\end{array}\right)\mathrm{,}{\mu }_{\text{r}}=\left(\begin{array}{ccc}{\mu }_{11}& {\mu }_{12}& {\mu }_{13}\\ {\mu }_{21}& {\mu }_{22}& {\mu }_{23}\\ {\mu }_{31}& {\mu }_{32}& {\mu }_{33}\end{array}\right)\mathrm{.}$ (9) Now, consider a material that consists of infinite columns in the x₃-direction modeled by a permittivity such that (x₁, x₂) ↦ ϵ(x₁, x₂). We assume an isotropic linear dielectric medium with permeability μ = μ₀ μr and permittivity ϵ = ϵ₀ϵ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 ${\epsilon }_{ij}=\{\begin{array}{ll}{ϵ}_r\mathrm{,}\hfill & i=j\hfill \\ \mathrm{0,}\hfill & i\ne j\hfill \end{array}\text{ and }{\mu }_{ij}=\{\begin{array}{ll}1,\hfill & i=j\hfill \\ 0,\hfill & i\ne j\mathrm{.}\hfill \end{array}$ (10) 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). $\begin{array}{l}\left({\nabla }_x^{\text{E}}{\epsilon }_{22}^{-1}{\nabla }_x^{\text{H}}+{\nabla }_y^{\text{E}}{\epsilon }_{11}^{-1}{\nabla }_y^{\text{H}}+{\omega }^2/c^2\right)H_z=i\omega J(z)\\ {\nabla }_y^{\text{H}}{H}_z=i\omega {ϵ}_0{\epsilon }_{11}{E}_x\mathrm{,}{\nabla }_x^{\text{H}}{H}_z=i\omega {ϵ}_0{\epsilon }_{22}{E}_y\end{array}$ (11) The terms ${\nabla }_x^{\text{E}}$ and ${\nabla }_y^{\text{E}}$ 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 ${\nabla }_x^{\text{E}}$ and ${\nabla }_y^{\text{E}}$ for this simple case are $2{\nabla }_x^{\text{E}}=\frac{1}{\Delta y}\left[\begin{array}{ccccccccc}-1& 1& & & & & & & \\ & -1& 1& & & & & & \\ U_x& & -1& 0& & & & & \\ & & & -1& 1& & & & \\ & & & & -1& 1& & & \\ & & & U_x& & -1& 0& & \\ & & & & & & -1& 1& \\ & & & & & & & -1& 1\\ & & & & & & U_x& & -1\end{array}\right]$ $2{\nabla }_y^{\text{E}}=\frac{1}{\Delta x}\left[\begin{array}{ccccccccc}-1& & & 1& & & & & \\ & -1& & & 1& & & & \\ & & -1& & & 1& & & \\ & & & -1& & & 1& & \\ & & & & -1& & & 1& \\ & & & & & -1& & & 1\\ V_y& & & & & & -1& & \\ & V_y& & & & & & -1& \\ & & V_y& & & & & & -1\end{array}\right]$ where $\Delta x(y)$ denotes the FDFD Yee grid size along the x(y) direction, $U_x=\exp (i{k}_x^{\text{bloch}})$, $V_y=\exp (i{k}_y^{\text{bloch}})$, and $k_{x(y)}^{\text{bloch}}$ denotes the phase across the periodic boundary in the x(y) direction. Similarly, the terms ${\nabla }_x^{\text{H}}$ and ${\nabla }_y^{\text{H}}$ 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). However, they are related through $2{\nabla }_x^{\text{H}}=-{({\nabla }_x^{\text{E}})}^{\text{HT}}\mathrm{,}{\nabla }_y^{\text{H}}=-{({\nabla }_y^{\text{E}})}^{\text{HT}}\mathrm{.}$ The “HT” superscript indicates a Hermitian transpose operation. For the TE mode, averaging is needed for ϵ₁₁ and ϵ₂₂. 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 ${\Xi }_{\text{gain}}(E)$ can be calculated using a periodic structure as follows. ${\Xi }_{\text{gain}}(E_y)=\frac{1}{\lambda }{\displaystyle {\int }_{\lambda /2}^{-\lambda /2}}{E}_y\left(x\mathrm{,}y;t=\frac{y}{\beta c}\right)\text{d}y$ Here, we focus only on the relationship between the permittivity distribution (ϵr) and the energy gain (${\Xi }_{\text{gain}}(E)$). In other words, this is exactly the same case for nonmagnetic materials mentioned earlier. As described in Sect. 2, we must compute ${\Xi }_{\text{gain}}(E)$ for ϵr as a derivative of $\frac{\partial {\Xi }_{\text{gain}}(E)}{\partial {ϵ}_{\text{r}}}=Real\left\{2\frac{\partial {\Xi }_{\text{gain}}(E)}{\partial E}\frac{\partial E}{\partial {ϵ}_{\text{r}}}\right\}\mathrm{,}$ (12) where ${\Xi }_{\text{gain}}(E)$ is equivalent to the objective function (ξ(xi)) mentioned previously. $\partial {\Xi }_{\text{gain}}(E)/\partial {ϵ}_{\text{r}}$ is an important parameter that indicates how to change ϵr to increase ${\Xi }_{\text{gain}}(E)$. Two optimization schemes for maximizing the acceleration gradient function ${\Xi }_{\text{gain}}(E_y)$ 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 $\text{maximum}\Gamma =\frac{{\Xi }_{\text{gain}}(E)}{E_0}\mathrm{,}$ where E₀ 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 $\begin{array}{ll}\text{maximum}\hfill & \chi =\frac{{\Xi }_{\text{gain}}(E)}{E_{\text{max}}},\hfill \\ \text{subject to}\hfill & E_{\text{max}}<E_{\text{material}\text{damage}}.\hfill \end{array}$ (13) 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 and Fig. 7. The parameter settings in Fig 6 are the same as those in [30], in which the purpose is to test the solver and optimization scheme involved in this study. To accomplish this, a plane wave (E₀ 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₂), the index is n = 1.45, that is, ε₁₁ = ε₂₂ = 2.1). The cross-sectional area utilized for the optimization design region is 5β λ × λ, and the particle injection channel region is 5β λ × $\frac{1}{5}\lambda$. An acceleration gradient function (${\Xi }_{\text{gain}}(E)$) of 0.32E₀ (Γ = 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] (the achieved gradient is 0.31E₀). As calculated based on [30], Fig. 6 corresponds to an acceleration gradient of 1.37 GV/m. It can be observed from 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 (Ey/$\left|E_0\right|$), it can be observed that for higher acceleration gradients, it is useful to use dielectric mirror resonance to enhance the field in the gap.

Fig. 6Full-screen View

(Color online) The parameter settings are the same as those used in Ref. [30], The material used in the optimization process is fused silica (SiO₂). A plane wave (E₀ 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, ε₁₁ = ε₂₂ = 2.1, the square grid size (Δ x(y)) used in the simulation is 20 nm, and the achieved gradient is 0.32E₀ (Γ = 0.32).

Fig. 7Full-screen View

(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 (E₀ 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.

The parameters in Fig. 7 are the same as those used in Ref. [29], and the purpose is to highlight the advantage of using the topology optimization method. A plane wave (E₀ 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., ε₁₁ = ε₂₂ = 2.4. The cross-sectional area utilized for the optimization design region is 4β λ × λ, and the particle injection channel region is 4β λ × $\frac{1}{4}\lambda$. The acceleration gradient function (${\Xi }_{\text{gain}}(E)$) of 3.1E₀ (Γ = 3.1) is achieved by optimizing the maximum conversion factor (Γ). In Ref. [29], the gradient achieved was 0.49E₀ (Γ = 0.49). Therefore, the topology-optimized structures achieve an approximately six-fold improvement in the acceleration gradient. It is worth noting from 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 (Ey/$\left|E_0\right|$) in the y-direction, as shown in 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. The parameters used are the same as those used in Ref. [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β λ × $\frac{1}{4}\lambda$. An acceleration gradient function (${\Xi }_{\text{gain}}(E)$) of 1.4E₀ (Γ = 1.4) is achieved by optimizing the maximum conversion factor (Γ). In Ref. [29], the gradient achieved is 0.49E₀ (Γ = 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), it can also be clearly observed that the optimized structure is very different from the intuitive physical structure. According to the normalized electric field (Ey/$\left|E_0\right|$), 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/ $\left|E_0\right|$) shown in Fig. 7a) and 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.

Fig. 8Full-screen View

(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 (Fig. 7b) and Fig. 8b)), the electric field (Ey/ $\left|E_0\right|$) 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].

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 Ex that the particle feels. However, the effect of Ex 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 Ex 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.