Complex modes and artificial magnetism in three-dimensional periodic arrays of titanium dioxide microspheres at millimeter waves

We characterize the modes with real and complex wavenumbers for both longitudinal and transverse polarization states (with respect to the mode traveling direction) in three dimensional (3D) periodic arrays of titanium dioxide ( TiO 2 ) microspheres in the frequency range between 250 GHz and 350 GHz. Modal results are computed by using a single magnetic dipole approximation (SDA) and an SDA model with correction (SDA-WC) that assumes the array to be embedded in a host with an effective permittivity computed through Maxwell Garnett formulas. Moreover, for the transverse polarization case, modal wavenumbers are computed also by fitting the full-wave simulation magnetic field (one point per unit cell) in a finite thickness structure, and their agreement and disagreement are discussed. The longitudinal polarization is not affected by the artificial correction introduced in the SDA-WC; in the transverse polarization case, instead, the correction is needed to obtain results in better agreement with the full-wave data fit. In the observed frequency range, there are one longitudinal mode and two transverse modes, one forward and one backward, where the forward one is “ dominant ” (i.e., it contributes mostly to the field in the array). Therefore, in the case of transverse polarization, we describe the composite material in terms of homogenized refractive index and relative permeability, comparing results from (i) modal analysis (with and without correction), (ii) Maxwell Garnett formulas, and (iii) Nicolson – Ross – Weir retrieval method from scattering parameters of finite thickness structures. The agreement among the different methods justifies the performed homo- genization procedure in the case of transverse polarization. We show that artificial magnetism is generated from a nonmagnetic composite material. © 2012 Optical Society of America OCIS codes: 160.3918, 160.1245,


INTRODUCTION
The lack of strong magnetism in natural materials has motivated in recent years the use of metamaterials to generate artificial magnetism (i.e., the relative permeability tensor is different from the unity tensor) from nonmagnetic constitutive materials, especially at high frequencies where natural magnetism disappears.
Several metamaterial configurations have indeed been proposed to overcome this natural limitation throughout the frequency spectrum. Artificial magnetism has been generated through the split ring resonator (SRR) structure [1], initially introduced at microwave frequencies, and then extended to infrared frequencies by scaling the dimensions of the SRR [2,3] [further miniaturization would not be effective because of the growth of the kinetic (internal) inductance [3]]. Another way to generate artificial magnetism is through pair-based metamaterials, whose working principle is based on the excitation of an antisymmetric resonance associated to an equivalent current loop, such as staples [4], strips [5,6], dogbones [7][8][9][10], and metallic nanospheres/nanoshells [10][11][12], based on the original design with double bars in [13]. Alternatively, one can choose to use magneto-dielectric particles (with high values of permittivity and/or permeability) embedded in a dielectric matrix [14][15][16], or two sets of particles with different sizes and/or materials (one designed to resonate at the electric resonant mode, the other one designed to resonate at the magnetic resonant mode), to produce negative permeability and/or an isotropic double negative material [17][18][19][20][21][22][23][24][25][26][27][28]. Similar results may be achieved by using arrays made of nanoshell particles [25,[29][30]. Also, the packing of plasmonic nanoparticles in an engineered fashion to create nanoclusters [31][32][33] allows for the generation of artificial magnetism.
A comprehensive way to understand and classify collective resonances in composite materials is by modal analyses of arrays periodic in three dimensions [34][35][36][37][38][39][40][41]. The approach described in the present paper allows for the tracking and especially for the characterization of the evolution of modes varying frequency. The numerical procedure adopted in this paper for evaluating the complex zeros of the dispersion relation uses the Ewald representation [34,42] for the dyadic periodic Green's function (GF) to represent the field in threedimensional (3D) periodic arrays, which has been adapted from [41] to the magnetic dipole case. The Ewald representation, besides providing analytic continuation to the complex wavenumber space, results in two series with Gaussian convergence where only a handful of terms are needed.
In this paper, we choose to use titanium dioxide (TiO 2 ) as the constituent material for particles on several grounds: first, it has useful properties in the 0.1-1 THz band among real materials, combining very high real permittivity values (close to 100) with low losses (loss tangent of a few percents) [43]; second, for practical reasons, because it is low-cost, nontoxic, easily chemically processable, and widely available, and therefore makes an excellent candidate for fabricating reallife metamaterials. In particular, we aim at showing that a composite material made of nonmagnetic TiO 2 microspheres (i.e., with unity relative permeability) allows for the generation of artificial magnetism at millimeter wave frequencies. This means that, when homogenized, the metamaterial properties are described by an isotropic effective permeability μ eff ≠ 1. As will be shown in Section 4, dealing with a metamaterial with finite thickness, despite the presence of two modes for transverse polarization, only one is "dominant" in the sense that it contributes mostly to the field in the array. Therefore, in case of transverse polarization, we describe the composite material in terms of homogenized effective refractive index and permeability. These parameters are obtained by using four methods as detailed in Section 5, and the agreement among them confirms the validity of the homogeneous treatment of the composite material. We extend and provide the analytical formulation based on GFs for magnetic dipoles (analogously to what was done for electric dipoles in [41]) and propose a corrected method based on single dipole approximation (SDA) to obtain results in better agreement with full-wave ones [SDA with correction (SDA-WC); see Subsection 3.A]. Furthermore, the field fitting procedure introduced in [41] for a 3D array of plasmonic nanospheres is improved by considering here a superposition of two sets of traveling waves with different wavenumber leading to the reconstruction of the dispersion diagram with two modes, one forward and one backward (see Section 4).
The structure of the paper is as follows. We discuss in Section 2 how to model the scattering arising from each microsphere under magnetic dipole approximation, and we summarize the rigorous representation of the field in 3D periodic arrays using the dyadic GF in the particular case of magnetic dipoles. The detailed computation of the dispersion diagrams for longitudinal and transverse polarization states (with respect to the mode traveling direction) for a set of array parameters is then reported in Section 3. The magnetic field calculated via a full-wave simulation based on the finite element method [High Frequency Structure Simulator (HFSS) by Ansys Inc.] on systems made of 11 layers of arrayed TiO 2 microspheres, stacked in the direction of propagation of the normally incident plane wave illumination, is fitted via a Bloch mode analysis in Section 4, justifiying the treatment of the composite material as a homogenized material, with effective refractive index and permeability. The retrieval methods for the computation of these parameters are summarized in Section 5. Last, in Section 6, these methods are adopted to compute the effective refractive index and permeability, and their agreement and disagreement are discussed. Also, these parameters are compared against parameterization of the periodicity of the 3D array. In addition, the evaluation of the Ewald representation for the dyadic GF for 3D periodic arrays of magnetic dipoles is summarized in Appendix A.

FORMULATION
The monochromatic time harmonic convention, exp−iωt, is assumed here and throughout the paper, and is therefore suppressed hereafter. In the following equations, bold letters refer to vector quantities, a caret "^" on top of a bold letter refers to unit vector quantities, and a bar under a bold letter refers to dyadic quantities.

A. Microsphere Modeling
Limiting the scattering from each microsphere to dipole-like scattering, the induced electric and magnetic dipole moments of a spherical particle are where α ee and α mm are the isotropic electric and magnetic polarizabilities of the microsphere, and E loc and H loc are the local electric and magnetic fields acting on the microsphere. According to Mie theory, the electric (a 1 ) and magnetic (b 1 ) scattering coefficients for a spherical particle are [44] a 1 mψ 1 mkrψ 0 and the electric and magnetic polarizabilities are [36,39] α ee 6πiε 0 ε h k 3 a 1 ; α mm 6πi where r is the microsphere radius, ψ 1 ρ ρj 1 ρ sin ρ∕ρ − cos ρ and ξ 1 ρ ρh 1 1 ρ −i∕ρ − 1e iρ are the Riccati-Bessel functions [45] [a prime in Eq. (2) refers to the first derivative of the function with respect to its argument], k ω ε h p ∕c 0 k 0 ε h p is the host medium wavenumber, where k 0 ω∕c 0 denotes the free space wavenumber and c 0 is the speed of light in free space, and ε h and ε m are the relative permittivity of the host medium and of the TiO 2 microspheres, respectively, so that m ε m ∕ε h p is the relative refractive index of the inclusions.
In what follows, we consider the frequency range between 200 GHz and 500 GHz, where, referring to the experimental data in [43], we assume the permittivity of titanium dioxide to be ε m 94 i2.35 (measured at 500 GHz). The value Imε m 2.35 represents a worst case scenario in the analyzed frequency range 200 GHz and 500 GHz, assumed to take into account extra losses that often arise in microparticle synthesis (due to lower cristallinity, nonhomogeneity, roughness, etc.). Finally we take ε h 1 for the host, and therefore k k 0 , and r 52 μm for the sphere radius. Under these assumptions, to establish which of the two dipolar scattering effects (i.e., electric or magnetic) is stronger, we observe the behavior of the magnitude of the Mie scattering coefficients in Eq. (2) versus frequency, reported in Fig. 1(b). It can be observed that the magnitude of b 1 is stronger than the magnitude of a 1 in the frequency range between 250 GHz and 350 GHz (with the peak at around 294.5 GHz), thus in this range each TiO 2 microsphere experiences a "magnetic-like resonance" and can be modeled as a single magnetic dipole (i.e., electric effects are negligible) by using an SDA method [39,41,44].

B. SDA for 3D Periodic Arrays of Microspheres Modeled as Magnetic Dipoles
Consider a 3D periodic array of microspheres as in Fig. 1(a), immersed in a homogeneous background, with relative permittivity ε h , for which each microsphere is placed at position r n r 0 d n , where n n 1 ; n 2 ; n 3 is a triple index, and d n n 1 ax n 2 bŷ n 3 cẑ, with n 1 ; n 2 ; n 3 0; 1; 2; …, r 0 x 0x y 0ŷ z 0ẑ , and a, b, and c are the periodicities along the x-, y-, and z-directions, respectively [39,41].
According to Subsection 2.A, each TiO 2 microsphere can be modeled as a single magnetic dipole using the SDA in the frequency range between 250 GHz and 350 GHz, i.e., ja 1 j ≪ jb 1 j (dimensionless), as in Fig. 1(b), which in turn leads to α ee ∕ε h ε 0 ≪ α mm (expressed in m 3 ). Accordingly, the SDA is a good approximation when particle dimensions are much smaller than the wavelength, and when the edge-to-edge spacing d between spheres is larger than the spheres' radius r (i.e., d ≥ r). However, even for smaller distances, the SDA may provide satisfactory approximated results [33], though in general, for a spacing between the spheres smaller than their radius (i.e., 0 < d < r), more accurate results may involve multipole field contributions [44,[46][47][48].
Suppose now that the array in Fig. 1(a) supports a mode with wavevector k B k xx k yŷ k zẑ . Consequently, each microsphere will have a magnetic dipole moment equal to m n m 0 expik B · d n . Then, the local magnetic field H loc in Eq. (1) acting on a microsphere at position r 0 is produced by all the spherical microparticles of the array except the considered microparticle plus the external incident field to the array (i.e., we are considering all the mutual magnetic couplings). Accordingly, whereG ∞ r 0 ; r 0 ; k B · m 0 is the magnetic field produced by all the other microspheres but the one at position r 0 , and G ∞ r 0 ; r 0 ; k B represents the regularized GF. The latter is defined asG ∞ r; r 0 ; k B G ∞ r; r 0 ; k B − Gr; r 0 , which is not singular at r r 0 , where is the magnetic-field dyadic GF for the phased period array of dipoles, and Gr m ; r n is the dyadic GF of the background medium (see Appendix A for more details on the computation ofG ∞ r; r 0 ; k B through the Ewald method). If the array is placed in a homogeneous medium, the GF (field produced at r m by a unit magnetic dipole at r n ) simply is where r mn jr m − r n j and r mn is the vector from the source at r n to the observer at r m . Substituting then the expression for the local field at position r 0 given in Eq. (4) into Eq. (1), we get which leads to the linear system Mode analysis in the 3D periodic array is then performed by computing the complex zeroes of the homogeneous linear system in Eq. (8), after having imposed no excitation source (i.e., H inc r 0 0). This leads to the computation of the complex zeroes of the determinant of Ak B (see [41] for more details on the properties of Ak B in 3D periodic arrays).

MODES WITH REAL AND COMPLEX WAVENUMBERS IN 3D PERIODIC ARRAYS OF TIO 2 MICROSPHERES
In this section, we show modes with real and complex wavenumber for traveling modes along the z-direction (i.e., k B k zẑ , with k z β z iα z ) with dipole moments polarized along the z-direction (longitudinal polarization, "L-pol," with respect to the mode traveling direction) or either the x-or y-direction (transversal polarization, "T-pol," with respect to the mode traveling direction) in a 3D periodic array of TiO 2 microspheres embedded in free space (i.e., k k 0 , ε h 1). In this paper we consider arrays with cubic lattice (i.e., a b c), and with the parameters reported in Table 1 for different volumetric filling factor f vol 4πr 3 ∕3 a 3 (for simple cubic lattices, the dense packing maximum filling factor is f vol;max π∕6 ≈ 52%), in the frequency range between 250 GHz and 350 GHz, in which the magnetic effects are dominant with respect to the electric ones, as explained in Subsection 2.A.

A. Transversal Polarization (T-pol)-Structure II
The dispersion diagrams for modes for T-pol for structure II are reported in Fig. 2. Only modes with α z ≥ 0, i.e., those whose power flow is toward the positive z-direction, are shown. Three different curves are shown: (i) solid curves are computed by using the SDA method described in Subsection 2.B for the 3D periodic array embedded in free space (i.e., ε h 1); (ii) dashed curves are computed by using the SDA method assuming the same 3D periodic array to be embedded in a host medium with a relative effective permittivity ε h ε eff ≈ 2.38 (i.e., artificially including in the SDA an effective host permittivity due to the high value of the permittivity of the TiO 2 microspheres) computed through Maxwell Garnett formulas as explained and shown in Subsection 5.B [this method will be referred to as SDA with correction (SDA-WC) throughout the paper]; (iii) dotted-circled curves are computed by fitting the magnetic field calculated via an HFSS full-wave simulation as explained in Section 4. Note the small frequency shift between the modal analyses and the fitting from the full-wave simulation. Moreover, it can be observed that the use of the SDA-WC method provides us with results in better agreement with the fitting curve than the SDA result. This result could be further improved by modeling the particle with dual (electric and magnetic) dipole approximation, as will be investigated in future work. At low frequencies, mode 1 (blue line) is characterized by a phase constant larger than the free space wavenumber and by a small attenuation; increasing the frequency, the dispersion curve bends exhibiting large phase constant. Further increasing frequency, mode 1 experiences a bandgap with a strong attenuation; finally, at higher frequencies mode 1 reenters a propagation band with small attenuation. Mode 2 (green line) at low frequencies is characterized by a large phase constant at the edge of the Brillouin zone and a large attenuation constant. Increasing the frequency, the attenuation constant decreases (remaining yet larger than that of mode 1 in the bandgap) in a small frequency region, with phase constant approaching very small values. At higher frequencies, mode 2 is characterized by very small phase constant and large attenuation constant. Other modes with normalized attenuation constant larger than α z a∕π 1.5 are present but not reported here since they dramatically decay as α z ≫ k 0 . The occurrence of two transversely polarized ("doubly degenerated," that is, admitting any transverse polarization) modes with moderately low attenuation constant α z reveals a certain degree of nonlocality in the 3D array, i.e., the presence of spatial dispersion. This was noticed in [41,49,50] for a 3D array of plasmonic nanoparticles. In [50], spatial dispersion was noticed also for high dielectric spheres (though the longitudinal case was not analyzed) and this depends on how dense the array is. Though in this study trends are similar to those in [41], the k − k z dispersion diagrams are different, because of different losses, permittivity frequency-dispersion, and especially, the fundamental difference in the exploited resonance in the spheres. We have also checked that a slightly smaller level of losses for TiO 2 , whether one takes a frequency-dependent value or the highest one at 500 GHz, as discussed in Subsection 2.B, does not quantitatively affect the results shown here and in the next sections. In particular, the only change would be in the range 0.235 < k 0 a∕π < 0.245 for the T-pol case and 0.275 < k 0 a∕π < 0.285 for L-pol in Subsection 3.B, where results would approach further the lossless case discussed in [41].
In Fig. 3, we show the evolution of the modal wavenumber k z β z iα z for varying frequency. Namely the trajectories of the complex propagation constants are tracked in the complex plane β z , α z in Fig. 3(a), normalized to the period a, and in Fig. 3(b), normalized to the free space wavenumber k 0 . Notice that a mode whose power flow is toward either the positive or negative z-direction, for which either α z > 0 or α z < 0, is forward (from the phase progression point of view) when β z α z > 0, whereas it is backward when β z α z < 0. By looking at Figs. 2 and 3, the forward mode 1 and the backward mode 2 could be guided in the structure. As a final remark, in the lossless case, these two modes assume the same modal wavenumber k z in the range 0.24 < k 0 a∕π < 0.245, as was noticed in Figs. 2 and 3 of [41] , and they may become degenerate. Note that in the lossy case, these two modes do not assume the same modal wavenumber k z for any frequency, as can be observed in Figs. 2 and 3 in this work.

B. Longitudinal Polarization (L-pol)-Structure II
The dispersion diagrams for structure II in Table 1 are shown in Fig. 4 for both the real and the imaginary parts of the wavenumber k z β z iα z in the case of L-pol (longitudinal with respect to the mode traveling direction z) computed with both the SDA and the SDA-WC methods. Only the modes with α z ≥ 0, i.e., those with power flow toward the positive z-direction, are shown. Notice again that other modes with normalized attenuation constant larger than α z a∕π 1.5 are present but not reported here since guided modes can travel a significant distance only when their attenuation constant is small, or α z ≪ k 0 . The small value of the slope of the curve in Fig. 4(a) reveals a certain degree of nonlocality, and thus spatial dispersion [49]. A similar, but not identical, result was observed in [41] for the case of a 3D crystal of plasmonic nanoparticles. However, mainly due to the fact that here we use different materials (different losses and material dispersion) and exploit a different resonance (due to b 1 instead of a 1 ), the slope of the curve in Fig. 4(a) is slightly larger than the one observed in [41], revealing a slightly larger spatial dispersion. The evolution of the modal wavenumbers varying frequency is shown in the complex k z plane in Fig. 5. At low frequencies, mode 1 (blue curve) is characterized by very small phase constant and large attenuation constant. Increasing frequency, the phase constant increases, whereas the attenuation constant decreases. Note that since the intermicrosphere longitudinal coupling is due to the near field coupling of the magnetic field, as in [51,52], the effect of the permittivity of the environment is weak, explaining why the two curves shown are almost superimposed.

FIELD FITTING
In general, considering only modes with a moderately low attenuation constant α z , a plane wave impinging on the finite thickness composite material in the inset in Fig. 6 could excite two modes with T-pol and one mode with L-pol, with different amplitudes (as also briefly discussed in [38,41] for a 3D array of plasmonic nanoparticles). For a normal incident plane wave, only modes with T-pol could be excited. Our aim is to show that, even though the analyzed structure has a certain degree of spatial dispersion (as previously described in Section 3, and in [38,41,49,50]), mode 1 (T-pol) in Subsection 3.A is "dominant" (i.e., it contributes mostly to the field in the array). Accordingly, a wave in the composite material, in the case of transverse polarization, could be described with good approximation as a TEM wave in a homogeneous material, which can in turn be represented by an effective refractive index and permeability. Though, strictly speaking, the effective permeability characterizing the crystal as a homogenized material has a small degree of nonlocality (spatial dispersion) [38,49].  From the HFSS full-wave simulation of 11 layers of TiO 2 microspheres (see inset in Fig. 6), illuminated by a normally incident plane wave traveling toward z, with magnetic field polarized along y, we extract the y component of the magnetic field (1 point per layer, at the center of each sphere) for structure II in Table 1. As stated in Subsection 3.A, the forward mode 1, and the backward mode 2 could be excited for some frequencies. Therefore, we assume that the total magnetic field could be represented as the superposition of two direct ("") and two reflected ("−") waves, pertaining to two modes (T-pol) with complex wavenumber k z;f and k z;b , traveling along the z directions as follows: where n 1; …; 11, and A , A − , B , B − , k z;f , and k z;b are all complex valued unknowns, with Rek z;f > 0 (i.e., forward, mode 1) and Rek z;b < 0 (i.e., backward, mode 2), assuming power flow toward the positive z-direction (i.e., Imk z;f , Imk z;b > 0). The H y n field obtained via HFSS is shown in Fig. 6 (black crosses) at 290 GHz (i.e., the frequency for which Imk z;b is minimum; for other frequencies Imk z;b is even larger). Though by inspection the HFSS data in Fig. 6 seems to have a single exponential decay, we perform a curve fitting by using Eq. (9), and we report the field fitting curves (both magnitude and phase) in comparison to the extracted full-wave simulation field in Fig. 6 (blue curves). It can be observed that the fitting curves are in good agreement with the fields from the HFSS simulation; thus we can conclude that the field can be represented as in Eq. (9) (note, however, that A − and B − are found to be several orders of magnitude smaller than A and B , with A orders of magnitude larger than B ). Moreover, we have reconstructed the dispersion diagram for T-pol in the frequency range between 250 GHz and 350 GHz from the information of the fitted wavenumbers k z;f and k z;b using Eq. (9) as the blue and green dotted-circled curves in Fig. 2 (where the circles represent the analyzed frequencies), which were found to be in good agreement with the ones computed by using the SDA-WC method in Subsection 3.A. We have also observed that, in the analyzed frequency points, A ≫ A − , B ≫ B − , thus there is no significant reflection at the two metamaterial-air interfaces, and that A ≫ B , thus mode 1 is "dominant" (i.e., it contributes mostly to the field in the array). Therefore k z;f could be also retrieved by performing a single exponential fitting of the data, which would lead to the same k z;f shown in Fig. 2.

EFFECTIVE MEDIUM THEORY: RETRIEVAL METHODS
In Section 4 we have observed that mode 1 (T-pol) in Subsection 3.A is "dominant"; thus the composite material,  in the case of transverse polarization, may be treated as an effective medium.
A. Retrieval from Mode Analysis: SDA and SDA-WC Using the modal results computed in Section 3, it is possible to retrieve the effective refractive index as n eff k z ∕k 0 , where k z is the wavenumber of the "dominant" mode (i.e., mode 1 for T-pol).
The effective relative permeability from the SDA method can be computed as μ eff n eff 2 ; in the SDA-WC case, instead, μ eff n eff ∕ ε eff p 2 , where ε eff is as defined in the next subsection.

B. Maxwell Garnett Formulas
In general, Maxwell Garnett theory [53,54] can be applied to retrieve the relative effective permittivity and permeability of a composite medium as where N D f vol ∕V N , and V N 4πr 3 ∕3 is the microsphere volume. Here, α ee and α mm are the electric and magnetic Mie polarizabilities as in Eq. (3), and radiative losses are subtracted in Eq. (10) since the microspheres are on a regular lattice with a < λ 0 ∕2 [55]. Then, the effective refractive index can be calculated as n eff ε eff μ eff p . Note that these formulas will be applied to arrays of microspheres embedded in free space, as described in Section 3.
We show in Fig. 7 the relative effective permittivity retrieved by using Eq. (10) for the three arrays in Table 1. Note that there is not a resonant effect; however, the effective value differs largely from the one of free space (i.e., ε h 1): it is observed that the effective permittivity is almost constant in the analyzed frequency range for structures II and III. In particular, for structure II, this constant effective value has been assumed to be ε eff ≈ 2.38 (as previously mentioned in Subsection 3.1) and included in the SDA-WC method (note also that the imaginary part is quite small, so the permittivity has been assumed to be real in the SDA-WC). As a last remark, we want to point out that the formulas in Eq. (10) work well for diluted systems, i.e., for f vol < 0.4 [56]. Therefore, to apply the SDA-WC method for structure I (for which f vol ≈ 0.49 > 0.4), one should use the effective permittivity computed, for example, through the method discussed in the next subsection.

C. Nicolson-Ross-Weir (NRW) Retrieval Method
The transmission (T) and reflection (R) coefficients for a stack of N 11 layers (Fig. 6) computed by the full-wave HFSS simulation are here used to retrieve the effective refractive index of the composite material by using the NRW method [57][58][59][60][61]. Treating the composite slab as a uniform continuous medium with the same thickness t, according to NRW, the complex effective refractive index and impedance can be retrieved by where q is an integer to be determined, and t Na, with N denoting the number of layers and a the separation between two contiguous layers. We address the reader to [60] for guidelines on how to choose q and in Eq. (11). Then, it is possible to retrieve the relative effective permeability as μ eff n eff Z eff ε 0 ∕μ 0 p .

ARTIFICIAL MAGNETISM PROPERTIES FROM NONMAGNETIC MATERIALS
In this section we show the comparison of the effective parameters results retrieved by using the methods described in Section 5. Note that the NRW result based on reflection and transmission computed through an HFSS full-wave simulation is the most accurate, and its computational burden is the heaviest. The second most accurate result is computed by using the SDA-WC method (again, this result could be further improved by modeling the particle with dual electric and magnetic dipole approximation as mentioned in Subsection 3.A).

A. Comparison of Retrieval Methods-Structure II
We start by comparing the relative effective permeability in Fig. 8 and the effective refractive index in Fig. 9 for structure II in Table 1. Note how the different methods are in good agreement with each other. In particular, referring to Fig. 8, large values of Reμ eff (both positive and negative) can be generated (with large Imμ eff in the bandgap region outlined  Table 1 for T-pol retrieved through Maxwell Garnett formulas in Eq. (10).  Table 1 for T-pol retrieved through different methods.  Table 1 for T-pol retrieved through different methods.  Table 1 for T-pol according to NRW-HFSS retrieval method.  Table 1 for T-pol according to NRW-HFSS retrieval method.
in Subsection 3.A); moreover, notice how the permeability retrieved through the SDA method matches the full-wave NRW result: indeed, in the SDA case, the computed refractive index shown in Fig. 9 corresponds to the magnetic part n eff m of the total refractive index, which can be expressed as n eff n eff e n eff m ε eff p μ eff p . Notice how using the SDA-WC method provides a better result for the refractive index than the SDA (with respect to the NRW result), while still providing a fair estimation of the effective relative permeability in Fig. 8. Table 1 We continue by comparing the relative effective permeability in Fig. 10 and the effective refractive index in Fig. 11 for the three structures in Table 1 (parametric analysis with respect to the spatial period a of the array) retrieved through NRW. Observe how an increase in the periodicity impacts in a decrease of the losses (i.e., imaginary part) of the system. Also, observe the frequency blue shift for increasing periodicity.

CONCLUSION
This work provides a description of modal analysis in 3D periodic arrays of titanium dioxide (TiO 2 ) microspheres, aiming at showing artificial magnetism (i.e., relative permeability ≠ 1) generation from nonmagnetic materials, including the arise of negative permeability frequency bands. We observe the presence of one longitudinal mode and two transverse modes, one forward and one backward. In case of T-pol, mode 1 has been found to be dominant, allowing for the description of the composite material in terms of homogenized effective refractive index and relative permeability. We have shown that the single magnetic dipole approximation provides accurate results as long as the dipoles scatter in a homogeneous environment where the host permittivity includes that of the microspheres and it is calculated by using Maxwell Garnett mixing formulas. These parameters have been retrieved by four different methods, in good agreement with each other, confirming the homogenization procedure to be possible. Artificial magnetism (including negative permeability) from nonmagnetic materials has been proven at millimeter and submillimeter wave frequencies. Our results are corroborated by a very recent experimental work [62] where some degree of artificial magnetism in composite materials made of TiO 2 microparticles has been demonstrated in the low-THz regime.

APPENDIX A: DYADIC EWALD GREEN'S FUNCTION FOR 3D PERIODIC ARRAYS OF MICROSPHERES MODELED AS MAGNETIC DIPOLES
The dyadic form of the GF in Eq. (4) representing the magnetic field H loc at a general observation point r due to the 3D lattice of magnetic dipoles m n but the one at r 0 is G ∞ r; r 0 ; k B k 2 I ∇∇ G ∞ r; r 0 ; k B : Here G ∞ r; r 0 ; k B is the regularized scalar GF, not singular at r r 0 , as in Eq. (11) in [41]. Note the difference between Eq. (A1) in this paper and Eq. (10) in [41] due to the magnetic dipole formulation.
The Ewald representation for the regularized dyadic GF is thenG ∞ r; r 0 ; k B G spectral r; r 0 ; k B G spatial r; r 0 ; k B ; (A2) and G spectral r; r 0 ; k B k 2 I ∇∇G spectral r; r 0 ; k B ; (A3) G spatial r; r 0 ; k B k 2 I ∇∇ G spatial r; r 0 ; k B ; (A4) where representations for the terms in Eqs. (A3) and (A4), in the specialized case r r 0 0, are given in [41].