Effective Relative Permittivity Determination of 3D Printed Artificial Dielectric Substrates Based on a Cross Unit Cell

This paper proposes closed-form analytical models for the determination of the effective relative permittivity of 3D printed artificial dielectric substrates based on a cross unit cell. The parallel plate capacitor approach is used to describe the real physical shape of the unit cell allowing to include anisotropic properties as well. A detailed comparison of the analytical models and effective medium approximations is carried out for air host material and inclusion materials with relative permittivities in the range from 2.5 to 1000 and the inclusion volume fraction from 0.1 to 1. It is observed that the proposed models predict the effective relative permittivity with much better accuracy than frequently used effective medium theorybased formulas and due to their closed-form expressions provide faster calculations than numerical methods. The proposed models were verified experimentally, achieving a very good agreement with simulations.


Introduction
In recent years, artificial dielectric structures have received a lot of attention from scientists and engineers trying to describe and measure their electromagnetic (EM) properties [1][2][3]. With the rise of additive manufacturing (AM) technologies, three-dimensional (3D) printing became more suitable for manufacturing prototypes and middle volume of units than conventional fabrication methods. Further, it also provides an extra degree of freedom in the design of desired EM properties, generates less waste of used material, and decreases manufacturing costs in general [4].
Artificial dielectrics are commonly divided into several categories, e.g., metamaterials which are usually made from heterogenous dielectric layers with a predefined arrangement of metallic inclusions [5], all-dielectric metamaterials made from high permittivity ceramics [6], or conventional dielectric materials with controlled effective permittivity which can be regulated by the inclusion and host material ratio [7]. Typical examples of those structures include perforated dielectric layers [8] or differently shaped inclusions created by AM technologies in the host material with cubical, rectangular, triangular, or hexagonal lattices at which the effective permittivity can be locally varied. Further, such structures introduce some degree of anisotropy depending on their geometry [9]. However, for an analytical description of the effective permittivity of those structures the effective medium theory based on inclusion filling factor approximations [2], [10] is usually used, which can result in an inaccurate solution in the case of anisotropic or more complexly shaped structures. In those cases, a precise description of the structures including their geometry is required. Typical examples of such problems are found in [9,11,12], where the authors investigate the anisotropy effects of various 3D printed unit cell configurations. They have found that unit cells containing symmetrical shapes like cubes, spheres or prisms can be considered as nearly isotropic. However, unit cells containing nonsymmetrical shapes like rods, hexagons, or triangles show a different degree of anisotropy depending on their mutual arrangement. The consequences are clearly observed on the 3D printed dielectric resonator antenna in [13], whose impedance matching depends on the antenna's orientation what is attributed to different dielectric substrate anisotropy. Further, by adding covering layers on an otherwise uncovered 3D printed structure, i.e., creating a dielectric substrate/composite, increases its effective permittivity and uniaxial or biaxial anisotropy [14][15][16]. Consequently, the precise characterization and modelling of such structures' material parameters is essential for accurate microwave designs. Here, precise geometric models allow an accurate determination of the effective permittivity in the main principal axes of the EM wave polarization direction, e. g. in the design of the dielectric resonator antennas [13], [17], circular polarized antennas [16], [18], the phase shifters [19], or the substrate integrated waveguide (SIW) filters with locally controlled permittivity [20].
Recently, several works focusing on analytical models of the effective relative permittivity have been published [21][22][23][24][25][26][27][28]. In [21], a capacitor network model for a bianisotropic sample based on a rectangular unit cell was presented. For its derivation, vertical and horizontal capacitor discretization schemes were used. That work states that the vertical discretization scheme is more suitable to describe the effective relative permittivity tensor of low permittivity inclusions (ε r = 2.8), however, a combination of both horizontal and vertical discretization schemes provides better estimation for higher permittivity inclusions (ε r = 9.6). In [22], an electrostatic capacitor model for the evaluation of the effective relative permittivity for small (µm-dimensions) metal dummies in a CMOS chip was developed. An artificial dielectric layer with an effective relative permittivity above 1000, which is described by a combination of a capacitor and transmission line model, was reported in [23]. Further, an electric network model describing the effective relative permittivity of cubical, elliptical, and spherical inclusions was introduced in [24]. Due to geometric limitations of the used spherical particles, only inclusion volume fractions ≤ 0.5 have been studied. However, inclusions of different shapes, i.e., ones with straight walls, can reach an inclusion volume fraction up to 1.
The base work for derivation of capacitor network models ( [21,22,24], this work) was originally presented by Patil et al. in [25], [26], where equivalent capacitor network models for spherical and cylindrical inclusions were deduced and applied even for extremely high permittivity (ε r = 3800) particles. The study [25] further claims that capacitor models of cubes with both vertical and horizontal discretization schemes provide a good prediction if the permittivity of each phase (inclusion or host medium) is isotropic. The study in [26] extends models in [25] to involve frequency dependent behavior by exploiting Debye and Cole-Cole dielectric relaxation models. To validate this statement, a study of cubic inclusions in an insulating press paper [27] showed a good data prediction. However, the statement related to discretization schemes was later disproved by Iglesias et al. in [28] on example of symmetric and isotropic cubic inclusions when the inclusion's symmetry does not ensure the discretization scheme selection independently [25]. Therefore, it can be concluded that the problem of the precise analytical models' development of given shapes based on its geometry is not so straightforward and depends on many variables as the inclusions' shape, symmetry, size, orientation, polarization or discretization scheme of the equivalent circuits. Further, the classic effective medium approximation formulas are very limited to describe those effects and might be strongly inappropriate. Note that many of the previously mentioned publications ( [1, 5, 6, 13-15, 17-21, 24, 27]) neither provide an error evaluation of the achieved results nor an analytical description of the structures' anisotropy.
In this paper, we propose a quasi-static analytical model for determining the effective relative permittivity and anisotropy of dielectric substrates, compounded of cross unit cells which are used in 3D printing technology. The models are based on vertical and horizontal equivalent parallel plate capacitor discretization schemes [25], enabling to involve their geometry and uniaxial anisotropic behavior, which is problematic to achieve by effective medium approximation (EMA) formulas based on an inclusion filling factor only. In addition, we study the discretization schemes' suitability for low and high permittivity inclusions. This knowledge can further serve for example as an input for machine learning optimization techniques [29], or to speed up design process of gradientindex-based structures such as dielectric lenses [30], which require many parametric calculations. The models' properties are compared to EMA formulas, numerical, and experimental results at 7 GHz (center frequency of WR137 waveguide).
The paper is organized as follows: Section 2 contains basic information about the cross unit cell structures and Section 3 presents parallel plate capacitor models of the cross unit cell with its equivalent circuit model containing vertical and horizontal discretization schemes and equations describing the effective relative permittivity. Section 4 presents analytical approximation formulas for the cell structures based on the effective medium theory. In Sec. 5, the numerical analysis of the cross unit cell is presented, and the results of the analytical models are compared to numerical ones. Section 6 provides experimental results and discussion on the accuracy of the model. Finally, Section 7 concludes the paper and suggests future work on this topic.

Dielectric Substrate Structures
An artificial 3D printed dielectric substrate can be modelled as a periodic structure, consisting of the predefined unit cells with specific dimensions which are covered by two lids [14]. However, the basic unit cell can be symmetrical (isotropic) or asymmetrical (anisotropic) [31]. Here, we focus on an isotropic unit cell (without lids) in which a uniaxial anisotropy can be introduced by adding covering lids. The cross unit cell with the lids can be exploited for the creation of a low profile dielectric substrate for planar microwave structures if the cell is organized in a square lattice. The selected cross unit cell with and without lids is shown in Fig. 1. Let us consider that the inclusion is the cross and lids (in Fig. 1 depicted by yellow color) and can be 3D printed by materials of a wide range of relative permittivity, e.g., 2.5-2.7 (PLA) [32], 5.6-5.8 (LMO-BF12 ceramics) [33], 9.5-10 (Alumina) [34], [35], 13.45 (Ag 2 Mo 2 O 7 ceramics) [36], or 48 (BVLMO/NMO ceramics) [37]. The air of the cell is considered as the host material (in Fig. 1 depicted by cyan color). Further, we will assume the subwavelength operating regime of the structure.

Analytical Models
The effective relative permittivity of the uniaxial structures is described by a diagonal tensor (1) [8]: The two components denoted as ε r,eff,|| and ε r,eff,⊥ represent effective relative permittivity in parallel and perpendicular direction, respectively. For the parallel case, the EM wave comes from the ±z direction (the vector of electric intensity is parallel with the z-axis) and is polarized along the x-y plane. For the perpendicular case, the EM wave comes from the ±x or ±y direction (the vector of electric intensity is parallel orthogonal to the z-axis) and is polarized along the y-z or x-z plane, respectively (see the coordinate system in Fig. 2).
For determining the effective relative permittivity, the parallel plate capacitor method [25] is exploited. The whole structure (dielectric substrate) is decomposed into a periodic array of unit cells. Each unit cell is modeled by sub-parallel plate capacitors, mutually connected seriesparallel, depending on the unit cell's geometrical configuration and discretization scheme. An example of the cross unit cell (Fig. 1) divided into sub-capacitors is shown in Fig. 2. The advantage of this approach is that it involves a real physical shape of the structure which can further consider its anisotropy.
All sub-capacitors are generally characterized by the area A, the layer thickness h, and the relative permittivity ε r (inclusion ε i ; host material ε h ) according to (2): The equivalent capacitor circuits of the structure for a vertical and horizontal discretization scheme are shown in Figs. 3 and 4, respectively. The unique discretization schemes are based on Patil et al. [25]. The strategy for the vertical discretization scheme (along z-axis) is connecting capacitors which are not interrupted by an inclusion in parallel (C lid , C 1 ) and the structure's remaining capacitors in a serial-parallel way. For the horizontal discretization scheme (along y-axis), the capacitors are firstly connected Vertical discretization along z-axis Capacity of cross sub-cells

Equation
No .
In Tab. 1, C i,j,o represents the capacitance of the i-th sub-cell, in j discretization with the orientation o in parallel ∥ or perpendicular ⊥ direction, ε 0 is the permittivity of the vacuum, ε h and ε i are the relative permittivities of the host material and the inclusion, respectively, ε r,eff represents the effective relative permittivity, l c is the length of the unit cell, d lid is the thickness of the lid, w h and w i are the particular widths of the host and the inclusion sub-cells, respectively ( Fig. 2), forming the capacitors.

Effective Medium Approximation (EMA) Formulas
The effective relative permittivity of heterogeneous materials is often approximated by various mixing models based on the Maxwell Garnett (MG) mixing rule [10]. Other, commonly used formulas include Looyenga (LO), Bruggeman (BG), Hashin-Sthrikman (H-S) [2], [10], [39] which can provide more accurate results than the M-G formula in some cases. Both the parallel and perpendicular values of the effective permittivity must fit into upper ε U eef and lower ε L eef bounds of the Wiener limit (38)-(39) defined as [8], [38]: where ε r,eff is the effective relative permittivity, ε i is the relative permittivity of the inclusion, ε h is the relative permittivity of the host material, and V i is the inclusion volume fraction.
Maxwell Garnett generalized formula for two material mixture (40) is given as [10]: Bruggeman and Looyenga approximation formulas (41)- (42) are further respectively defined as [2]: Tab. 2. Comparison of effective medium approximations used for 3D printed structures (ε i,max represents the maximum inclusion´s relative permittivity and δ r,max is the maximum relative error given in the reference; the references marked with N/A do not provide relative error evaluation).
A recent study of inclusions' distribution patterns [39] has shown that the effective relative permittivity strongly depends on the shape of the inclusion, but it was done only for two-dimensional (2D) structures using the finite difference numerical method. Some of these investigated shapes can be potentially used for 3D printed substrates, however, they require more precise characterization, especially for higher inclusion volume fractions and higher material permittivity contrast, i.e., ratio between maximum and minimum permittivity value. The comparison of effective medium approximations used for various inclusion patterns applicable for 3D printing and their deviation between theoretical and experimental (reference) results is presented in Tab. 2.
In the next section, we will present our analysis based on the parallel plate capacitors method which can potentially provide higher accuracy than effective medium approximations.

Numerical Analysis
Let´s now carry out the detailed analysis of the pro- parallel or perpendicular direction, we can calculate the value of the effective relative permittivity from the propagation constants and the free-space wavelength (43): where β is the determined phase constant and λ 0 is the freespace wavelength.
The simulation results are eigenfrequencies for the phase shift from 0° to 180° along corresponding directions. The propagation constants associated with the eigenfrequencies are then evaluated and used to calculate the effective relative permittivity of homogenized dielectrics by (43).
In the case of the E-field aligned parallel with the z-axis for ε eff,⊥ determination, we set the top and bottom boundary conditions to be perfect electric conductors (PEC) and lateral boundary conditions to be perfect magnetic conductors (PMC). On the front and rear side of the unit cell we applied periodic boundary conditions with gradually increasing values of phase difference from 0° to 180°. To determine ε eff,∥ , we rotated the structure by 90° around the x-axis. Considering maximum resolution of the 3D printer [40] and subwavelength size of the inclusion, we chose the length of the unit cell to be 0.5 mm and the thickness of the lids to be 50 µm The simulations were performed for the relative permittivities of the inclusion of 2.5, 10, 100 and 1000 (we chose the first value according to our material characterization of the Prusament PLA Jet Black [41], the second value is assumed to correspond with the PREPERM ® 3D ABS1000 filament [42], the third and fourth value is considered for method evaluation and possibly as for future composite based on BaTiO 3 ceramics [43][44]). The maximum investigated frequency corresponds to the size of the unit cell λ m /10, where λ m is the wavelength in the material of the given effective relative permittivity. Going beyond those frequencies, the whole structure will start to behave as a photonic crystal, creating bandgap where incoming EM waves cannot properly propagate [45], which is visible from a dispersion diagram shown in Fig. 6 and the obtained results are no longer valid.
The choice of discretization scheme for the given polarization direction can be additionally deduced from the electric energy density shown in Fig. 7. In the case of the low inclusion volume fraction, the vertical discretization is more suitable since the dominant effect is along the z-axis, while the horizontal discretization (along x-axis) gives inaccurate results. However, in the case of the high inclusion volume fraction, the combination of both vertical and horizontal discretization provides better accuracy because the electric field density is distributed almost uniformly across the entire volume of the structure.
The relative error evaluation is depicted in Figs. 9, 11, 13. The maximum relative error of the proposed models in vertical discretization for the relative permittivity of 2.5, 10, 100 and 1000 are 1.09%, 5.79%, 17.36% and 20.01%, respectively. However, all effective medium approximations show higher relative errors, especially for higher values of the relative permittivity, e. g., the frequently used Maxwell Garnett approximation shows the maximum errors of 6.09%, 43.3%, 90.8%, and 99.4% for permittivity of 2.5, 10, 100 and 1000, respectively. The summary of maximum relative errors for each approximation method at corresponding inclusion volume fraction (V i ) is presented in Tab. 3.
Considering achieved results, we can conclude that the proposed analytical models have better or comparable precision in the whole inclusion volume fraction range than the conventional effective medium theory formulas. It can be observed that the cross structure without lids can be precisely described by a vertical discretization scheme (Fig. 3) for relative permittivities in the range 2.5-10 with the relative error under 5%. However, the horizontal discretization scheme is not applicable in such a range of relative permittivities due to its large relative error over 130%. On the other hand, if the relative permittivity of the material is high (> 100), the vertical discretization can be combined with horizontal discretization to provide better results for higher inclusion volume fraction (> 0.75). In this case, the relative error of mean results of both vertical and horizontal discretization remains below 11%. Further, frequently used Maxwell Garnett approximation formula can easily fail to predict the correct value of effective relative permittivity even for symmetrical inclusions which differ from cubes or spheres with relative error up to 99%. Other studied effective medium approximation formulas could be exploited for higher inclusion volume fraction (> 0.5).
For the cross structure with lids (dielectric substrate) and both the parallel and perpendicular polarization directions, it can be observed that vertical and horizontal discretization schemes exhibit similar relative errors of effective relative permittivity as in the case of the structure without lids. The mean results of both vertical and horizontal discretization are then applicable for parallel polari-         zation direction and inclusion volume fraction above 0.6, while for perpendicular polarization direction, they do not show an advantage over results obtained from vertical discretization scheme. Based on our analysis, we can claim that our analytical are more accurate than effective medium approximations.
To further investigate the relative error of the effective relative permittivity depending on the relative permittivity and discretization schemes, we chose the inclusion Inclusion permittivity ε r = 2.5

Approximation
Maximum error at given (V i )  volume fraction of 0.9 which shows the highest relative error between individual discretization schemes. The results of this analysis for the cross cell structures are depicted in Fig. 14. It can be observed that in the case of the cross cell without lids, the relative error of vertical discretization scheme slightly grows up to 7.08% while the relative error of horizontal discretization scheme decreases to the value of 7.44% for the relative permittivity of 100. This trend continues even to higher relative permittivity of 1000. However, the mean value of both the vertical and horizontal schemes starts to provide more accurate results from the relative permittivity of 20. In the case of a dielectric substrate, the vertical discretization scheme is more suitable for determining the effective relative permittivity in the perpendicular polarization direction, while the mean value of both the vertical and horizontal schemes is more accurate in the parallel polarization direction.
Dielectric substrates based on a symmetrical isotropic cross unit cell exhibit some degree of anisotropy in the parallel and perpendicular polarization direction. Due to this, we have further evaluated the parallel plate capacitor method and its capability to describe this effect. The anisotropy degree can be calculated by (45) [9] and the results are shown in Fig. 15   The proposed models are capable to approximate the anisotropy effect with an acceptable accuracy, while the effective medium theory is not applicable in this case. The relative error of anisotropy prediction increases linearly with the inclusion volume fraction. This behavior can be observed for both the vertical and horizontal discretization scheme. For lower inclusion volume fractions (< 0.4), where the anisotropy degree achieves relatively large values, and thus it has a larger potential for applications requiring anisotropic behavior, e. g., circularly polarized antennas, the parallel plate capacitor method can be exploited with relative errors (calculated in the same way as (44)) under 30%. The summary of maximum dielectric anisotropy and the maximum relative error for the given inclusion volume fraction is presented in Tab. 4.

Experimental Verification
To verify the proposed models experimentally for cross cell-based substrates, the transmission/reflection method [46] was exploited at a frequency of 7 GHz. For the measurement in the WR137 waveguide, samples of the dielectric structures without lids were printed on an SLA 3D printer (FormLabs 2) with 25 m layer height. Structures with lids were printed on an FDM 3D printer (Prusa i3 MK3S) (Fig. 16) with 100 m layer height and a nozzle diameter of 250 m, as they can be printed without supporting structures because the lids serve as supporting walls.
The closed structures must be printed by FDM technology because in practice, the surplus UV resin (SLA printer) remains preserved inside the printed structure which creates cells with different inclusion volume fractions than designed. The substrate structures with inner dimensions smaller than 0.25 mm were also printed by an SLA 3D printer because of the FDM printer nozzle  diameter limitation. For printing, the size of the unit cells and lids was finally scaled up to 3 mm and 0.3 mm, respectively, due to better printability. The printed samples are depicted in Fig. 17.
The relative error of the proposed analytical models' effective relative permittivity (15), (19), (32), (37) As multiples of the unit cells dimensions do not exactly fit the waveguide's dimensions, additional numerical corrections were performed in CST Studio Suite (the frequency-domain solver). The fundamental unit cell of each configuration was repeated in x, y and z directions to fit the waveguide and the edge cell was partially truncated. After printing, some unwanted artifacts like polymerized resin in the gaps between individual unit cells, non-ideal perpendicularity, strings or small air holes on the unit cell walls were observed. The printed samples were then inserted into the waveguide and had to be tightly fitted. The resulting very small air gaps between samples and the waveguide walls could be taken as a source of possible discrepancy between simulated and measured results. For example, the simulated gap size of 0.3 mm along the wider side of the measured sample can cause an error of 1.8% (modeled in CST Studio Suite).
The measured relative permittivities of the grey UV resin "V4 FLGPGR04" and the PLA material "Prusament PLA Jet Black" are 2.9 and 2.5, respectively. The measured results of the printed samples (Fig. 17), analytical model (vertical discretization) and corrected values of the numerical model in the waveguide are depicted in Fig. 18.
The relative error between our analytical models (vertical discretization), simulations and measured results is shown in Fig. 19. For comparison, we favored the vertical discretization due to its high accuracy for low permittivity materials and omitted the comparison to the horizontal model because of its large relative error. The     anisotropy degree of cross cell based dielectric substrate and its relative error (model to simulation/measurement) is captured in Fig. 20. The overall comparison summary between analytical, numerical and measured results is presented in Tab. 5.
Obviously, the results of the measurement are in good agreement with the proposed model and waveguide corrected simulation results. Including the inaccuracy caused by non-ideally printed samples and the waveguide measurement method, the maximum relative error of the cross cell based substrate's effective relative permittivity (our Exact value is not provided; relative error is approximated from corresponding curves by (40). ** Simulated results only.
Tab. 6. Measured effective relative permittivity and its maximum relative error related to previous works. model relative to measurement) is 2.08% and 2.85% in the parallel and perpendicular direction, respectively. The maximum relative error of the cross cell based substrate's effective relative permittivity (our model relative to simulations) is 0.96% and 1.87% in the parallel and perpendicular direction, respectively.
In the case of the cross unit cell structure without lids, the maximum relative error (our model relative to measurement and simulations) is 5.63% and 1.42%, respectively, which can be caused by excessive resin in the structure gaps resulting in higher inclusion volume fraction than originally desired, which is typical for inclusion volume fraction higher than 0.6 for this unit cell configuration. Also, the estimation of dielectric anisotropy is verified with a maximum relative error (our model relative to measurement and simulations) of 25.4% and 52%, respectively. The maximum measured anisotropy achieved 11.9%. The comparison of measured effective relative permittivity and its maximum relative error with previous works is presented in Tab. 6.
We conclude that our proposed models exhibit very good maximum relative error considering a wide range of effective permittivities which varies in simulations from 2.5 to 1000. However, because of the materials' availability, the maximum measured permittivity was 2.9.

Conclusions
In this paper, we proposed a new model which solves a problem of precise and fast effective permittivity determination in artificial, uniaxial dielectrics in 3D printing technology which can be easily parametrized. A detailed study of a parallel plate capacitor method to model structures based on a cross unit cell was presented. We verified different discretization schemes and proved that vertical discretization is more accurate to describe effective relative permittivity for our targeted low (2.5) to middle (10) permittivity materials for inclusion volume fractions between 0 and 1. In this case, we achieve a maximum relative error of 5.8%, which outperforms frequently used effective medium theory approximations (Maxwell Garnett, Bruggeman, Looyenga). A combination of mean values obtained by vertical and horizontal discretization can be generally more accurate for cross cell based structures with high permittivity (> 100) materials for inclusion volume fractions above 0.6, or more specifically for high inclusion volume fractions (0.9) and relative permittivities above 12. The proposed model's performance was further verified by means of describing anisotropy behavior of the structures, showing good agreement for highly anisotropic cases with inclusion volume fractions below 0.4. Further, good agreement between the proposed analytical models and experimental results was achieved. The proposed analytical equations can be used for the design of artificial dielectric substrates based on the cross unit cell. Future investigations will be focused on models' accuracy limitations associated with higher inclusion fractions in terms of effective relative permittivity and dielectric anisotropy and further extension of the given models to precisely model dielectric losses.