Open Access
Add to BibSonomy
Abstract
In this paper, novel optimization methodologies of sub-relativistic guided interaction structures for dielectric laser particle acceleration (DLA) are presented. In particular, we focus on co-propagating geometries based on slot waveguides in continuous wave (CW) operation, where the particle flow and the direction of propagation of the accelerating field are co-linear. Since the velocity of sub-relativistic particles varies along the acceleration path, proper tapering of the waveguide geometry is required to achieve an extended acceleration region, and, thus a large energy gain. The design of an optimal taper ensuring particle-wave synchronicity and maximum energy gain is pursued through a physics-based approach, and these results are compared, for validation, with the outcomes of a downhill simplex method searching algorithm. Additionally, the application of a simplified 2D model of the accelerating slot waveguide is investigated and profitably used to get qualitative results useful for fast structure optimization. Indeed, this approach can hold significant potential for the development of novel accelerating structures, as it enables a thorough and fast exploration of the design space.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Scientific and industrial applications of modern particle accelerators push the research towards more compact and cheaper devices. In this framework, the accelerating electric field gradients of classical radio-frequency accelerators are mainly limited by breakdown of metals and by high peak powers at RF frequencies which are necessary for efficient acceleration [1]. However, recent progress in material science and integrated optics allowed to conceive miniaturized laser-driven interaction structures based on dielectrics [2]. Such a novel class of linear accelerators has recently attracted a lot of interest due to the larger breakdown threshold of dielectrics with respect to metals, and thanks to the availability of both powerful and stable laser sources and nano-fabrication technologies already developed for telecommunications [3].
In this paper, we focus on guided structures that exploit a well-confined accelerating mode which is perfectly superimposed and co-propagating with the particle beam. Co-propagating structures can exhibit better performance with respect to side-pumped geometries. Nevertheless, they need to be carefully engineered to maintain synchronicity between the particle beam and the guided electromagnetic field along the whole interaction length. In this respect, we compare two novel methods to design an optimal continuously tapered interaction structure. The first proposed design methodology relies on a physics-based reasoning that represents a sort of “ground-truth”; then, that method is compared with a searching algorithm based on the simplex method [4].
We also discuss the accuracy of an effective 2D electromagnetic model of the 3D accelerating structure, highlighting the clear advantages of 2D models for optimal design in terms of computational efficiency. In particular, it is recommended to use the developed 2D method for quick checks of design ideas or for exploring the design space, as it provides qualitative/approximate consistent results that require fine-tuning and 3D simulations before proceeding with any device fabrication.
Last but not least, we advocate continuous wave (CW) operation, which is an uncommon choice in this field because of the technology challenges.
The manuscript is organized as follows. In Section 2 the tapered slot waveguide under exam is presented, and the 3D-to-2D model reduction is detailed. Then, the physics-based optimization method is described in Section 3 and validated by using the optimization algorithm in Section 4. Finally, Section 5 is devoted to investigate the thermal and the time-domain behaviour of the optimized structure.
2. Slot waveguide accelerating structure
The proposed optimization approaches are discussed here with reference to the slot waveguide design illustrated in [5]. The accelerating structure consists of two high-index silicon rectangular ribs (with cross-section dimensions $a$ and $b$) sitting on a low-index silica substrate and separated by a sub-wavelength air gap of width $d$ (see Fig. 1). The geometrical parameters for an operating wavelength of $\lambda _0=2~{\mu \mathrm{m}}$ are summarized in Table 1.
Fig. 1. Slot waveguide cross section (not in scale). Silicon refractive index is $n_\mathrm {Si}=3.45$, whereas SiO$_2$ refractive index is $n_{\mathrm {SiO}_2}=1.44$.
The modes of the structure are obtained through full wave simulations at the 2D port of the 3D waveguide depicted in Fig. 1. By applying symmetry considerations, it is possible to find that the two solutions can be labeled as even and odd by analyzing symmetry of the transverse field $E_x$ with respect to $x$ [6, p. 387] [7]. In Fig. 2(a) the even solution has zero longitudinal $E_z$-field in the gap, whereas the odd solution has an appreciable longitudinal $E_z$-field in the gap, see Fig. 2(b). For the purpose of comparison, the input power of the full wave electromagnetic simulation has been set to $P_\mathrm {in}=227$ W to yield an accelerating field $E_z(x=0)=0.33$ GVm$^{-1}$ as in [5]. Notably, in actual implementations such a power level should be reduced to align with practical constraints. This will not affect the performance of the accelerating structure that can be guaranteed by extending its acceleration length.
Fig. 2. Modes profiles along $(x,y=b/2,z=\overline z)$. (a) Fundamental (even) mode and (b) accelerating (odd) mode supported by the structure in Fig. 1.
In order to exploit the longitudinal component of the waveguide mode as an accelerating field for particles along the channel, in the case of sub-relativistic particles, the waveguide parameters (e.g., the width $a$) need to be adjusted to guarantee wave-particle synchronicity (i.e. matching between phase and particle velocity) as the particle speed increases. To perform the adjustments and fulfill such a requirement, it is necessary to resort to some kind of numerical optimizations. This means that several full-wave 3D simulations of non-uniform (in $z$) slot waveguides should be performed, with the inherent drawbacks in terms of computation burden. In this scenario, approximated 2D models can be useful in order to speed up the computation time while guaranteeing a reasonable accuracy of the results.
2.1 2D simplified model
The modes supported by the guiding structure depicted in Fig. 1 can be analyzed either using rigorous 3D full-wave simulations or a 2D approximation approach, such as the effective index method (EIM) [8, §2.2]. One can choose either approach by making a trade-off between the desired level of solution accuracy and the acceptable computational speed. Methods with reduced dimensionality can indeed guarantee a trade-off between computational efficiency and accuracy, providing approximated solutions and deeper physical insights. This makes them a valuable tool for optimization algorithms requiring a fast and efficient exploration of a large parameter space, in particular when dealing with complex geometries. At the same time, rigorous 3D simulations can be used for fine tuning of the devices before fabrication.
In this contribution, we explore the EIM as a 3D-to-2D model reduction tool. In particular, we assess the accuracy of the approximation and examine its impact first on the evaluation of the effective index and then on particle dynamics results in terms of energy gain and interaction impedance.
EIM techniques have been widely applied in the past in integrated optics (see for instance [9]), also for band-gap guiding structures [10]. Nevertheless, this analysis provides novel insights into the potential of this approach in the specific framework of particle acceleration.
The application of the EIM allows to reduce the actual 3D structure composed of silicon and silica into an equivalent 2D structure composed of materials with properly derived effective indices. In order to determine the new effective material, two modal simulations were performed as schematized in Fig. 3. In the first simulation, a 3-layer structure along $y$ composed of air-silica-air was considered and, by solving for the $E_x$ polarization, an effective index $n_\mathrm {eff1}=1.17926$ was obtained. The structure in the second simulation was composed of 4 layers, i.e., air-silicon-silica-air, and the evaluated effective index for $E_x$ polarization was $n_\mathrm {eff2}=2.59511$. Based on these results, the equivalent 2D model can be built by properly juxtaposing the effective media along the $x$-axis, see Fig. 4.
It is worth to note that, while the effective index $n_\mathrm {eff2}$ is well defined primarily by the precise dimension $b$, the effective index $n_\mathrm {eff1}$, obtained for the first region, is somewhat more arbitrary as it depends on the parameter $h$ that may vary. Nevertheless, it was noticed that the qualitative agreement between the 2D and 3D models remains satisfactory also when varying $h$. However, the fabrication process may require a completely different stack-up, such as placing a thicker silica layer on top of silicon and different strategies may be adopted to fix $n_\mathrm {eff1}$. For instance, one may explore the introduction of ‘empirical’ effective index values for the slot region to enhance the 2D model accuracy.
Comparison of mode profiles of the 3D and the equivalent 2D structures can be used for a first assessment of the EIM approach. In Fig. 5, we compare the mode profiles evaluated at the 2D port and at the 1D port of the 3D-full and 2D-reduced models, respectively. As it can be seen, the modes fit very well and their striking profile similarity is evident, in particular for the accelerating mode.
Fig. 5. Comparison between 3D and 2D equivalent mode solutions. (a) $E_x$ and (c) $E_z$ component of the even mode; (b) $E_x$ and (d) $E_z$ component of the odd mode.
It is worth noting that in order to reach the desired value of the accelerating field $E_z(x=0)=0.33 \; \textrm{GV/m}$ inside the gap, an input power per unit length $P_\mathrm {in}'=0.75$ GWm$^{-1}$ for the 2D model should be fixed. Alternatively, one could consider a quasi-2D model by properly setting the $\Delta y$ parameter of the EIM-based reduced model in Fig. 4 so that $P_\mathrm {in}'=P_\mathrm {in}/\Delta y$. For the case at hand, $\Delta y=0.303$ μm and the quasi-2D simulation requires top and bottom Perfect Magnetic Conductor (PMC) boundaries which are compatible with the $E_x$ polarization and provide $y$-invariant solutions.
Qualitative agreement of the results between 3D and 2D model has been shown in Fig. 5, nevertheless a quantitative estimation of the accuracy of the 2D reduced model requires the evaluation of carefully selected metrics. To this aim, we performed a comparison in terms of effective index, $n_\mathrm {eff}$. Table 2 shows the effective indices relative to the two models. Despite the differences, which could be significant in integrated optics, these findings suggest that 2D reduced models may have the potential to predict the performance of the actual 3D model, as it will be shown in the following Section. In particular, its usefulness and significance lies in its ability to facilitate rapid exploration of the design space and can be used to test new geometries and identify a suitable working point.
Table 2. Quantitative comparison between the 3D study and the reduced 2D analysis.
In a practical implementation, several significant physical constraints should be considered when determining the ‘operating point’. These factors include: (a) laser pump availability; (b) the potential to up-shift the operating wavelength and the exploration of alternative materials/integrated photonic platforms to mitigate the two-photon absorption (TPA) effect in silicon [11]; (c) the feasibility of using different lower index materials to accommodate larger core sizes and reduce peak optical intensities.
In Table 2 we also report the interaction impedance, $Z_c$, a popular metric for DLAs. The interaction impedance is a local property that pertains to a specific cross-sectional area at a given $z$-position and is defined as:
where $E_\mathrm {acc}$ is the average axial electric field given by $E_\mathrm {acc}=\frac {1}{\delta l}\int _{-\delta l/2}^{\delta l/2}E_z(z)dz$, $E_z(z)$ being the axial accelerating field; $\delta l$ is the waveguide length and $P_\mathrm {in}$ is the power coupled into the structure.For a 2D equivalent model we can modify the definition (1) to consider the input power per unit length $P_\mathrm {in}'$:
Alternatively, one can exploit the previously defined fictitious equivalent thickness $\Delta y$ that, for the same input power $P_\mathrm {in}=P_\mathrm {in}'\Delta y$, ensures the same $E_z=0.33$ GV/m and thus the same $Z_c$. For the model at hand, $Z_c$ is equal to 1.92 k$\Omega$ in both the 3D and 2D$\times \Delta y$ case.
3. Physics-based optimization approach
The interaction impedance $Z_c$ is often employed in the literature to compare different DLA structures. However, it is not a suitable figure of merit (FoM) for designing an extended sub-relativistic accelerating structure where the waveguide geometry must be finely tuned to guarantee synchronicity between particles and accelerating field. In this respect, a physics-based optimization strategy can rely on the maximization of a more suitable FoM for the structure under exam, i.e. the energy gain between output and input section of the accelerator.
By assuming that the particle path in the gap is located at $y=\frac {b}{2}$, i.e. the generic point of the path is $P\equiv (0,\frac {b}{2},z)$ with $z \in [0,L]$, $L$ being the overall length of the device, the energy gain can be defined as:
where $E_z(z,t)$ is the accelerating field along the $z$-direction of the reference frame, and $q$ is particle charge.The particle moves along this path with a $z$-dependent velocity, $v_\mathrm {part}(z)$. Hence, from:
it arises that:As it was argued in the Introduction, the synchronicity between the particle beam and the electromagnetic field guided along the path can be guaranteed by properly engineering the transverse geometry of the waveguide. In our case, the ribs composing the slot waveguide can be tapered in order to fulfill this requirement. Notably, the effective index of such a “tapered” waveguide is also $z$-dependent, so the accelerating field in (3) can be expressed as [12, Eq. (2.1), p. 32]:
wherein, indicating with $k_0$ the wave-number of vacuum, we have:The integral in (3) can be evaluated along the particle path by using the previous equation and the time expression given by Eq. (5) to get the final formulation for the FoM:
The proposed physics-based procedure aims at finding the optimal geometrical properties of the accelerating structure in order to maximize the output energy gain (9). In particular, even if the illustrated technique is quite general, optimization will be performed by varying only the width $a$ of the silicon ribs due to fabrication constraints.
The design procedure is based on the following statement: if the wave phase velocity $v_\mathrm {ph}$ equals the particle velocity $v_\mathrm {part}$ for each $z \in [0,L]$, the term within the square brackets in (9) becomes zero, the cosine function one, and the accelerating field seen by the particle is maximized (i.e., for the particle, Eq. (6) simplifies to $E_z=E_0$). Thus, we envisage a simple but effective way to impose the desired synchronicity condition, that is:
relying on the derivation of an optimal profile $a(z)$ of the taper for the guiding structure.The optimal design can be implemented by following three steps that are detailed in the next subsections.
3.1 Step 1: Beam dynamics (BD) evaluation
Given the initial values for particle energy and accelerating field $E_0$, the particle velocity evolution $v_\mathrm {part}(z)$ is calculated by solving the relativistic dynamics under the assumption (10), as detailed in Appendix A. By following [5], we considered a 80 keV input energy, corresponding to an initial normalized velocity $\beta =v_\mathrm {part}/c_0=0.5$, and $E_0=0.33$ GV/m. Due to assumption (10), the cosine term in (9) equals one, i.e. the particle experiences the maximum possible accelerating field, $E_0$. In this step, for the relativistic dynamics calculations, an average constant value can be used for $E_0$, since in common structures it changes slowly with the considered $\overline z$-section. This approximation can also be readily refined because, at each $\overline z$-section, $E_0(\overline z)$ depends on the cross-sectional parameters (in our case on $a(\overline z)$) in a known manner. Hence, the current value of $E_0$ could be employed to eliminate the approximation.
The behavior of the normalized particle velocity along the path, $\beta (z)$, under the above settings is shown in Fig. 6(a) for $L=30\,{\mathrm{\mu} \mathrm{m}}$. Notably, this best case scenario predicts an output energy equal to 90.14 keV, corresponding to an energy increase of about 12%.
Fig. 6. Theoretical BD result in a slot waveguide with initial particle energy equal to 80 keV and $E_z(\bar z=0,t)=E_0=0.33$ GV/m, in case of: (a) shorter ($L=30\,{\mathrm{\mu} \mathrm{m}}$) and (b) longer ($L=300 \, {\mathrm{\mu} \mathrm{m}}$$L=300 \, {\mathrm{\mu} \mathrm{m}}$) length.
The evolution of $\beta (z)$ along the particle path can be directly translated into the required effective index for each waveguide section, as $n_\mathrm {eff}(z)=1/\beta (z)$. Figure 6(a) also shows the optimal profile of $n_\mathrm {eff}(z)$ (see the right axis). It can be observed that, as the velocity increases the effective index should decrease in order to keep matching between phase and particle velocity.
In longer structures deviations from approximately constant acceleration and linear normalized velocity $\beta (z)$ are expected due to the relativity constraint. Therefore, we repeated the analysis for a larger length. The obtained curve in Fig. 6(b) shows an energy increase of more than 100%, and a saturation effect for the normalized velocity $\beta (z)$ increase becomes visible.
3.2 Step 2: Derivation of a calibration curve
The sub-relativistic $z$-dependent normalized velocity $\beta (z)$ determines the required effective index profile $n_\mathrm {eff}(z)$ that is also shown in Fig. 6 (see right secondary axis) and it is utilized to extrapolate the optimal tapering profile of the structure at hand. As a matter of fact, given the particle velocity $\beta (z=\bar z)$ at each $\bar z$-section of the waveguide, it is possible to determine the corresponding width of the silicon ribs $a(z=\bar z)$ that guarantees the required phase velocity of the accelerating field.
In this respect, we carried out a parametric study in COMSOL Multiphysics by considering a generic slot waveguide cross-section with $a \in [0.3 - 0.8] $ $\mathrm{\mu}$m and steps of $0.01$ $\mathrm{\mu}$m. Notice that the computational domain was reduced by simulating only half structures by taking advantage of mirror symmetry of the geometries under exam. Hence, the odd accelerating mode can be selected by placing PMC boundary condition at $x=0$ (the even mode can be studied by using a Perfect Electric Conductor (PEC) boundary condition).
This modal analysis was performed first for the reduced 2D model depicted in Fig. 7(a), which was derived by applying the EIM as in Fig. 3. In this case, the effective index $n_\mathrm {eff}$ as a function of $a$ was evaluated by solving for the 1D-port modes. Then, the study was performed again for the complete 3D model shown in Fig. 7(b) by evaluating the effective index from calculation of the 2D-port modes.
Fig. 7. Schematic views of the models simulated in COMSOL. (a) 2D reduced model. (b) Actual 3D model. Half structures can be studied by using PMC boundaries.
In Fig. 8(a) we show the calibration curves for both the reduced model (2D) and the actual model (3D). Notice that both these curves are less precise for small values of $a$ (i.e. $a \lesssim 0.45$) as the waveguide approaches the cut-off condition, while they fit very well in the remaining part.
Fig. 8. (a) Calibration curve for the reduced 2D model (continuous line) and for the actual 3D model (dashed line). (b) Optimal taper profile for a slot waveguide derived from the reduced 2D model (continuous line) and from the actual 3D model (dashed line). The shaded area indicates the region where the 3D waveguide is approaching the cutoff.
The curves $n_\mathrm {eff}(a)$ plotted in Fig. 8(a) can be inverted to get $a(n_\mathrm {eff})$ and these functions can be used to find the optimal width of the silicon ribs, $a(z)$, for each $z$-section along the waveguide by using the expression:
The optimal taper profiles for the 2D and the 3D accelerating structures are depicted in Fig. 8(b). As it can be observed, the designed tapers exhibit the same qualitative trend, but for a slight offset, until the rib width approaches the cut-off value, i.e., $a=0.45\, {\mathrm{\mu} \mathrm{m}}$, which corresponds to $L=120\,{\mathrm{\mu} \mathrm{m}}$.
The effect of this discrepancy on the acceleration properties will be investigated in the following subsection, with particular emphasis on the evaluation of the selected FoM.
3.3 Step 3: Testing of the taper
The propagation of the accelerating field in the tapered slot waveguide was simulated to assess the performance in terms of energy gain. This step must be performed as a numerical validation of the physics-based optimization procedure presented in the first two steps.
The slot waveguide in Fig. 7 has been re-designed with a fixed (minimum) slot size, $d$, and a tapered rib width $a$. The profile of the rib width, $a$, shown in Fig. 8(b) in the interval $[0,120]\,{\mathrm{\mu} \mathrm{m}}$, has been modeled in COMSOL by using a polynomial in the variable $s \in [0,1]$:
in which $p_3=0$ in case of $L=30 \,{\mathrm{\mu} \mathrm{m}}$, while $p_3\ne 0$ in case of $L=120 \, {\mathrm{\mu} \mathrm{m}}$. The actual width, $a(z)$, is given by $a(z)=a\left (s=\frac {z}{L}\right )$.The $2^{nd}$ order polynomial coefficients for the optimal taper profiles of Fig. 8(b) ($L=30 \,{\mathrm{\mu} \mathrm{m}}$) are reported in Table 3. However, fitting of the optimal taper for the longer structure (e.g., $L=120 \,{\mathrm{\mu} \mathrm{m}}$) requires a $3^{rd}$ order polynomial (see Table 4).
Table 3. $2^{nd}$ order polynomial coefficients fitting the optimal taper profiles in Fig. 8(b) when $z \in [0, 30]$μm is considered.
Table 4. $3^{rd}$ order polynomial coefficients fitting the optimal taper profiles in Fig. 8(b) when $z \in [0, 120]$μm is considered.
COMSOL full wave electromagnetic simulations were used to validate the physics-based optimization procedure, both in the 3D and 2D case. First, the longitudinal electric field $E_z(z, t)$ was extracted from COMSOL and the behavior of its phase, $\angle {E_z}$, was compared with the expected ideal evolution:
Then, the particle BD (with initial energy of $80$ keV) was evaluated by using the accelerating field taken from the full-wave simulation. Notice that the input power was fixed, in agreement with the previous sections and for comparison with [5], to $P_\mathrm {in}=227$ W for the 3D model and to $P_\mathrm {in}'=0.75$ GWm$^{-1}$ for the 2D model.
As it can be seen in Fig. 9(a), when the optimal taper profile arising from the reduced 2D model is used for a 2D full-wave simulation, the particle is continuously accelerated and the output energy (i.e., the input energy plus the FoM) is $91.30$ keV, which is close to the theoretical value ($90.14$ keV). Furthermore, excellent agreement between simulated and theoretical phase of the accelerated field is clearly visible (see Fig. 9(b)). Conversely, when the same profile coming from the 2D model is used for tapering the actual 3D slot waveguide, the behavior with dashed line in Fig. 9 is obtained. In this case the particle is accelerated only in the first $\sim 11\,{\mathrm{\mu} \mathrm{m}}$ of the structure and then it is decelerated, since synchronicity between particle and accelerating field is not guaranteed, as it is possible to verify in Fig. 9(b) by referring to the green dashed line with circular markers. Instead, by repeating the analysis by using the optimal taper profile evaluated from the 3D model in the actual 3D structure, the phase synchronicity is restored and the particle is correctly accelerated along the entire waveguide (see Fig. 10), by reaching an output energy of $89.58$ keV.
Fig. 9. (a) Comparison of the normalized particle velocity achieved with the optimal tapering plotted with a continuous line in Fig. 8(b) ($L=30\, {\mathrm{\mu} \mathrm{m}}$), derived from the 2D model and implemented: in the reduced 2D model (continuous line) and in the actual 3D model (dashed line). (b) Comparison of the accelerating field phase between: theoretical derivation (continuous line), simulation in the reduced 2D model (triangular markers), simulation in the actual 3D model (circular markers).
Fig. 10. (a) Normalized particle velocity achieved with the optimal tapering plotted with a dashed line in Fig. 8(b) ($L=30 \,{\mathrm{\mu} \mathrm{m}}$), derived from the 3D model and implemented in the actual 3D model. (b) Comparison of the accelerating field phase between: theoretical derivation (continuous line) and simulation in the actual 3D model (triangular markers).
Notwithstanding the slight discrepancies between the effective indices of the 3D and the 2D structures, the acceleration properties are estimated in a satisfactory fashion. A careful analysis of Table 3 reveals that the optimal tapers evaluated from 3D and reduced 2D model are characterized by a different coefficient $p_0$, which is directly related to the width $a$ at $z=0$, whereas it is noteworthy that $p_1$, $p_2$ (and $p_3$) are similar. This suggests that there is an offset between the 3D and 2D calculations of the width $a$, which emerges also from the qualitative analysis of Fig. 8(b), that can be taken into account at the initial stage of the design.
As a final assessment, we performed a full wave electromagnetic simulation of the longer structure ($L=120 \,{\mathrm{\mu} \mathrm{m}}$) by using the optimized profiles in Fig. 8(b). The corresponding evolution of the particle velocity along the path is shown in Fig. 11(a) and the output energy is 122 keV which is still close to the theoretical value (120 keV). It is straightforward to notice that the particle is continuously accelerated within the structure, and the resulting velocity does not increase smoothly but with oscillating slope (acceleration) because the phase matching is sub-optimal. This phase mismatch is due to the fact that the optimal taper profile has been derived by assuming perfect synchronicity between particle and accelerating field, as well as constant electric field $E_0$ experienced by the traveling particle. However, in such a long waveguide the width of the rib is reduced of about 20% at the output of the structure with respect to input due to tapering. As a matter of fact, this leads to an increase in the field intensity within the gap, since by reducing $a$ the waveguide gets closer to cut-off by generating longer field tails in the air gap region. In a nutshell, when the tapered slot waveguide is very long, the accelerating field as seen by the particle cannot be considered constant along the extended structure (see Fig. 11(b)), thus the physics-based approach introduces an approximation which gets worse as long as length increases. However, it is still useful to identify a suitable working point and systematic errors can be rectified through a final adjustment using full wave simulations.
Fig. 11. (a) Normalized particle velocity achieved with the optimal tapering plotted in Fig. 8(b) ($L=120\,{\mathrm{\mu} \mathrm{m}}$), derived from the 2D model and implemented in the reduced 2D model. (b) Amplitude of the accelerating field along the optimized tapered waveguide.
The limitation highlighted above can be easily circumvented in at least two manners: (i) by resorting to a sequential optimization of shorter waveguides where we update the value $E_0$ in the design procedure (e.g., by considering the correct axial electric field $E_0$ for each section under exam); (ii) by deriving a calibration curve for $E_0(a)$ with the same procedure carried out for the calibration curve $n_\mathrm {eff}(a)$ of Fig. 8(a) and evaluating the correct $E_0(a(z))$ in the first BD step.
4. Nelder-Mead method optimization
A numerical optimization procedure based on the Nelder-Mead (NM) method was applied in order to have a benchmark for the evaluation of the effectiveness of the physics-based approach presented in the previous Section.
The NM method (also called downhill simplex method) [4] is an algorithm used to find the optimum (minimum or maximum) of an objective function or a FoM in a multidimensional space. It is a direct-search algorithm, i.e., it requires only function evaluations, and neither gradient nor derivative approximations of the FoM are computed. Hence, it can be applied to non-linear optimization problems for which derivatives may not be known. Its main drawback is related to the slow convergence rate.
The NM method in $n$-dimensions is based on a set of $n + 1$ test points arranged as a simplex, the geometrical figure consisting, in $N$ dimensions, of $N+1$ vertices (e.g., a triangle in two dimensions, a tetrahedron in three dimensions). Once the FoM in each vertex of an initial simplex is computed, the NM method takes a series of iterative steps: reflections of the vertex with the worst FoM value through the opposite face of the simplex; expansions or contractions of the simplex in one or all directions; and so on, until the termination criterion is reached (e.g., when the vector distance in the current step is fractionally smaller in magnitude than some tolerance).
It is important to note that the NM method is a local search algorithm, i.e., it finds the FoM local optimum closest to the starting simplex.
In this Section, the same FoM in Eq. (3) considered for the physics-based approach in Section 3 is optimized by means of the NM method. More in detail, the optimization aims at finding the polynomial coefficients describing the taper, and just $p_0$ is considered known and fixed to the value obtained with the previous method. Conversely, the others parameters $p_1$, $p_2$ and $p_3$ (only for the case with $L=120 \,{\mathrm{\mu} \mathrm{m}}$) are optimized within the given ranges reported in Table 5 and 6, defining the search space.
Table 5. $L=30$μm waveguide: design variables range
Table 6. $L=120$μm waveguide: design variables range
The NM algorithm we used for the optimization is implemented in the Scipy library [13]. Since for each FoM evaluation a full-wave simulation is required, the strategy illustrated in Fig. 12 is applied to allow the NM method to automatically perform the whole optimization procedure. In particular, the NM algorithm chooses the coordinates $p_1$, $p_2$ and $p_3$ (only for the case with $L=120 \,{\mathrm{\mu} \mathrm{m}}$) of a vertex and, by means of a batch script, calls COMSOL to run the full-wave simulation of the model in Fig. 7(a). At the end of the numerical analysis, COMSOL returns the accelerating field along the $z$-axis in a file which is used by another batch script to compute the FoM (3). The energy gain is then given back to NM and the procedure is repeated until convergence is achieved.
The optimization was performed starting from a random simplex in the search space. The NM optimal $a(z)$ profiles are shown with a dashed blue line in Fig. 13(a) for $L=30 \,{\mathrm{\mu} \mathrm{m}}$ and Fig. 13(b) for $L=120 \,{\mathrm{\mu} \mathrm{m}}$, respectively, superimposed with the curves obtained through the physics-based approach (continuous red line). The corresponding optimal polynomial coefficients, as well as the optimized output energy, are summarized in Tables 7 and 8 and have been reached after $219$ and $311$ iterations, respectively. In both cases the two optimization methods gave very similar results for the initial part of the taper, while the curves become different when the output is approaching (see Fig. 13). This result is interesting as it highlights the fact that while the physics-based approach is based on the assumption of constant $E_z$ field experienced by the particle, the numeric optimization aiming at maximizing the FoM takes into account the field variation along the $z$-direction intrinsically. As such, this allows to reach a higher output energy ($133.4$ keV vs $122$ keV) and to furnish a more stable and monotonic acceleration to the particles, as confirmed by the plot in Fig. 14.
Fig. 13. Comparison between the optimal taper profiles obtained with the physics-based approach (continuous line) and from the Nelder-Mead numerical approach (dashed line) in case of (a) $L=30 \,{\mathrm{\mu} \mathrm{m}} $ and (b) $L=120 \,{\mathrm{\mu} \mathrm{m}}$.
Fig. 14. Comparison of the obtained acceleration profile and output energy with the different optimization strategies.
Table 7. Coefficients of the 2$^\mathrm {nd}$ order polynomial describing the NM optimal waveguide tapering when $L=30$μm.
Table 8. Coefficients of the 3$^\mathrm {rd}$ order polynomial describing the NM optimal waveguide tapering when $L=120$μm.
5. Thermal and time domain considerations
5.1 Thermal behavior
DLAs provide highly efficient acceleration in CW regime, therefore the solution of the steady-state problem, in which a CW signal is injected into the slot waveguide, is considered in this section. The objective of this analysis is to find the steady-state temperature of the optimized slot waveguide when the input power is equal to $227$ W (which is required to obtain the accelerating field gradient of $0.33$ GV/m, see the above sections).
To this aim, a 3D COMSOL model of the structure with the optimal-valued taper has been employed to solve the stationary heat transfer equation coupled to frequency-stationary Maxwell equations. A loss tangent similar to that of low-loss communication optical fibers (attenuation $0.3$ dB/km) [3,14] has been utilized to take into account volume losses in the slot waveguide. Possible losses introduced by surface roughness and scattering [15] or TPA [11] are not included in this analysis.
The Heat equation (Fourier’s Law) is solved with a heat source $Q=Q_0$ generated by the electromagnetic power loss density within the silicon ribs and silica substrate. For silicon, the heat capacity is $C_{Si}$= 700 J/(kg$\cdot$K) and the thermal conductivity is $k_{Si}$=130 W/(m$\cdot$K), whereas for silica the heat capacity is $C_{SiO_2}$= 703 J/(kg$\cdot$K) and the thermal conductivity is $k_{SiO_2} = 1.38\,\mathrm{W}/(\mathrm{m} \cdot \mathrm{K})$. A temperature node on the bottom of the silica substrate has been used to impose a fixed temperature $T=T_0=20^\circ$ C. Notice that this condition can be readily achieved by using a heatsink that is connected to a Peltier junction. On the surfaces of the waveguide, a “Surface-to-Ambient Radiation” boundary condition $q=\sigma \epsilon (T_\mathrm {amb}^4-T^4)$ has been set for the Radiative Heat Flux $q$, where $\sigma$ is the Stefan-Boltzmann constant, $\epsilon$ is the surface emissivity equals to 0.7 [16], and T$_\mathrm {amb}$ is the room temperature taken as the external reference temperature (see Fig. 15 (a)).
Fig. 15. Thermal model: Simulated geometry and boundary conditions (a), temperature distribution in the volume without (b), and with (d) the air box over the cladding; temperature profile along $(x,y=b/2,z=\overline z)$ (c).
The temperature distribution in the volume without (Fig. 15 (b)) and with (Fig. 15 (d)) the air over-cladding shows that the temperature keeps almost constant and equal to the prescribed value (set by the Peltier thermocouple). However, a minor increase of the temperature is observed in correspondence of ribs and gap, as it is shown in Fig. 15 (c). These simulations confirm that, in principle, thermal effects do not hinder CW operation of an optical accelerator [14]. As expected, losses caused by absorption in the material are negligible.
5.2 Time-domain considerations
It is well-known in the accelerator community that metallic structures operating at room temperature are limited in repetition rate to tens of Hz. Conversely, low-loss DLAs can potentially work at MHz repetition rates, or even in CW, without resorting to expensive superconductive technologies. This is the reason why we focused our treatment to CW operation, therefore group velocity and dispersion are not issues. However, some considerations related to the performance of guided-wave DLAs when pumped with short and even ultrashort pulses can be useful when working in different scenarios.
The temporal behavior of tapered slot waveguides can be investigated by evaluating group velocity and dispersion of the guiding structures. In this scenario, the low group velocity which is achievable in highly confining structures can be a limiting factor, whereas the effect of dispersion appears to be negligible. In order to quantify this phenomenon, we calculated the transit time of a single electron along the designed DLAs, $T_\mathrm {part}$, by integrating the particle velocity profiles in Fig. 6. It results that $T_\mathrm {part}$ is equal to $194$ fs for the $30$-$\mathrm {\mu }$m-long waveguide, whereas the transit time is $0.74$ ps for the $120$-$\mathrm {\mu }$m-long DLA. Instead, the group delay experienced by the pulsed accelerating field has been determined by resorting to modal simulations along the waveguides, where we neglected material dispersion, with the goal of estimating the group velocity profile as a function of the parameter $a$. Then, the group delay as a function of $z$ is evaluated as:
We call $T_g=t_g(z=L)$ the group delay at the output of a waveguide, which should be compared with the electron transit time $T_\mathrm {part}$. In our case, $T_g$ is equal to $454$ fs for the $30$-$\mathrm {\mu }$m-long waveguide, whereas the output group delay is $1.73$ ps for the $120$-$\mathrm {\mu }$m-long DLA.
These results suggest some important qualitative considerations. First, notice that particles are systematically way faster with respect to the pulsed accelerating field, therefore during propagation along the gap of the slot waveguide particles tend to progressively lose alignment with the peak of the pulse. In details, the accelerating field is delayed by about $260$ fs with respect to the electron at the output of the $30$-$\mathrm {\mu }$m-long waveguide. As a consequence, the effect of group velocity can significantly degrade the acceleration performance if the pulse duration is not larger or comparable to this delay. At the output of the $120$-$\mathrm {\mu }$m-long waveguide temporal shift between electromagnetic field and particle is about $1$ ps, therefore it is clear that a comparable longer input pulse is required for efficient particle acceleration in such a long structure.
A quantitative analysis of temporal effects related to particle acceleration in DLAs operating in pulsed regime requires full-wave time-domain simulations which are beyond the scope of the present work. In particular, it would be interesting to understand what happens when shortening the pulse duration with fixed energy. Indeed, shorter pulses are characterized by larger peak intensity of the accelerating field. However, in this scenario, longer sections of the waveguide do not contribute to the acceleration process due to temporal delay between pump and particle. We thus envisage that a trade-off between pulse duration and peak power could be identified.
6. Conclusion
In this paper, novel optimization methodologies of DLAs have been presented, with particular focus on slotted-waveguide structures for sub-relativistic acceleration regimes in CW operation, tapered in order to guarantee synchronicity between the phase velocity of the accelerating wave and the particles velocity as they travel along the accelerating channel, gaining energy in the process. The optimal tapered profile of the slotted waveguide, which ensures maximum energy gain, has been obtained through a physics-based approach, consisting in the derivation of a calibration curve that relates the particle velocity, $\beta (z)$ with the optimum rib width, $a(z)$, ensuring the matching between the phase and the particle velocity along the accelerating channel. All the calculations have been done by employing a 3D structure model and a simplified 2D model, evaluated through the effective index method, that ensures lower simulation complexity and thus fast simulation time (about 2 seconds for the 2D model with respect to about 15 minutes required for the full 3D model). The availability of this fast 2D model-based approach enables the adoption of iterative methods for optimization; indeed, a numerical procedure based on the Nelder-Mead method has been applied to validate the physics-based approach and results in terms of calibration curve profile are in excellent agreement between the two methods.
In perspective, the presented simplified 2D approach could allow the optimization, also through inverse design methods [17], of more complex structures, such as Photonic Crystal (PhC) [18] based DLAs. The reduced 2D model enables rapid exploration of the design space and opens the path to optimization of unconventional, non-intuitive geometries. Moreover, fast direct solvers and inverse design methods can be applied to optimize the transition from sub-relativistic [5] to relativistic [18] sections/stages of particle accelerators.
Appendix A: Relativistic dynamics
Starting from paper [19], in which the acceleration of a single particle (in particular a proton) is considered, we exploit a different way to integrate the equation of motion [19, Eq. (3)].
In the case at hand the particle moves along the $z$-direction, so all equations concerning its relativistic dynamics are scalar. For the sake of simplicity, from now on the subscript $z$ is dropped and all the variables are intended as $z$-components.
The electron relativistic mass, $m$, is related to the rest mass, $m_0$, by the expression $m=\gamma m_0$, where $\gamma$ is the Lorentz factor:
The momentum $p$ of the electron is given by: By inverting Eq. (16) for the particle velocity, it is possible to get: Since $v=\frac {d z}{dt}$ we have: When the equation of motion is also considered: where $q$ is the electron charge and $E(z,t)$ is the accelerating field, a system of two equations is derived:Funding
INFN INFN V National Committee through the MICRON grant; Ministero dell'Università e della Ricerca as part of the PRIN 2022 projects DOSE (2022MANFK5) and GRACE6G (2022H7RR4F).
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. E. I. Simakov, V. A. Dolgashev, and S. G. Tantawi, “Advances in high gradient normal conducting accelerator structures,” Nucl. Instrum. Methods Phys. Res., Sect. A 907, 221–230 (2018). [CrossRef]
2. R. J. England, D. Filippetto, G. Torrisi, et al., “Laser-driven structure-based accelerators,” arXiv, arXiv 2203.08981 (2022). [CrossRef]
3. R. J. England, R. J. Noble, and K. Bane, “Dielectric laser accelerators,” Rev. Mod. Phys. 86(4), 1337–1389 (2014). [CrossRef]
4. J. A. Nelder and R. Mead, “A simplex method for function minimization,” Comput. J. 7(4), 308–313 (1965). [CrossRef]
5. Z. Zhao, T. W. Hughes, S. Tan, H. Deng, N. Sapra, R. J. England, J. Vuckovic, J. S. Harris, R. L. Byer, and S. Fan, “Design of a tapered slot waveguide dielectric laser accelerator for sub-relativistic electrons,” Opt. Express 26(18), 22801–22815 (2018). [CrossRef]
6. A. W. Snyder and J. D. Love, Optical waveguide theory, (Chapman & Hall, 1983), Vol. 175.
7. A. W. Snyder and W. R. Young, “Modes of optical waveguides,” J. Opt. Soc. Am. 68(3), 297–309 (1978). [CrossRef]
8. K. Kawano and T. Kitoh, Introduction to Optical Waveguide Analysis: Solving Maxwell’s Equation and the Schrödinger Equation (John Wiley & Sons, 2004).
9. K. S. Chiang, “Analysis of optical fibers by the effective-index method,” Appl. Opt. 25(3), 348–354 (1986). [CrossRef]
10. M. Qiu, “Effective index method for heterostructure-slab-waveguide-based two-dimensional photonic crystals,” Appl. Phys. Lett. 81(7), 1163–1165 (2002). [CrossRef]
11. L. Zhang, A. M. Agarwal, L. C. Kimerling, and J. Michel, “Nonlinear Group IV photonics based on silicon and germanium: from near-infrared to mid-infrared,” Nanophotonics 3(4-5), 247–268 (2014). [CrossRef]
12. T. P. Wangler, RF Linear Accelerators (John Wiley & Sons, 2008).
13. P. Virtanen, R. Gommers, and T. E. Oliphant, “Scipy 1.0: fundamental algorithms for scientific computing in python,” Nat. Methods 17(3), 261–272 (2020). [CrossRef]
14. V. Karagodsky, A. Mizrahi, and L. Schächter, “Thermal scaling laws of the optical Bragg acceleration structure,” Phys. Rev. Spec. Top.--Accel. Beams 9(5), 051301 (2006). [CrossRef]
15. X. Li, X. Feng, X. Xiao, K. Cui, F. Liu, and Y. Huang, “Experimental demonstration of silicon slot waveguide with low transmission loss at 1064nm,” Opt. Commun. 329, 168–172 (2014). [CrossRef]
16. “Laser Heating of a Silicon Wafer,” COMSOL Gallery, Application ID: 13835, https://www.comsol.com/model/laser-heating-of-a-silicon-wafer 13835.
17. N. V. Sapra, K. Y. Yang, D. Vercruysse, K. J. Leedle, D. S. Black, R. J. England, L. Su, R. Trivedi, Y. Miao, O. Solgaard, R. L. Byer, and J. Vučković, “On-chip integrated laser-driven particle accelerator,” Science 367(6473), 79–83 (2020). [CrossRef]
18. G. Torrisi, G. Mauro, L. Celona, D. Mascali, S. Gammino, G. Sorbello, C. De Angelis, and A. Locatelli, “Numerical study of photonic-crystal-based dielectric accelerators,” in Journal of Physics: Conference Series, 1350 (IOP Publishing, 2019), p. 012060.
19. L. Pálfalvi, J. Fülöp, G. Tóth, and J. Hebling, “Evanescent-wave proton postaccelerator driven by intense THz pulse,” Phys. Rev. Spec. Top.--Accel. Beams 17(3), 031301 (2014). [CrossRef]
