Additive Manufacturing of Linear Continuous Permittivity Profiles and Their Application to Cylindrical Dielectric Resonator Antennas

The utilization of additive manufacturing (AM) to engineer the permittivity profile of dielectric resonator antennas (DRAs) is considered. For the first time, the capabilities of AM are exploited to create continuously swept permittivity profiles and applied to cylindrical DRAs. The spatial variant lattices (SVL) synthesis algorithm is implemented to create the desired permittivity profiles from a single material, and resulting geometries are manufactured using a high-permittivity material in a fused deposition modeling AM process. Three individual antennas for global navigation satellite system bands are designed and manufactured, two inhomogeneous DRAs with continuous permittivity profiles along the radial and vertical axis, and one homogeneous DRA for comparison. The manufactured antennas are characterized by impedance, realized gain, and axial ratio. Experimental results agree well with simulations and show increased impedance-, gain-, and axial-ratio bandwidths for both inhomogeneous antennas compared to the homogeneous one.


I. INTRODUCTION
D IELECTRIC resonator antennas (DRAs) are viable alternatives to traditionally employed radiating elements such as wire and patch antennas thanks to their high radiation efficiency due to the absence of conduction losses [1]. Depending on the design approach DRAs provide appealing gain, bandwidth, and polarization characteristics due to their versatility in shape, dielectric characteristics, and feeding mechanism [2], [3]. Furthermore, DRAs offer great potential for the miniaturization of antenna elements as their size is inversely proportional to the square root of the resonating body's relative permittivity. Thorough reviews of DRA technology highlighting historical developments and recent advancements are presented in [4], [5]. Multi-permittivity and inhomogeneous variations of the dielectric resonator are thought to be a proper way to achieve wideband operation [6], [7], [8], multiband behavior [9] or increase the spurious-free window [10]. These inhomogeneous DRAs are usually manufactured by either milling perforations into the resonator or stacking (gluing) different materials together. Alternatives to these methods are desired since the milling process is limited concerning complex internal features, and the utilization of glues is cumbersome and introduces an additional possible uncertainty source. Accompanying the developments in DRA technology, advances in additive manufacturing (AM) have exploded during the last decade. This is especially interesting concerning DRAs since AM allows the manufacturing complex dielectric structures with relatively low-cost machines and processes. Several works successfully utilized AM of dielectrics to create bulk homogeneous DRAs with various materials and shapes [11], [12], [13]. Furthermore, the AM capabilities have been exploited to create complex periodic distributions of high-permittivity material to engineer the permittivity of a DRA [14]. Additionally, multi-permittivity (inhomogeneous) DRAs have been created for bandwidth improvement [15], [16], [17], [18]. Additive manufacturing of multi-permittivity devices has also been successfully employed to create flat lenses [19]. However, so far, additively manufactured devices with inhomogeneous permittivity have only been manufactured with discrete sections of a specific permittivity. The potential of continuously swept permittivity distributions in the device has yet to be explored in literature. This work aims to address this literature gap by utilizing spatial variant lattices (SVL), as introduced in [20], to create periodically structured geometries with continuously varied volumetric fill-fraction. The goal is to show the design, simulation, and manufacturing of swept permittivity profiles and investigate their effect on the performance of cylindrical DRAs. For this purpose, linearly swept profiles with increasing permittivity along the radial and vertical axis are studied. Earlier studies on DRA with piece-wise constant permittivity profiles inspire this choice for the selected permittivity profiles, and this work focuses instead on demonstrating the feasibility of precisely manufacturing a continuous swept profile, thus enabling future research about an ad-hoc-optimization of a DRAs permittivity profile. The antennas discussed in this work are intended for applications with global navigation satellite systems (GNSS) and optimized for operation at the Galileo E1 band center frequency of 1.575 GHz. However, the proposed approach for creating continuously swept permittivity profiles can be scaled to any desired frequency band, given that the manufacturing process can produce structures sufficiently small with respect to the wavelength. The rest of the paper is structured as follows. Section II explains the DRA concept, the selected permittivity profiles, and the utilized face-centered cubic (FCC) unit cell design based on spatial harmonics. Furthermore, the approach to varying the effective permittivity throughout the DRAs resonating body is explained, and an example of a lattice with a swept volumetric infill is presented. Section III presents the optimization results to find the appropriate dimensions for DRAs with linearly swept effective permittivity profiles along their radial and vertical axis for operation at the E1 band (1.575 GHz). Insight into the manufacturing process and measurement results are presented in Section IV. Furthermore, a discussion and a conclusion are given in Sections V and VI, respectively.

A. ANTENNA GEOMETRY AND PERMITTIVITY RANGE
The individual cylindrical DRAs with radius a, height h on a cylindrical ground plane with radius a GP , as depicted in Figure 1, are fed via two 50 coaxial probes of length l feed with an offset of a feed from the DRA center. The probes are excited with 90 • phase offset (Port 1: 0 • ; Port2: −90 • ) for right-hand circular polarization (RHCP). The antennas are manufactured via a fused deposition modeling (FDM) process utilizing the ABS1500 filament from Avient (former Preperm) with a nominal relative permittivity of ε r,n = 15. However, the relative permittivity of printed substrates heavily depends on the FDM process parameters like layer height, line width, extrusion multiplier, etc., as demonstrated in [21]. Therefore, the effect of the later utilized printing parameters on the relative permittivity of a printed bulk ABS1500 sample was investigated with a material characterization setup based on rectangular waveguides as described in [22]. The experiment showed a reduction in relative permittivity from the nominal value to ε r,max = 13.7, which is further considered to be the maximum achievable permittivity in our design. The swept permittivity profiles are generated via the method introduced in [20], based on an FCC unit cell generated via spatial harmonics as described in [24]. The effective permittivity of the unit cell can be adjusted with the threshold values th, which controls the volumetric infill, as further explained below.

B. UNIT CELL DESIGN
In contrast to other works where unit cells are parametrically defined [14], [23], this work takes a different approach and defines unit cells via spatial harmonics along the reciprocal lattice vectors of the desired crystal geometry as described in [24]. Face-centered cubic unit cell symmetry was chosen in this work as it allowed for easier manufacturing. The approach via spatial harmonics allows the straightforward application of the algorithm discussed in [20], which is later used to introduce the swept permittivity profile. The unit cell is defined with its lattice constant a and primitive lattice vectors of the FCC symmetry t 1 = a 2 [0 1 1] T , t 2 = a 2 [1 0 1] T and t 3 = a 2 [1 1 0] T which point from the current lattice node to its three nearest neighbors, as depicted in Figure 2. For each set of two primitive lattice vectors, a perpendicular reciprocal lattice vector can be found via (1) where the magnitude describes the periodicity of the lattice. The lattice geometry in the respective Bravais lattice is based on a superposition of N spatial harmonics with individual grating vectors where A is a complex valued function in terms of spatial coordinates x, y, and z. By choosing the individual grating vectors, g i such that they are oriented along the reciprocal lattice vectors T i as defined in (3), and their periodicity is an integer multiple of the lattice periodicity, the superposition of spatial harmonics will obtain the respective crystal symmetry and periodicity. Before generating the unit cell geometry a min-max feature scaling is applied on the real part of the superimposed spatial harmonics A as an intermediate step withÃ now representing a real function in terms of spatial coordinates x, y, and z with values between 0 and 1. The overall unit cell UC is computed, for some threshold th ∈ [0 1], via where th represents an the value of an isosurface throughout the unit cell. Regions whereÃ ≥ are further considered to consist of material 1 (in our case the 3D printed substrate) and regions whereÃ < th are considered to consist of material 2 (in our case free space). Figure 3 shows individual unit cells with FCC symmetry with different thresholds created with N = 3 equal amplitude spatial harmonics and grating vectors equal to the FCC reciprocal lattice vectors. The benefit of utilizing spatial harmonics to create the unit cell geometry is that one can easily vary the fill fraction δ V by adjusting the inverse proportionally related threshold value th. Furthermore, the utilization of spatial harmonics allows the straightforward and computationally efficient application of the algorithm introduced in [20] to create swept permittivity profiles, as is further discussed below in Section II-E.

C. EFFECTIVE PERMITTIVITY
The relationship between the unit cell threshold th and its effective permittivity is not trivial and such relationships are usually studied with effective media theories. The Maxwell Garnett Approximation (MGA) [25] has been especially prominent for modeling the effective permittivity of additive manufactured dielectrics. However, several studies have shown that the MGA lacks accuracy when applied to high volumetric infill fractions and moderate dielectric contrasts [26]. The plane wave expansion method (PWEM) [27] is a suitable alternative to the MGA to compute the effective permittivity of additive manufactured dielectric unit cells [22], [28]. The method discretizes Maxwells equations with periodic boundary conditions on a plane wave basis. Given the permittivity distribution of the unit cell ε r,uc (r) and a Bloch wave vector β the PWEM computes the modes (eigenvectors) that fit the phase boundary condition imposed by the chosen Bloch vector β and their respective wavenumbers k 0 (eigenvalues). A complete analysis of the crystal with this method requires a dense sampling of Bloch vectors in the Brillouin zone of the crystal (Figure 2), which is computationally demanding. However, we are only concerned with the effective permittivity in the long wavelength limit, which corresponds to low-magnitude Bloch wave vectors (Bloch vectors around the key point of symmetry). We can extract the effective refractive index n eff of the fundamental mode via which for purely dielectric media (μ r = 1), is the square root of the effective permittivity ε r,eff . This method is now used to compute the effective permittivity of FCC unit cells made from a material with permittivity ε r,max with different thresholds th as depicted in Figure 3. The effective permittivity as a function of the threshold th is depicted in Figure 4. The lattice is connected to its neighboring unit cell up to a maximum threshold of th max = 0.625, and for convenient manufacturing, a maximum threshold of th max = 0.6 was chosen which corresponds to a minimum manufacturable effective relative permittivity of ε r,min = 3, which is further considered to be the minimum achievable permittivity in our design.

D. PERMITTIVITY PROFILES
Three individual cylindrical DRA versions with different permittivity profiles ε r (r), based on the geometry depicted in Figure 1 are considered in this work. A homogeneous DRA is manufactured with the maximum achievable relative permittivity while two inhomogeneous DRAs are engineered to obtain a linearly swept permittivity profile from ε r,min to ε r,max . The individual DRAs will further be referred to as: • DRA-H -Homogeneous DRA with ε r (r) = ε r,max for comparison purposes; • DRA-IHZ -Inhomogeneous DRA with a linearly swept permittivity profile along the vertical (z)-axis, ε r (r) = k z z + d; • DRA-IHR -Inhomogeneous DRA with a linearly swept permittivity profile along the radial (r)-axis, ε r (r) = k r r + d = k r x 2 + y 2 + d; with The individual permittivity profiles of the three considered DRA configurations are depicted in Figure 5.

E. SPATIAL VARIED LATTICES
The idea of utilizing a spatial variation of dielectric infill to introduce inhomogeneous effective permittivity for design purposes is not new. A common planar design approach is utilizing space-filling curves with modulated trace thickness [29] to achieve spatial-dependent material responses. However, space-filling curves are challenging to adapt to 3D geometries. 3D spatial variation has been introduced with voxelated approaches where parametrically defined unit cells are treated like individual voxels, which have a practical permittivity value assigned that is engineered by changing the geometric parameters of the unit cell. Another popular method is to create sections in a lattice that will be printed with different materials or infill fractions, as shown in [19].
Here we utilized the less known but powerful algorithm to synthesize spatially variant lattices (SVL) as introduced in [20]. The details of the SVL algorithm are not repeated here. Roughly explained, the approach decomposes a given unit cell into spatial harmonics and applies spatial variations like rotations or changes in the lattice period onto these individual harmonics. Later the whole lattice is assembled by the superposition of the individual harmonics, which results in a continuous and smooth structure with minimum deformation to individual unit cells. The SVL algorithm, although computationally expensive, offers a powerful tool to influence dielectric crystal lattices in numerous ways. Due to the computationally demanding nature of the SVL algorithm we utilize the spatial harmonics approach to design a unit cell rather than decomposing an arbitrary unit-cell geometry in a large number of spatial harmonics. The spatial harmonics design approach enables us to create a fully 3D periodic unit cell via only N = 3 spatial harmonics resulting in a manageable computation time when applying the SVL algorithm. Figure 6 depicts a cuboid consisting of 1×1×10 unit cells, where the SVL algorithm is employed to continuously sweep the fill fraction of the FCC unit cell along the z-axis.

III. DESIGN AND SIMULATION
While the homogeneous DRA dimensions can be found quickly utilizing the design equations in [30], no such equations are available for inhomogeneous DRAs. Therefore, the dimensions of the two inhomogeneous antenna variations are found via optimization in the commercial EM solver package Ansys high-frequency structure simulation software (HFSS). Since HFSS does not allow the definition of inhomogeneous material, the model for the inhomogeneous DRAs is built with distinct layers or shells of homogenous material to approximate the desired permittivity profile in the resonating body, as depicted in Figure 5. While the range of the permittivity profile is fixed form ε r ∈ [3 13.7] as explained above, the model geometry is subject to an optimization problem to match the antenna to the center frequency of the E1 GNSS band (1.575 GHz). The optimization is set up to minimize the reflection S 11 and transmission S 21 at the E1 center frequency by changing the geometric parameters a DRA , h DRA , a feed and ł feed . It has been observed that the

TABLE 1. Geometric parameters of DRA antennas and permittivity function εr(r)
after optimization with the goal obtain a minimum in the reflection coefficient at the design frequency 1.575 GHz.
results obtained from the simulations sufficiently converged around 20 shells and layers. Therefore, the optimization routine was carried out with 20 individual shells/layers, and for the final simulations, this number was increased to 50. While Table 1 lists the optimized geometry parameters for each of the three antenna variations with their respective permittivity function ε r (r), the resulting reflection coefficients, realized gain, and the axial ratio at θ = 0 • are plotted as a function of frequency in Figures 7 and 8, respectively. In Figure 7 one can observe that the reflection coefficients show respective minima approximately at the E1 center frequency, as was the goal of the optimization. Furthermore, both inhomogeneous DRAs extend the impedance bandwidth with respect to the homogeneous one. An interesting observation is that the coupling between the two feed ports is reduced for the inhomogeneous DRAs at the design frequency, which implies improved realized gain and axial ratio values, which is confirmed by the plots in Figure 8.

IV. MANUFACTURING AND MEASUREMENTS
The geometric mesh of the continuously varied lattice for the inhomogeneous DRAs is obtained directly from the

FIGURE 8. Simulated RHCP (solid) and LHCP (dashed) realized gain (top) and axial-ratio (bottom) at θ = 0 • as a function of frequency for the homogeneous (blue), radial inhomogeneous (red), and vertical inhomogeneous (yellow) cylindrical dielectric resonator antennas.
SVL algorithm and saved as a standard tessellation language (STL) file. Individual STL files are prepared for printing via the PrusaSlicer v2.5.0 which translates the antenna geometries into G-code for the respective manufacturing machine. Models are printed with 100% infill, employing a heavily customized 3D printer based on the E3D tool changer with the ABS1500 filament from Preperm. We utilize a 0.25mm nozzle at an extrusion temperature of 250 • Celsius with a heated bed at 110 • Celsius and a layer height of 0.15mm. 1 The finished antennas with installed coaxial probes and mounted on individual aluminum ground planes are depicted in Figure 9. Scattering parameters between the two feed-ports were measured with a vector network analyzer between 1 and 2 GHz and are compared to the simulation results in Figure 11. For the measurement of the DRAs radiation patterns, the individual antennas are fed via a 90 • hybrid (MiniCircuits ZX10Q-2-19+) for RHC polarization and placed in an MVG Starlab near-field chamber as shown in Figure 10. Measured realized gain in RHC and LHC polarization values at θ = 0 • , as well as the axial ratio at θ = 0 • , are plotted as a function of frequency in Figure 12 and compared to simulated results. The simulated and measured realized gain and axial ratio patterns at the design frequency are compared for the DRA-H, DRA-IHR, and DRA-IHZ in Figures 13, 14 and 15 respectively. Simulated and measured numerical values of impedance bandwidth, axial ratio bandwidth, and realized gain are compared for each of the three antennae in Table 2. From Figure 11 one is able to observe some discrepancies between simulation and measurement. First of all, the homogeneous antenna obtains a broader impedance bandwidth 1. The printing parameters are the same used for manufacturing the bulk sample which was measured to obtain a permittivity of ε r,max = 13.7, as explained in Section II-A.  than predicted from simulations. Furthermore, the reflection coefficients for the DRA-IHR and DRA-IHZ antennae do not obtain their respective minima at the design frequency but are shifted to higher frequencies. Although there are some shifts in the center frequency, the predicted effect of improved isolation between the feed ports is confirmed for both the DRA-IHR and DRA-IHZ antennas.
The measured RHC realized gains of the DRA-H, DRA-IHR, and DRA-IHZ antennas in Figure 12 show slightly reduced values compared to the simulated ones, which can be explained via the insertion loss of the utilized 90 • hybrid that was employed to achieve RHC polarized radiation. Additionally, the realized gain of the DRA-IHR antenna shows a more significant reduction in realized Gain between 1.35 and 1.55 GHz compared to simulation which most probably due ratio curves is shifted to higher frequencies similar to the shift observed in the feed coupling |S 21 | in Figure 11.
The realized gain and axial ratio patterns at 1.575 GHz depicted in Figures 13, 14, and 15 are not perfectly symmetric as the azimuth region ϕ ∈ [0 • 90 • ] experiences higher RHC realized gain and improved axial ratio values. This effect is due to the position of the two coaxial feed probes at ϕ = 0 • and ϕ = 90 • as shown in the schematic of the antenna in Figure 1. This distortion can be improved by further increasing the number of feed points of the individual antennas. Overall the shape of the realized gain and axial ratio patterns predicted by simulations is adequately reproduced in the measurement. The inhomogeneous DRAs show significantly improved axial ratio values, not only at θ = 0 • but also over a larger angular area compared to the homogeneous DRA due to the reduced feed coupling in the inhomogeneous resonator bodies.

V. DISCUSSION
Generally, the chosen simulation approaches deliver a reasonable agreement between simulations and measurements. This fact is essential to highlight due to two reasons. First, the modeling of the relationship between the volumetric fillfraction and effective permittivity of the unit cell is done with the plane wave expansion method, which assumes a plane wave propagating through the infinite lattice, while the DRA feeding can be considered to be more of a pointlike source. The point-like source excites a broad spectrum of homogeneous and inhomogeneous spectral harmonics of the electromagnetic field. Nevertheless, the results are in reasonable agreement confirming the validity of the PWEM modeling approach. Another reason is that the numerical simulation model of the DRA is not considered with the final geometry of the DRA but with individual layers and shells of homogeneous dielectrics. The layered/shelled approach was chosen because the spatial harmonic design approach cannot be replicated in HFSS. However, still, a high number of layers/shells is necessary for the model to adequately represent the intended permittivity profile, leading to a rather complex mesh and significantly increased simulation times. The observed improvement in the operational bandwidth of DRAs with heterogeneous dielectric can be attributed to a decreased coupling between the two coaxial feeds of the DRAs. The physical reason for the mitigated feed-coupling seems to be a decreased dielectric permittivity in the DRAs resonating body between both feeding points. One drawback of the inhomogeneous DRAs is their increased complexity in manufacturing and, for the DRA-IHZ only, the increased volume compared to the homogeneous version. Another potential issue is that this approach is challenging to scale to higher frequencies. The unit-cell size needs to be small compared to the design wavelength for the resulting volume to act as an effective medium but at the same time big enough to be compatible with the employed manufacturing process.

VI. CONCLUSION
This work accomplishes two goals. For the first time in literature, it describes the design and additive manufacturing of an artificial dielectric volume with a continuous linear variation of its relative permittivity and applies such permittivity profiles to improve the performance of probe-fed dielectric resonator antennas. This is achieved by utilizing a dielectric unit cell with face-centered cubic symmetry that is created via the superposition of spatial harmonics. The individual spatial harmonics are used in the powerful SVL algorithm to synthesize a volume with continuously changing volumetric fill-fraction that is then manufactured with a single highpermittivity material in a fused deposition modeling process. Two inhomogeneous DRAs with linear varied permittivity profiles along their radial and vertical axis, and one homogeneous DRA for comparison. The choice of the manufactured permittivity profiles in this work was driven by simplicity in their description rather than optimization constraints. Nevertheless, both inhomogeneous DRAs show mitigated coupling between the feed ports leading to improved realized gain and axial ratio values at the design frequency. These benefits are demonstrated with both numerical simulations and measurements, validating the design and simulation approach proposed. The presented approach to continuously linearly sweep the permittivity profile of dielectric volumes is applicable to DRAs. However, the same technique can be applied to design DRAs with inverse linear or even nonlinear profiles, which shall be addressed in future research. Furthermore, this approach can be applied in numerous other ways, for example, to engineer flat lenses or to introduce tapers in dielectric waveguides without changing their geometry.