The design of photonic devices is usually done through analytical modeling or variation in geometry and material parameters to obtain the required functionalities. Here, we report the use of topology optimization to obtain a bi-focal lens that concentrates the electromagnetic field at different spatial positions depending on the wavelength. Numerical inverse design is carried out to obtain the permittivity layout, satisfying this objective. The resulting device is then 3D printed using two low-loss dielectrics, and experimental field mapping at microwaves demonstrates the ability to enhance the field locally at distinct locations for two separate frequencies.

The development of photonic and electromagnetic devices with the characteristic size of the order of the wavelength has historically relied on trial and error approaches. Starting from a given physical process, a small set of key parameters is tuned to achieve an acceptable level of matching with a predefined figure of merit required by the application. This intuition based method has helped to develop a diverse and extensively used collection of designs, taking advantage of photonic resonances, dispersion engineering, waveguiding, or antenna radiation principles, enabling an increasingly finer control of light across the electromagnetic spectrum.

In the past two decades, topology optimization (TO)1 has become a widely used tool in computational nanophotonics2 and has allowed the inverse design of a broad range of devices, such as invisibility cloaks,3,4 illusion devices,5 metasurfaces,6–8 photonic crystals,9,10 and metamaterials.11,12 Broadly speaking, density based TO is an inverse design procedure that can produce highly optimized structures serving a dedicated objective. One of its main advantages is to offer unparalleled design freedom since the material distribution is updated locally (at the pixel or voxel level) inside the domain of interest. On the other hand, fabrication constraints often limit this versatility, and several auxiliary tools can be included to tackle these issues,13 for instance, for imposing minimal length scales or ensuring the connectivity of the resulting layout. Since the number of degrees of freedom is usually prohibitively high to obtain the gradient using naive finite differences, adjoint sensitivity analysis14 is an indispensable part of all inverse design algorithms. Recent advances in additive manufacturing,15 whereby a component is built up layer by layer, have allowed us to open up the design domain significantly and relax fabrication constraints. Techniques such as 3D-printing can now routinely produce components with intricate complexities, enabling the realization of optimized materials and devices with improved performances.16–19 

A lens is usually a transmissive electromagnetic device that focuses or collimates radiation. At radiofrequencies, the lens focuses the incoming radio waves onto the feed antenna in a receiving setup and, conversely, collimates the wave radiated by small feed into a beam for the transmissive mode. Microwave lenses are used for a variety of practical applications, including high gain antennas, passive beam formers, and free-space material characterization systems. The ability to control the near and far-field characteristics depending on the wave parameters, such as frequency, polarization, or incidence angle, is of paramount importance in these technologies.

We present the numerical simulation, topology optimization, manufacturing, and experimental validation of the bi-focal lens that can concentrate the electric field at two different points according to the excitation frequency. Full wave electromagnetic simulations from a two dimensional finite element model are used to compute the response of the considered structure excited by a Transverse Electric (TE) polarized line source. A gradient-based topology optimization is performed to find the layout of a dielectric device made of two low-loss dielectric materials that maximize a figure of merit for the local field enhancement. The obtained design is then 3D-printed by a fused filament fabrication (FFF) based method. Finally, spatially resolved near field measurements in a waveguide configuration are realized and demonstrate the dual focusing effect.

For numerical simulations and optimization, we assume a planar problem independent of the z-direction and consider a square design domain Ωdes of side length a = 90 mm centered at the origin of a Cartesian coordinate system, in which the permittivity distribution is parameterized by a continuous scalar density function ρ ∈ [0, 1]. The feature size needed in our design should be of the order of the operating wavelength λ: at 10 GHz, we constrain the minimum feature size to be λ/10 ≃ 3 mm. This can be easily achieved with state of the art fused filament based 3D printers. Thus, a filtered density ρ̃ is used to avoid these smaller features and pathological “chessboard” patterns (regions of alternating elements with different materials ordered in a checkerboard-like fashion) by solving the following Helmholtz equation on Ωdes20 with homogeneous Neumann conditions:

R22ρ̃+ρ̃=ρ,
(1)

with R a real positive parameter characterizing the filter radius (R = a/20 in the rest of this paper). A threshold projection is used to obtain a binary design consisting of two materials only,21 

ρ̂(ρ̃)=tanhβν+tanhβ(ρ̃ν)tanhβν+tanhβ(1ν),
(2)

where the level of projection ν = 1/2 and β is a real positive parameter and is increased during the course of the optimization. Finally, the permittivity is defined with the solid isotropic material with the penalization (SIMP) method22 as

ε(ρ̂)=(εmaxεmin)ρ̂+εmin,
(3)

where ɛmin = 3 (1 − j tan δ) and ɛmax = 6.5 (1 − j tan δ) are the permittivities of the two dielectric materials, respectively, and δ = 0.004 is the loss tangent. The gradient-based optimization is initialized with a density ρ0 and performed for 20 iterations or until convergence on the objective or design variables. This step is then repeated for n global iterations setting β = 2n, where n is an integer between 0 and 7, restarting the algorithm with the optimized density obtained at the previous step as an initial guess.

The design domain is embedded in free space and is excited by an electromagnetic line source located at rs = (−ds, 0), where ds = 260 mm. The electric field is assumed to have components only along z (TE polarization), i.e., E = Eez, and Maxwell’s equations, assuming non-magnetic materials, can be combined into the following propagation equation:

L(E,ε,k0)=2E+k02εE=δ(rs),
(4)

where k0 = 2π/λ0, λ0 is the wavelength in free space, and δ is the Dirac distribution.

Our objective is to produce a focal point at two different locations depending on the excitation frequency in the X band of the spectrum: at the top of the metalens r1 = (0, d) for a frequency f1 = 9 GHz and at the bottom r2 = (0, −d) for a frequency f2 = 10 GHz. More precisely, the problem consists in maximizing the amplitude of the electric field at two locations and can be stated as

maxρ(r)Φ(E)=E1(r1)+E2(r2),s.t.LE1(r),ερ̂(r),k1=δ(rs),LE2(r),ερ̂(r),k2=δ(rs),ε(ρ̂)=(εmaxεmin)ρ̂+εmin,ρ̂(ρ̃)=tanhβν+tanhβ(ρ̃ν)tanhβν+tanhβ(1ν),R22ρ̃(r)+ρ̃(r)=ρ(r),0ρ(r)1,
(5)

with k1 = 2πf1/c, k2 = 2πf2/c, and c the velocity of light in vacuum. We choose the focal point to be at around one-half of the wavelength at 10 GHz away from the side of the lens, i.e., d = 60 mm. The optimization process consists in maximizing the figure of merit Φ, starting with a homogeneous density ρ = 0.5. At each iteration, the density is then filtered, projected, and interpolated to obtain the current permittivity distribution, and the two electromagnetic problems (for the two frequencies) are solved. After that, the sensitivities (the gradient of the objective function with respect to the density) are calculated by the adjoint method,23 providing a direction of improvement, and the density is then updated accordingly.

All the numerical results are obtained using open source libraries with bindings for the python programming language using a custom code.24 Equations (4) and (1) are solved by the finite element method (fenics25), the sensitivity analysis is then obtained by the adjoint method (with automatic differentiation through dolfin-adjoint26), and the method of moving asymptotes (MMA27) is used for the gradient-based optimization (nlopt package28).

The optimization history is shown in Fig. 1, where we plot the field enhancement e1=E1(r1)/E10(r1) (at the top of the lens at 9 GHz) and e2=E2(r2)/E20(r2) (at the bottom of the lens at 10 GHz), where E10 and E20 correspond to the electric field in free space at frequencies f1 and f2, respectively, as a function of the iteration number. Both values increase and reach a local maximum around e1 = e2 = 3.9 at the end of the optimization process, showing the numerical convergence of the algorithm. Note that the discontinuities seen correspond to an increase in the projection parameter β in Eq. (2), allowing us to obtain a binary design from a continuous density (see insets in Fig. 1). The final permittivity layout shows well defined regions and contains no single-pixel features, with a minimum size corresponding to the filter radius used in Eq. (1). In addition, using two solid dielectric materials implies that we do not have to take into account the connectivity of each phase, which would be needed if one of the materials was air to avoid freestanding features that are to be avoided for a practical fabrication of the lens. Simulated field maps [Figs. 2(a) and 2(c)] show a field enhancement at the top for f = 9 GHz and bottom for f = 10 GHz, as required by the figure of merit. This field enhancement comes from the excitation of modes in the dielectric lens with resonant frequencies close to the target frequencies, resulting in constructive interferences adding up to increase the electric field at the corresponding prescribed focal point. This is in contrast with lenses operating by refraction in the geometrical optics approximation.

FIG. 1.

Optimization history. The insets show the filtered and projected density ρ̂ at various steps of the optimization.

FIG. 1.

Optimization history. The insets show the filtered and projected density ρ̂ at various steps of the optimization.

Close modal
FIG. 2.

Experimental validation of the bi-focal lens: simulated (a) and measured (b) electric fields at 9 GHz with focusing on top of the metalens and simulated (c) and measured (d) electric fields at 10 GHz with focusing at the bottom. A cylindrical wave emanating from a line source on the left located at rs = (−260, 0 mm) is incident on the structure. (e) 2D mapping measurement setup. (f) CAD model (top) and photograph (bottom) of the 3D printed lens.

FIG. 2.

Experimental validation of the bi-focal lens: simulated (a) and measured (b) electric fields at 9 GHz with focusing on top of the metalens and simulated (c) and measured (d) electric fields at 10 GHz with focusing at the bottom. A cylindrical wave emanating from a line source on the left located at rs = (−260, 0 mm) is incident on the structure. (e) 2D mapping measurement setup. (f) CAD model (top) and photograph (bottom) of the 3D printed lens.

Close modal

The optimized planar permittivity distribution is extruded in the z-direction by 14 mm and printed using a Raise3D Pro 2 that is a dual extrusion, fused filament fabrication (FFF) based 3D printer. The dual focus lens was designed around available ceramic loaded Acrylonitrile Butadiene Styrene (ABS) filaments produced by PREPERM; the specific filaments used in the design were “PREPERM 3D ABS300” and “PREPERM 3D ABS650” that have advertised the relative permittivity values of 3.0 and 6.5, respectively, and the loss tangent values of 0.0046 and 0.0034, respectively. Before printing the lens, the filaments were characterized, and the printing parameters were optimized by printing sample tiles (45 × 45 × 0.5 mm3) and measuring them in a 10 GHz split post dielectric resonator. This was necessary as the material flow rate; the extrusion overlap and the printing speed influence the density of the printed part (and hence effective permittivity).29 0.6 mm diameter nozzles were used on both extruders in order to reduce the chance of blockages while retaining a reasonable print definition. After optimization, the measured sample tile for the ABS300 filament had a relative permittivity of 2.94 and a loss tangent of 4.422 × 10−3, and the sample tile printed from the ABS650 filament had a relative permittivity of 6.79 and a loss tangent of 3.601 × 10−3. The two parts of the dual focus lens were then printed together in one process with the optimized printing parameters [see Fig. 2(f)].

A near-field scanning system operating at frequencies between 6 and 12 GHz is used to map the electric field [see Fig. 2(e)]. It is composed of two parallel conducting plates, each with a diameter of 1 m and spaced 16 mm apart, ensuring single mode operation for the frequency band of interest. The source is an x-band waveguide (single-mode operational frequency range: 8.2–12.4 GHz and cutoff frequency: 6.557 GHz) fixed to the rotating bottom plate and produces a cylindrical wave from its end, similar to a line source excitation.5 Microwave absorbers are added at the boundaries of the system in order to reduce scattering and reflections from the edges. Holes have been drilled along the radius of the fixed top plate with a spatial resolution of 5 mm. A monopole probe connected to the first port of a vector network analyzer (VNA) measures the amplitude and phase of the S-parameters with respect to the source waveguide connected to the other port. By varying the position of the probe along the radius and rotating the bottom plate, we obtain a 2D polar map of the near field. A computer controls and synchronizes the motor and the VNA through a LabView code.

Measured field maps for the target frequencies are shown in Fig. 2. One can clearly see the focal point at the two prescribed positions, at the top for f = 9 GHz [Fig. 2(b)] and at the bottom for f = 10 GHz [Fig. 2(d)], and prove our ability to measure spatially resolved complex electric fields with sub-wavelength resolution. In both cases, measured results are in good agreement with FEM simulations [Figs. 2(a) and 2(c)], demonstrating experimentally the frequency dependent local field enhancement. Discrepancies are attributed to a residual air gap of around 2 mm between the top of the device and the top plate of the waveguide, which is needed to be able to map the field on top of the lens with a short dipole antenna.

The field enhancement frequency spectra are displayed in Fig. 3. We identify two resonant peaks corresponding to the desired frequencies: at 9 GHz for the focal point on top of the lens and 10 GHz for the bottom point. Measured data (red and blue markers) are in good agreement with numerics for the spectral position and quality factor of the resonances, with a slightly lower field enhancement (3.4 at 9 GHz and 3.1 at 10 GHz), confirming the resonant nature of the field enhancement. These experimental results demonstrate local field control with a multi-frequency optimized all dielectric metalens.

FIG. 3.

Field enhancement as a function of frequency. Red and blue solid lines correspond to simulated results for the top and bottom focal points, respectively. Red circles and blue squares are the measured values for the top and bottom focal points, respectively. Dashed lines are a moving average for experimental results using a widow with five frequency points.

FIG. 3.

Field enhancement as a function of frequency. Red and blue solid lines correspond to simulated results for the top and bottom focal points, respectively. Red circles and blue squares are the measured values for the top and bottom focal points, respectively. Dashed lines are a moving average for experimental results using a widow with five frequency points.

Close modal

To conclude, we presented an inverse design methodology to obtain a dielectric permittivity layout acting as a frequency and spatially dependent focusing lens. Our approach allows us to find the material distribution within a design domain maximizing the electric field at predefined locations for targeted frequencies. The convergence of the method was illustrated, and the obtained structure made of two solid materials was 3D printed and its response was characterized by spatially resolved field measurements, showing an optimized near field shaping. Field enhancements of 3.4 and 3.1 at 9 and 10 GHz for focusing on the top and bottom of the metalens, respectively, were experimentally obtained. The method, combining inverse design and one-step multi-material 3D printing, allows for rapid prototyping of tailored photonic devices, including fabrication constraints. Although the bi-focusing effect has been demonstrated here at microwave frequencies, the methodology is general and suitable for the optical domain. The optimization route used here will be of interest for the design of other devices and materials with enhanced characteristics able to control the propagation and localization of electromagnetic waves in the near field, such as plasmonics, near field imaging, or second harmonic generation enhancement.

This work was funded by the Engineering and Physical Sciences Research Council (EPSRC), UK, under Grant Nos. EP/N010493/1 [SYnthesizing 3D METAmaterials for RF, microwave and THz applications (SYMETA)] and EP/R035393/1 [software defined materials for dynamic control of electromagnetic waves (ANIMATE)].

All authors declare that they have no conflicts of interest.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
M. P.
Bendsoe
and
O.
Sigmund
,
Topology Optimization: Theory, Methods, and Applications
(
Springer Science & Business Media
,
2013
).
2.
S.
Molesky
,
Z.
Lin
,
A. Y.
Piggott
,
W.
Jin
,
J.
Vucković
, and
A. W.
Rodriguez
, “
Inverse design in nanophotonics
,”
Nat. Photonics
12
,
659
670
(
2018
).
3.
J.
Andkjær
and
O.
Sigmund
, “
Topology optimized low-contrast all-dielectric optical cloak
,”
Appl. Phys. Lett.
98
,
021112
(
2011
).
4.
B.
Vial
and
Y.
Hao
, “
Topology optimized all-dielectric cloak: Design, performances and modal picture of the invisibility effect
,”
Opt. Express
23
,
23551
(
2015
).
5.
B.
Vial
,
M. M.
Torrico
, and
Y.
Hao
, “
Optimized microwave illusion device
,”
Sci. Rep.
7
,
3929
(
2017
).
6.
Z.
Lin
,
V.
Liu
,
R.
Pestourie
, and
S. G.
Johnson
, “
Topology optimization of freeform large-area metasurfaces
,”
Opt. Express
27
,
15765
15775
(
2019
).
7.
J. A.
Fan
, “
Freeform metasurface design based on topology optimization
,”
MRS Bull.
45
,
196
201
(
2020
).
8.
R.
Pestourie
,
C.
Pérez-Arancibia
,
Z.
Lin
,
W.
Shin
,
F.
Capasso
, and
S. G.
Johnson
, “
Inverse design of large-area metasurfaces
,”
Opt. Express
26
,
33732
33747
(
2018
).
9.
J. S.
Jensen
and
O.
Sigmund
, “
Systematic design of photonic crystal structures using topology optimization: Low-loss waveguide bends
,”
Appl. Phys. Lett.
84
,
2022
2024
(
2004
).
10.
J. S.
Jensen
and
O.
Sigmund
, “
Topology optimization of photonic crystal structures: A high-bandwidth low-loss T-junction waveguide
,”
J. Opt. Soc. Am. B
22
,
1191
(
2005
).
11.
A. R.
Diaz
and
O.
Sigmund
, “
A topology optimization method for design of negative permeability metamaterials
,”
Struct. Multidiscip. Optim.
41
,
163
177
(
2009
).
12.
S.
Nishi
,
T.
Yamada
,
K.
Izui
,
S.
Nishiwaki
, and
K.
Terada
, “
Isogeometric topology optimization of anisotropic metamaterials for controlling high-frequency electromagnetic wave
,”
Int. J. Numer. Methods Eng.
121
,
1218
1247
(
2020
).
13.
R. E.
Christiansen
and
O.
Sigmund
, “
Inverse design in photonics by topology optimization: Tutorial
,”
J. Opt. Soc. Am. B
38
,
496
509
(
2021
).
14.
J. S.
Jensen
and
O.
Sigmund
, “
Topology optimization for nano-photonics
,”
Laser Photonics Rev.
5
,
308
321
(
2010
).
15.
T. D.
Ngo
,
A.
Kashani
,
G.
Imbalzano
,
K. T. Q.
Nguyen
, and
D.
Hui
, “
Additive manufacturing (3D printing): A review of materials, methods, applications and challenges
,”
Composites, Part B
143
,
172
196
(
2018
).
16.
H. F.
Álvarez
,
D. A.
Cadman
,
A.
Goulas
,
M. E.
de Cos Gómez
,
D. S.
Engstrøm
,
J. C.
Vardaxoglou
, and
S.
Zhang
, “
3D conformal bandpass millimeter-wave frequency selective surface with improved fields of view
,”
Sci. Rep.
11
,
12846
(
2021
).
17.
A.
Goulas
,
S.
Zhang
,
J. R.
McGhee
,
D. A.
Cadman
,
W. G.
Whittow
,
J. C.
Vardaxoglou
, and
D. S.
Engstrøm
, “
Fused filament fabrication of functionally graded polymer composites with variable relative permittivity for microwave devices
,”
Mater. Des.
193
,
108871
(
2020
).
18.
A. W.
Powell
,
J.
Ware
,
J. G.
Beadle
,
D.
Cheadle
,
T. H.
Loh
,
A. P.
Hibbins
, and
J. R.
Sambles
, “
Strong, omnidirectional radar backscatter from subwavelength, 3D printed metacubes
,”
IET Microwaves, Antennas Propag.
14
,
1862
1868
(
2020
).
19.
H.
Giddens
and
Y.
Hao
, “
Multibeam graded dielectric lens antenna from multimaterial 3-D printing
,”
IEEE Trans. Antennas Propag.
68
,
6832
6837
(
2020
).
20.
B. S.
Lazarov
and
O.
Sigmund
, “
Filters in topology optimization based on Helmholtz-type differential equations
,”
Int. J. Numer. Methods Eng.
86
,
765
781
(
2011
).
21.
F.
Wang
,
B. S.
Lazarov
, and
O.
Sigmund
, “
On projection methods, convergence and robust formulations in topology optimization
,”
Struct. Multidiscip. Optim.
43
,
767
784
(
2010
).
22.
M.
Bendsøe
and
O.
Sigmund
, “
Material interpolation schemes in topology optimization
,”
Arch. Appl. Mech.
69
,
635
654
(
1999
).
23.
M. B.
Giles
and
N. A.
Pierce
, “
An introduction to the adjoint approach to design
,”
Flow, Turbul. Combust.
65
,
393
415
(
2000
).
24.
B.
Vial
(
2021
). “
Gyptis
,” Zenodo. .
25.
M.
Alnæs
,
J.
Blechta
,
J.
Hake
,
A.
Johansson
,
B.
Kehlet
,
A.
Logg
,
C.
Richardson
,
J.
Ring
,
M. E.
Rognes
, and
G. N.
Wells
, “
The FEniCS project version 1.5
,”
Arch. Numer. Software
3
,
923
(
2015
).
26.
S.
Mitusch
,
S.
Funke
, and
J.
Dokken
, “
dolfin-adjoint 2018.1: Automated adjoints for FEniCS and firedrake
,”
J. Open Source Software
4
,
1292
(
2019
).
27.
K.
Svanberg
, “
A class of globally convergent optimization methods based on conservative convex separable approximations
,”
SIAM J. Optim.
12
,
555
573
(
2002
).
28.
S. G.
Johnson
, “
The NLopt nonlinear-optimization package
,” http://github.com/stevengj/nlopt.
29.
A.
Goulas
,
S.
Zhang
,
D.
Cadman
,
J.
Järveläinen
,
V.
Mylläri
,
W.
Whittow
,
J. C.
Vardaxoglou
, and
D.
Engstrom
, “
The impact of 3D printing process parameters on the dielectric properties of high permittivity composites
,”
Designs
3
,
50
(
2019
).