1. Introduction

Improvement upon current computing architecture is required in order to enhance computation power, time and resources for industries such as medicine, catalysis, and finance. As compared to digital computers, optical computers offer the capabilities of parallel processing, high speed, lower power consumption, and operation at multiple frequencies [1,2] thus exploiting a photonic system can make these devices a prime candidate to realize this possibility [3], especially in the least exploited THz regime. In order to achieve an integrated analog optical computing system, it is necessary to focus the incident wave to the next operational computing component which can be achieved by the use of computational metamaterials. Numerous systems, designed for free space, have been proposed but have been difficult to integrate and bulky in design [4–7]. One such metamaterial known as the gradient index (GRIN) lens, can alleviate these issues through being an on-chip device. A prototypical GRIN lens consists of a series of micro-layered structures with holes of different geometries to vary the index of refraction [8,9] that manipulate the propagation of electromagnetic waves [10]. Previous works of GRIN lens metamaterial structures have theoretically and experimentally demonstrated to operate in the THz regime [11–13] with strong focusing capabilities [14–16] and tunability [17–19].

Recently, GRIN lenses have been integrated as components of larger photonic devices, consisting of a multitude of subblocks, to realize a new optical computing technology, the quantum algorithm emulator (QAE). QAEs simulate quantum search algorithms with classical waves via the superposition principle and interference phenomena which is integral to rapid searching and solving difficult problems that would be time and power consuming on digital computers [4,20,21]. QAEs have only been explored in the microwave region where the measured electric field amplitude is the probability amplitude of the equivalent quantum state [21]. Zhang et al. proposed a dielectric device, made of Veroclear810, consisting of an oracle subblock, two Fourier transform (FT) subblocks, and a phase plate subblock [4,21]. The oracle subblock imprints a spatially phase dependent profile on the incoming wave while the FT subblock and phase plate subblock converts the phase difference for the oracle subblock to amplitude information [21]. Through this device they were able to simulate Grover’s Algorithm and show that the number of iterations performed on the device was consistent with the efficiency of the quantum search algorithm.

Metamaterials have reached a high degree of maturity and have recently emerged as a new approach for quantum computing. In this context, Wei et al. proposed a quantum searcher of an on-chip silicon device that consisted of four metasurfaces: an oracle metasurface, two metalenses, and a middle metasurface where different spatial positions of the incident wave showed repeatability in the distribution of the output wave intensity [4]. Cheng et al. experimentally verified a Deutsch-Jozsa (DJ) algorithm with a millimeter scaled all-dielectric device [20].

Here, we report a new design for a quantum algorithm emulator in the THz frequency regime based on a simple platform and optimized by machine learning. The investigated device is composed of a microstructured oracle subblock made of silicon substrate and a FT subblock made of Kapton polyimide. The Kapton film chosen acts as a 2D photonic crystal (PhC) slab and proven to exhibit strong electric field confinement and interaction with THz waves with minimal absorption loss [22]. In this work, numerical simulations enhanced by machine learning (ML) were applied to optimize the hole radii and thickness of the FT subblock to achieve an optimal distribution of wave intensities. The structural design of the QAE is first presented along with the numerical analysis showing the initial optimized output. Lastly, the process and results for the ML process are discussed. Our aim is that the simple design of our device and its compact size may strongly relax the constraints for high frequency domains (e.g., IR and visible) which are expected to play a larger role in photonics and quantum technology as a whole.

2. Numerical evaluation of QAE

The block diagram of the QAE is shown in Fig. 1 (top panel). A schematic of the DJ algorithm is shown in Supplement 1, Sec. A. The whole device consists of two functional subblocks, oracle and Fourier transform, respectively, brought into optical contact. Figure 1 (bottom panel) represents the schematic view of the metamaterial-based quantum emulator. The oracle block is made from a 500-µm-thick silicon (Si) substrate. It modulates the electric field profile of the incident THz light by assigning a phase shift of 0 or $\pi$ on each spatial position along $y$-axis. This phase shift is introduced by physically varying the radius of the holes array drilled in the silicon substrate, or oracle block, and by consequence also varies the effective electric permittivity $\epsilon _e$ along the y-axis thereby encoding the function $f(y)$.

 

Fig. 1. Schematic of the Quantum Algorithm Emulator. (Top panel) Block diagram of the DJ algorithm and (Bottom panel) metamaterial-based multi-layer configuration of the DJ algorithm emulator. The input THz wave of width D is modulated by the oracle block of varying hole diameters and then transformed through the Fourier block to show the output signal of either a constant function (top output) or balanced function (bottom output).

Download Full Size | PDF

The Fourier transform block is made from 127-µm-thick polyimide film with various hole sizes acting as an all-dielectric GRIN metalens as shown in Figs. 2(a) and (b). Each hole has a period of 70 µm with varying radii of 10, 15, 20, 25, and 30 µm, respectively from the center to the edge of the device. Due to the rotational symmetry imposed by the geometry of the design, the structure has an identical response for both linearly TE-polarized and TM-polarized waves. Also, note that the use of flexible substrates provides an unprecedented route to achieve frequency tunability due to modifications in the profile and the periodicity of the structures when the substrates are manipulated mechanically [22].

 

Fig. 2. (a) Structural design of the unit cell with the relevant geometrical dimensions: $p_{x}$ = $p_{y}$ = 70 µm, $t$ = 127 µm and $r$ varies 10, 15, 20, 25 and 30 µm, respectively. (b) Schematic view of the GRIN lens. The aperture size of the lens is $\sim$ 1 mm. (c) Index profile of the structure at 0.8 THz. The blue solid line is a parabolic fit to the refractive index profile. (d) Simulated electric field intensity on the focal plane at 0.8 THz. The input electric field intensity is also plotted for comparison. (e) Simulated normalized electric field distribution of the GRIN lens at 0.8 THz, for linearly TE-polarized radiation.

Download Full Size | PDF

The polyimide film is treated as a dielectric with $\varepsilon = 3.3 + i0.05$. For the GRIN lens design, we chose a radially symmetric refractive index gradient following a parabolic index profile $n(r)= n_0\text {sech}(\alpha {r})$ with $\alpha =1/r_0 + \cosh ^{-1} (n_0/n_r0)$ [Fig. 2(c)]. For this purpose, we introduced a spatial variation of the refractive index by arraying unit cells of different radius such that the refractive index gradually decreased from the center of the GRIN lens. The dimensions of each cell are modified to fit the refractive index profile on the device. The design of the final GRIN lens structure is shown in Fig. 2(b). The effective permittivity of the oracle block can be expressed as: $\epsilon _e(y)=(3\lambda _0/d_0)^{2}$; $\Delta \phi =0$, $\epsilon _e(y)=(2.5\lambda _0/d_0)^{2}$; $\Delta \phi =\pi$, where $d_0$ is the length of the oracle block and $\lambda _0$ the working wavelength.

To examine the performance of the designed metalens, we plotted in Fig. 2(d) the cross-section of the normalized intensity profiles in the plane $y$ = 0 (axial $x$–$z$ plane) around the focal point of the metalens at 0.8 THz. The input electric field intensity is also plotted for comparison. Shown in Fig. 2(e) we present the normalized local electric field distribution of the GRIN lens at 0.8 THz, for linearly TE-polarized radiation. In it, we aim to clearly demonstrate the successful realization of the focusing property of our design. When the incident light irradiates the surface of the oracle block, the phase of the transmitted wave is modulated with a factor $k_{0}n(y)d_{0}$, where $k_{0}=2\pi /\lambda _{0}$ is the vacuum wave vector, $d_{0}$ is the thickness of the oracle block and $n(y)$ is the effective refractive index at position $(y)$. The detecting function $f(y)$ is encoded into the input states by assigning a phase modulation on each spatial position $(y)$ along the input transversal direction. The refractive index of the oracle operator is designed to achieve either 0 or $\pi$ phase distribution depending on the value of the function $f$. The Fourier transform operator is used to evaluate the final results on the output signal. There are optical processes that can produce Fourier transform of field distribution, such as diffraction and Optical spatial filtering, respectively. Generally, the far field diffraction pattern is observed at infinity. By placing a lens after the diffracting aperture, the plane at infinity is imaged onto the focal plane of the lens. This explains why a lens can perform a Fourier transform.

To evaluate the performance of the device, we performed numerical simulations based on the finite-difference time-domain (FDTD) method. The length scale of the mesh was set to be less than or equal to ${\lambda _{0}}/{10}$ throughout the simulation domain, where ${\lambda _{0}}$ is the central wavelength of the incident radiation. The input and output ports are located at about 10${\lambda _{0}}$ from the device with open boundary conditions. The blue line plots in Fig. 3 show the electric field intensity at the focal plane of the GRIN lens with the detecting function being constant and balance, respectively, computed for the initial numerical analysis at 0.8 Thz. The maximum intensity at y = 0 position means the encoded function $f(y)$ is constant. However, if the center intensity is zero this indicates that the oracle subblock carries a balance function. Since the 0 and $\pi$ phase elements correspond with different effective permittivities, it makes the oracle subblock processes spatially varying impedance which results in different transmissions creating a non-symmetric intensity for the output of the electric field as shown in Fig. 3 (blue line, right panel).

 

Fig. 3. The electric field intensity at the focal plane of the GRIN lens for an intermediate working frequency of 0.8 THz, with the detecting function being constant (left panel) and balance (right panel), respectively. The blue and red lines are simulations of the output pre- (original) and post-application (optimized) of machine learning based inverse design.

Download Full Size | PDF

3. Machine learning optimization of QAE

To further engineer the structural parameters, a machine learning (ML) based optimization procedure is implemented into the design of the grin lens. ML is a powerful tool for optimization [23–28] and is expected to be an efficient complementary of electromagnetic wave numerical simulations [29–31]. Recently, one of the most popular ML techniques used in conjunction with EM field wave simulations is the inverse-design model [32–36]. Comparing to conventional methods, inverse-design traces the geometry configurations from the output performance such as resonant frequencies or S-parameters [28,37]. The overall procedure used in this study is shown in Fig. 4. First, the data for training the proposed neural network (NN) model is generated using the numerical simulation software (i.e., CST) based on the initial structure parameters as in Fig. 4(a). After which, the geometry configuration is randomly generated using MATLAB with the constraints on $r_{i}$ and h, respectively. The constraints for $r_{i}$ are $0<r_i<a/2$ and $r_{1}<r_{2}<r_{3}<r_{4}<r_{5}<r_{6}$ where $a = 70~\mu m$ is the characteristic length of the unit cell of each hole in Fig. 4(b). Additionally, the constraint of $h$ is $50~\mu m<h<300~\mu m$. To optimize the performance of the GRIN lens, we divided the holes on the surface into six groups from inner to outer ring with the radius of $r_{1}$, $r_{2}$, $r_{3}$, $r_{4}$, $r_{5}$, and $r_{6}$, respectively. In addition, the thickness of the grin lens h is also involved into the optimized procedure since the path of the wave’s propagation can also affect the performance. Additional details are provided in Supplement 1, Sec. B.

 

Fig. 4. The schematic of the proposed ML-based inverse-design of the Fourier subblock. (a) Unit cell used for generating the whole GRIN lens with varying hole diameters shown in (b). (c) Result of the balanced function which is used to apply inverse design to best fit the peaks (d) comparing the initial data of amplitude at $\pi$ and the full width at half maximum (FWHM) of $\pi$ the two peaks as the input parameters. The focused spot sizes are illustrated in the insets for the initial and optimized device.

Download Full Size | PDF

According to the Deutsch-Jozsa algorithm, the key point of the design is in fact the performance of the balanced function case since the output of all zero or all one case will absolutely be one peak. It is important that the proposed structure has an ability to distinguish the input status (i.e., constant or balance) from the output performance. Therefore, the proposed optimization procedure mainly focuses on the output spectrum [38,39]. Hence, we only need to optimize the performance of the balanced function shown in Fig. 4(c). Moreover, since the output of the half-zero half-one case is expected to have two peaks on two sides, the input for the inverse design model can be simplified as the features of the peaks. Here, we choose the amplitude $A_{pi}$ and the full width at half maximum (FWHM) $FWHM_{pi}$ of the two peaks as the input parameters, i.e.,

The data is then split into the training and validating set with the ratio of $80{\% }-20{\%}$. The training set is used to train the proposed NN model, while the validating set is used to check the performance of the model after each training round. Since the length of the input and output vectors are four and seven only, the artificial neural network (ANN) is well enough to achieve the optimal design. There are seven hidden layers in the proposed ANN model, each hidden layer has 16, 16, 32, 32, 16, 16, and 7 channels, respectively. All the layers except for the last one come with a rectified linear activation function (ReLU). The channel in the last layer is the output of the NN model representing the geometry configuration from Fig. 3 which is represented in the following equation:

After training with 100 iterations, the proposed model can precisely predict the geometry configuration from $(\bar v_{out})$. Hence, we can feed the model with the desired $(\bar v_{in})$ so that the desired optimal geometry configuration $(\bar v_{out})$ is achieved. However, if the values of $(\bar v_{in})$ exceed the performance limitation of the proposed GRIN lens, $(\bar v_{out})$ may not get the same performance as $(\bar v_{in})$. As a result, it is required to re-validate the connection between $(\bar v_{in})$ and $(\bar v_{out})$ via the numerical simulation software. After several trials, an optimal geometry configuration is achieved and validated by FDTD method. According to Fig. 3, the performance with respect to the output spectral balance is much better for the inverse-design than the initial design. The detail comparison of the initial and optimal design is represented in Table 1.

Table 1. Top: Comparison of initial numerical GRIN lens parameters to the machine Learning optimized parameters. Bottom: The performance comparison of the initial and optimal design for (i) all-zero all-one case and (ii) half-zero half one case as represented by the output wave characteristics.

View Table

From Table 1, the optimized design of the GRIN lens shows minimal increases of all radii except $r_{3}$ with changes being no more than 15%. However, the biggest change occurred with a factor of two increase of the thickness. This means the radii sizes were near optimal from the initial numerical analysis and that the amount of the computational metamaterial interacting with the incident wave is crucial in achieving the best output amplitudes and full width half maximum (FWHM) for both cases. The thickness of the GRIN lens needs to be at least one-working wavelength in such a way that the wave can better interact with the medium. The thicker the medium the better the interaction will be. However, if the GRIN lens is too thick, the focal point may fall inside the GRIN lens.

For the constant function case the amplitude has been increased by about 26 % and the FWHM has been decreased by 52%. For the balanced function, the first peak had its amplitude increased by 50% and its FWHM decreased by 26% while the second peak was already close to optimal. Knowing that the efficiency of this method is proportional to the ratio the FWHM of the incoming and outgoing waves [21] we can see from the spectral improvements that the optimization enhanced the device. For the output of either balanced or constant functions, the sharper peaks (higher amplitude and smaller FWHM) mean greater distribution of the wave intensity, stronger focus of the outgoing wave, and a better probability of handling a higher number of database inputs [4,21] as shown in bottom right of Fig. 4 and in the red plot of the right panel of Fig. 3. Cheng et al. also showed an enhancement of their balance function through increasing the number of phase elements in their oracle block.[20] Our design makes the fabrication process simpler for this improvement in output in that we only have to manufacture a thicker GRIN lens.

4. Conclusion

To summarize, we designed, with machine learning, an optimized all-dielectric metadevice as a part of a quantum algorithm emulator for simulating the DJ algorithm in the THz region. Initial structural optimization was constructed using numerical analysis based on FDTD to evaluate its initial performance. Using a ML-based optimization procedure, the original design of the structure was further engineered for an enhanced performance as shown by the two-fold increase in the thickness of the GRIN lens. The resultant optimization showed performance improvements in the amplitude and FWHM of the peaks for both the balanced and constant cases and shows promise for the emulation of quantum algorithms with current THz technologies.