Terahertz scattering by subwavelength cylindrical arrays

We demonstrate the use of a full-wave electromagnetic field simulator to verify terahertz (THz) transmission-mode spectroscopic measurements of periodic arrays containing subwavelength cylindrical scatterers. Many existing THz scattering studies utilize analytical solutions, which were developed for a single scatterer. For multiple scatterers, a scaling factor equal to the number of scatterers is applied, accounting for interference between far-field radiative contributions from those scatterers but not their near-field mutual coupling. Consequently, analytical solutions do not accurately verify measurements. Conversely, results from the fullwave electromagnetic field simulator elucidate our measurements well, and provide an important insight into how the scattering behavior of cylindrical scatterers is influenced by test conditions. © 2011 Optical Society of America OCIS codes: (300.6495) Spectroscopy, terahertz; (290.5820) Scattering measurements; (050.1755) Computational electromagnetic methods; (120.4290) Nondestructive testing. References and links 1. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley & Sons, 1983). 2. L. M. Zurk, B. Orlowski, D. P. Winebrenner, E. I. Thorsos, M. R. Leahy-Hoppa, and L. M. Hayden, “Terahertz scattering from granular material,” J. Opt. Soc. Am. B 24, 2238–2243 (2007). 3. A. Bandyopadhyay, A. Sengupta, R. B. Barat, D. E. Gary, J. F. Federici, M. Chen, and D. B. Tanner, “Effects of scattering on THz spectra of granular solids,” Int. J. Infrared Millim. Waves 28, 969–978 (2007). 4. Y.-C. Shen, P. F. Taday, and M. Pepper, “Elimination of scattering effects in spectral measurement of granulated materials using terahertz pulsed spectroscopy,” Appl. Phys. Lett. 92, 051103 (2008). 5. M. Franz, B. M. Fischer, and M. Walther, “The Christiansen effect in terahertz time-domain spectra of coarsegrained powders,” Appl. Phys. Lett. 92, 021107 (2008). 6. J. E. Bjarnason, T. L. J. Chan, A. W. M. Lee, M. A. Celis, and E. R. Brown, “Millimeter-wave, terahertz, and mid-infrared transmission through common clothing,” Appl. Phys. Lett. 85, 519–521 (2004). 7. M. A. Kaliteevski, D. M. Beggs, S. Brand, R. A. Abram, J. R. Fletcher, G. P. Swift, and J. M. Chamberlain, “Propagation of electromagnetic waves through a system of randomly placed cylinders: The partial scattering wave resonance,” J. Mod. Opt. 53, 2089–2097 (2006). 8. J. R. Fletcher, G. P. Swift, D. C. Dai, J. A. Levitt, and J. M. Chamberlain, “Propagation of terahertz radiation through random structures: An alternative theoretical approach and experimental validation,” J. Appl. Phys. 101, 013102 (2007). 9. Y. L. Hor, J. F. Federici, and R. L. Wample, “Nondestructive evaluation of cork enclosures using terahertz/millimeter wave spectroscopy and imaging,” Appl. Opt. 47, 72–78 (2008). 10. P. Y. Han, G. C. Cho, and X.-C. Zhang, “Time-domain transillumination of biological tissues with terahertz pulses,” Opt. Lett. 25, 242–244 (2000). 11. M. Reid and R. Fedosejevs, “Terahertz birefringence and attenuation properties of wood and paper,” Appl. Opt. 45, 2766–2772 (2006). #144099 $15.00 USD Received 16 Mar 2011; revised 25 Apr 2011; accepted 26 Apr 2011; published 9 May 2011 (C) 2011 OSA 23 May 2011 / Vol. 19, No. 11 / OPTICS EXPRESS 10138 12. J. Beckmann, H. Richter, U. Zscherpel, U. Ewert, J. Weinzierl, L.-P. Schmidt, F. Rutz, M. Koch, H. Richter, and H.-W. Hübers, “Imaging capability of terahertz and millimeter-wave instrumentations for NDT of polymer materials,” in Proceedings of 9th European Conference on Nondestructive Testing (NDT), R. Diederichs, ed. (The e-Journal of Nondestructive Testing, 2006). 13. K. Naito, Y. Kagawa, S. Utsuno, T. Naganuma, and K. Kurihara, “Dielectric properties of eight-harness-stain fabric glass fiber reinforced polyimide matrix composite in the THz frequency range,” NDT & E Int. 42, 441– 445 (2009). 14. K. Naito, Y. Kagawa, S. Utsuno, T. Naganuma, and K. Kurihara, “Dielectric properties of woven fabric glass fiber reinforced polymer matrix composites in the THz frequency range,” Compos. Sci. Technol. 69, 2027–2029 (2009). 15. X.-L. Tang, Y.-W. Shi, Y. Matsuura, K. Iwai, and M. Miyagi, “Transmission characteristics of terahertz hollow fiber with an absorptive dielectric inner-coating film,” Opt. Lett. 34, 2231–2233 (2009). 16. Y. Dikmelik, J. B. Spicer, M. J. Fitch, and R. Osiander, “Effects of surface roughness on reflection spectra obtained by terahertz time-domain spectroscopy,” Opt. Lett. 31, 3653–3655 (2006). 17. T. Kleine-Ostmann, C. Jansen, R. Piesiewicz, D. Mittleman, M. Koch, and T. Kürner, “Propagation modeling based on measurements and simulations of surface scattering in specular direction,” in Proceedings of Joint 32nd International Conference on Infrared and Millimetre Waves, and 15th International Conference on Terahertz Electronics, M. J. Griffin, P. C. Hargrave, T. J. Parker, and K. P. Wood, eds., vol. 1, pp. 408–410 (2007). 18. R. Piesiewicz, C. Jansen, D. Mittleman, T. Kleine-Ostmann, M. Koch, and T. Kürner, “Scattering analysis for the modeling of THz communication systems,” IEEE Trans. Antennas Propag. 55, 3002–3009 (2007). 19. J. Pearce, Z. Jian, and D. M. Mittleman, “Statistics of multiply scattered broadband terahertz pulses,” Phys. Rev. Lett. 91, 043903 (2003). 20. J. Pearce and D. M. Mittleman, “Using terahertz pulses to study light scattering,” Phys. B 338, 92–96 (2003). 21. Z. Jian, J. Pearce, and D. M. Mittleman, “Characterizing individual scattering events by measuring the amplitude and phase of the electric field diffusing through a random medium,” Phys. Rev. Lett. 91, 033903 (2003). 22. S. Mujumdar, K. J. Chau, and A. Y. Elezzabi, “Experimental and numerical investigation of terahertz transmission through strongly scattering sub-wavelength size spheres,” Appl. Phys. Lett. 85, 6284–6286 (2004). 23. X. J. Zhong, T. J. Cui, Z. Li, Y. B. Tao, and H. Lin, “Terahertz-wave scattering by perfectly electrical conducting objects,” J. Electromagn. Waves Appl. 21, 2331–2340 (2007). 24. R. A. Cheville, M. T. Reiten, R. McGowan, and D. R. Grischkowsky, Applications of Optically Generated Terahertz Pulses to Time Domain Ranging and Scattering (Springer-Verlag, 2003). 25. G. M. Png, R. J. Falconer, B. M. Fischer, H. A. Zakaria, S. P. Mickan, A. P. J. Middelberg, and D. Abbott, “Terahertz spectroscopic differentiation of microstructures in protein gels,” Opt. Express 17, 13102–13115 (2009). 26. C. A. Balanis, Antenna Theory: Analysis and Design (Harper & Row Publishers, 1982). 27. E. Hecht, Optics (Addison-Wesley Publishing Company, 2002). 28. H. C. van de Hulst, Light Scattering by Small Particles (John Wiley & Sons, Inc., 1957). 29. M. Kerker, Scattering of Light and Other Electromagnetic Radiation (Academic Press, New York, 1969). 30. E. F. Knott, J. F. Shaeffer, and M. T. Tuley, Radar Cross Section (SciTech Publishing Inc., 2004). 31. R. Piesiewicz, T. Kleine-Ostmann, N. Krumbholz, D. Mittleman, M. Koch, and T. Kürner, “Terahertz characterisation of building materials,” Electron. Lett. 41, 1002–1004 (2005). 32. M. Naftaly and R. E. Miles, “Terahertz time-domain spectroscopy of silicate glasses and the relationship to material properties,” J. Appl. Phys. 102, 043517 (2007). 33. C. A. Balanis, Advanced Engineering Electromagnetics (John Wiley & Sons, 1989). 34. G. Chattopadhyay, J. Glenn, J. J. Bock, B. K. Rownd, M. Caldwell, and M. J. Griffin, “Feed horn coupled bolometer arrays for SPIRE—design, simulations, and measurements,” IEEE Trans. Microwave Theory Tech. 51, 2139–2146 (2003). 35. L. Liu, S. M. Matitsine, Y. B. Gan, and K. N. Rozanov, “Effective permittivity of planar composites with randomly or periodically distributed conducting fibers,” J. Appl. Phys. 98, 063512 (2005). 36. C. Hafner, The Multiple Multipole Program (MMP) and the Generalized Multipole Technique (GMT), pp. 21–38 (Elsevier Science B.V., 1999).


Introduction
An "infinitely long circular cylinder" is an approximate model for a circular cylinder that has length l diameter d, l wavelength λ , and is orientated orthogonal to the direction of propagation of an incident electromagnetic beam.Many naturally occurring organisms and minerals, such as some viruses and asbestos, can be modeled accurately by infinitely long circular cylinders [1].Equally ubiquitous in everyday life are synthetic objects that too can be represented as infinitely long cylinders, such as fabric fibers or hairline scratches on hard surfaces.Arrays of infinitely long cylindrical objects can display distinct scattering and transmission charac-teristics, such as birefringence or dichroism, which allow for their detection and identification (e.g.birefringent proteins are luminous when dyed appropriately).The prevalence of infinitely long circular cylinders make them interesting and useful objects to study.A study of scattering from arrays of infinitely long circular cylinders in the terahertz frequency range (THz or T-ray, frequency ν = 0.1-10 THz, wavelength λ = 3 mm-30 μm) forms the basis of this paper.
Most scattering studies in the THz frequency range have been conducted on samples with dimensions that are either at the same order as the THz wavelengths or at subwavelength values.Two popular materials reported in existing literature are powders of different dimensions [2][3][4][5] and fibers from clothing [6][7][8].Other materials used in THz scattering studies include cork enclosures [9], pork-fat and onion cells [10], wood and paper [11], polymers [12], fiberglass cloth [13,14], and hollow fiber with a film of absorptive dielectric inner-coating [15].Surface roughness studies have also emerged from reflection mode THz scattering studies for communication purposes [16][17][18].Statistical studies of THz scattering in random media [19][20][21][22], computer models for THz scattering [8,23], and THz ranging akin to that used in radar ranging [24] are other areas of research.
Within the overall body of literature pertaining to cylindrical scatterers, few theoretical models have been presented to explain the fibrillar scattering phenomena observed.One notable exception is in [8], where a phase distribution function is used to model the random structure of various fabrics.Although more literature exists for theoretical modeling of spherical scatterers, the overall landscape of THz scattering studies (regardless of the samples' shapes) is sparse.This gap is the first of two motivations for this study.
The second motivation for this study stems from our past experiments involving soft cylindrical (fibrillar) protein microstructures [25].In order to model the fibrillar proteins, an assumption was made when using analytical solutions of Rayleigh scattering that neighboring fibrils interfere constructively in the scattering and transmission of the incident electromagnetic wave.Although the resultant model sufficiently explained the observed THz measurement, it can be improved.
From antenna theory it is known that for a linear antenna array with all elements in phase, the total electric field radiating in the broadside direction is the result of constructive interference and mutual coupling between adjacent antenna elements.This results in azimuthal and zenithal radiation patterns with main and side lobes [26].If an array of circular cylinders is only modeled using the superposition of analytical single element solutions, then the influence of mutual coupling is lacking.In this paper, we present a study utilizing a 3D full-wave electromagnetic field solver that can account for both interference (far-field) and mutual coupling (near-field).Our results are compared with THz transmission measurements of subwavelength cylindrical arrays made from fiberglass.We show that the simulation results accurately elucidate the observed transmission measurements.

Fibrillar samples and experimental results
The transmission mode time domain spectroscopy (TDS) THz system used in this paper is based on the generation of a THz pulse using an ultrafast (90 fs) near-infrared laser (Spectra Physics Tsunami) and a photoconductive antenna, with coherent detection via electro-optic (EO) sampling.Cylindrical (fibrillar) scatterers are made from fiberglass strands obtained from woven fiberglass cloth (plain weave E glass, weight 7.5 oz/yd 2 ) as shown in Fig. 1(a).
Each fiberglass strand is approximately 10 μm in diameter (λ min /3, λ max /300) and is 15 mm long (500λ min , 5λ max ).Since the length of each strand is 1500 times larger than the diameter, then each strand can be approximated by an infinite cylinder according to the definition given in Section 1.The strands are aligned as a near-periodic array as shown in Fig. 1(b).The term "near-periodic" is used as it is acknowledged that the strands' small cross-sectional dimensions  , with propagation vector k normal to the long axis of the cylinder (z-axis).The electric field vector E is parallel to the z-axis.After [1].make handling difficult, resulting in imperfections in the array structure.These imperfections are taken into consideration later in Section 5.
Four arrays are deployed in this study.Each array is clamped between two polyethylene frames during measurements as shown in Fig. 1(c) so that the array is only exposed to air in the region of the incident THz beam.One of the frames is covered with aluminum foil to prevent transmission of THz through the frames.As illustrated in Fig. 1(d), the array is placed orthogonally in the path of the vertically polarized incident THz beam.The arrays are placed either parallel ( ) or perpendicular (⊥) to the vertically polarized THz beam. Figure 1(e) shows how the parallel orientation relates to the azimuthal and zenithal directions of an infinite circular cylinder.
The THz frequency response of two of the four fiberglass array samples are presented in Fig. 2(a).Figures 2(b) and 2(c) present the measured optical properties of the two arrays.The term extinction coefficient (with symbolic notation α) is used instead of absorption coefficient to account for scattering.Each THz measurement is averaged between 16 to 25 scans.Results are verified using two different THz-TDS systems (at the University of Adelaide, and at the University of Leeds).
The samples perpendicular to the THz beam have very similar extinction coefficients.The Both parallel samples block more THz transmission than the perpendicular samples, suggesting agreement with the filtering capabilities of grid arrays, such as wire grid arrays in the transverse magnetic (TM) orientation [27].However, the two parallel samples have differing extinction coefficient profiles.Sample 1 not only appears to block more THz transmission, its extinction coefficient profile slopes downwards beyond 2 THz.This slope is not observed in sample 2. Visual checks of both samples do not reveal any significant physical variations in the two samples.Repeat measurements of all four samples in the both orientations reveal similar inconsistencies for the parallel orientation.The decrement in extinction coefficient for sample 1 (parallel) occurs within the estimated system bandwidth of ≈ 0.2-2.5 THz, hence it can be argued that the decrement is not an artifact.Further investigations are necessary to ascertain if the sloping feature is an artifact, and to answer the question as to why the parallel orientation appears more susceptible to inconsistencies.Outcome of investigations are presented in Section 6.
The parallel arrays have uneven effective n profiles that vary more over frequency (sample 1: n = 2.13 at 0.5 THz; 1.74 at 2 THz), indicating dispersion in the fiberglass strands.The difference in refractive index between the parallel and perpendicular orientations suggests that the strands are birefringent, and that the measured refractive indices are effective values that depend in geometry rather than solely bulk values.Furthermore, given the more significant increase in extinction for the parallel samples, the strands may also be dichroic.Alternatively, the increase in α could be due to strong THz scattering in this orientation.In order to elucidate if scattering causes the increase in α for the parallel orientation, and to investigate why the two parallel samples have differing transmission profiles, Mie scattering is investigated in the next Section.

Analytical solutions: Mie and Rayleigh scattering
Analytical solutions for Mie scattering-of which Rayleigh scattering is a special case where object size wavelength λ -can be found in numerous pieces of literature [1,28,29], hence they will not be elaborated on here.However a summary of salient equations is provided as follows.Referring to Fig. 1(e), when an electromagnetic wave E inc propagating in a medium with complex refractive index ni is incident normal (θ = 90 • ) to the long (z) axis of an infinite cylinder with complex refractive index nt , some of the power from the incident wave is scattered back by the target object ("backscatter").The measure of this backscattered power is the scattering cross section C sca (units: m 2 ), which is analogous to the monostatic radar cross section (or echo area) [30].The backscattering cross sections for the parallel and perpendicular cases are defined in Eqs. ( 1) and (2): where wavenumber k = 2π/λ , r = cylinder's radius, l = cylinder's length, n = n − iκ, κ = imaginary component of the refractive index or absorption index = αλ /(4π), m = nt / ni , J s (p) and Y s (p) are the Bessel functions of the first and second kind respectively with integral order s ∈ R, J s (p) is the first derivative of J s (p), and H (1)  s (p) = J s (p) + iY s (p) is the Hankel function.
In literature relevant to scattering, there is no distinction between the effective and the bulk refractive indices of a material, regardless of the material's shape or size [1,28,29].In [29], examples seem to suggest the use of the bulk refractive index.In the framework of the fullwave electromagnetic simulation, it is clear that the input material parameters have to be taken from the bulk material properties.The results of the simulation for the particular geometry then yield information about the effective material properties.

Comparison between experimental and analytical results
To explain the increase in extinction coefficients α with frequency as shown in Section 2, the scattering cross sections C sca, and C sca, ⊥ in Eqs. ( 1) and ( 2) are investigated.However, knowledge of n is needed for m in Eqs. ( 3)- (5).Since the measured effective material parameters n and α of the fiberglass strands (displayed in Fig. 2) are dependent on orientation, they are not suitable for use in n.One alternative is to substitute the n and α values of the fiberglass strands with those of bulk glass reported in literature.A second alternative is to use n and α of bulk fiberglass cloth.
The THz characterization of a wide range of bulk glass types, ranging from optical grade glass to generic window glass, has been reported in literature [31,32].However, the number of references to bulk fiberglass cloth is sparse [13].To supplement existing data, measurements of bulk fiberglass cloth are included in this study.The measured n and α values are: (at 0.2 THz) n bulk = 2.1, κ bulk = 0.096; (at 1 THz) n bulk = 1.95, κ bulk = 0.287.These values are within the range reported by other authors for bulk glass and bulk fiberglass cloth, hence are deemed acceptable.
Figure 3(a) presents the C sca, and C sca, ⊥ plots on a logarithmic scale when E inc is polarized parallel to the z-axis in Fig. 1.For the frequencies under investigation C sca, ⊥ C sca, .This means that a fiber oriented perpendicular to the incident E field has a much smaller scattering area than one that is oriented parallel to E. In accordance with the definition of backscattering cross section in Section 3, a smaller scattering area implies that a weaker signal is scattered back to the illumination source, resulting in a stronger transmitted signal.This trend is consistent with those in Figs.2(a   Although the trend in Fig. 3(a) is consistent with those in Figs.2(a) and 2(c), the analytical model does not reveal if the decrease in measured α from 2 THz onwards for the parallel case is an artifact.The relative difference between C sca, and C sca, ⊥ is also exaggerated in the analytical model.These gaps in the analytical model motivate the exploration of other modeling techniques.One possible technique is to account for constructive and destructive interference from the array elements by using an array factor; this is elaborated on in the next Subsection.Another technique is to use a numerical electromagnetic simulation tool as described in Section 5.

Modeling an array with analytical solutions
To model multiple scatterers, Mie scattering equations are usually multiplied by a scaling factor equal to the number of scatterers under test [3,25], assuming constructive interference from all scatterers in the far field.When a scaling factor is applied to both plots in Fig. 3(a), the resulting plots still do not resemble the measurements in Fig. 2(c).
One possible way of deriving a multiplication factor is to apply the mathematics used to model antenna arrays, such as the array factor (AF) for N isotropic sources [33]: where , d is the distance between two adjacent sources, β is the difference in phase excitation between adjacent sources, φ is the azimuth (with -90 • being the forward scatter direction, and 90 • being the back scatter direction), ε 0 is the electric permittivity in vacuo, μ 0 is the magnetic permeability in vacuo, and μ r = 1 is the relative magnetic permeability of fiberglass.Assuming that all the fiberglass strands are excited simultaneously by the incident THz plane wave, then β = 0.The plot of the total scattering cross section C sca, total = C sca × AF when diameter d = 1 μm ≈ λ min /100 is presented in Fig. 3(b), where the maximum frequency ν max = 2.788 THz.
The effect of multiplying C sca by AF becomes evident when Figs. 3(a) and 3(b) are compared, showing that AF scales C sca up by 2 orders of magnitude.The inclusion of AF for normal incidence is therefore no different from multiplying the Mie scattering equations by a scaling factor equal to the number of scatterers.This consideration is inadequate for modeling THz scattering for the test scenarios presented in this paper.In the next Section, an alternative model based on full-wave electromagnetic field simulation is presented.The full-wave solver produces results that are a better match with the experimental observations.

Full-wave electromagnetic field solver
The High Frequency Structure Simulator (HFSS, versions 10.1 and 11, produced by Ansoft) is an electromagnetic field simulator that utilizes the Finite Element Method (FEM) as a full-wave numerical solver of Maxwell's equations.The FEM solver finely meshes the virtual problem space to achieve accuracy, but this often leads to heavy utilization of available random access memory (RAM) on a computer.In order to significantly reduce memory use, symmetry and periodicity can be applied in HFSS.This is a practice that has been used by other authors investigating periodic structures in the microwave and THz frequency ranges [34,35].
Figure 4 presents a HFSS geometry modeling an infinite array of infinitely long fiberglass strands.The symmetrical boundaries and periodic boundaries allow the complex test environment containing many infinitely long fiberglass strands to be modeled simply as two 100 μm long strands (each 10 μm in diameter) seated between periodic master/slave side boundaries.To create the infinitely long cylinders (with respect to wavelength) along the z-axis, symmetrical boundaries are applied to the top and bottom airbox walls that touch the ends of the cylinders (symmetrical E boundaries for the parallel case; symmetrical H boundaries for the perpendicular case).Given that HFSS utilizes a 3D solver, the symmetrical boundaries also serve the dual purpose of reducing the problem space into a 2D one, simplifying and speeding up the simulation.It is noted here that other electromagnetic field simulators may be more suited for the 2D problem presented in this study, e.g.COMSOL (by COMSOL AB, also based on FEM but includes a 2D implementation) and the Multiple Multipole Program [36].The latter simulator is most suited for 2D problems, particularly periodic surfaces.
Two strands are used in order to account for the imperfection in aligning the fiberglass strands, such as a slight offset between adjacent strands as shown in Fig. 4(f).The airbox walls that are parallel to the y-z plane are assigned perfectly matched layers (PML) as absorbing boundaries in order to simulate the extension of the problem space into infinity.Since the far-field response is of interest in this study, the incident THz beam is modeled in HFSS as a plane wave as described in Fig. 1(e) instead of as a Gaussian beam.
One limitations inherent in HFSS is that the PML boundaries excludes the airbox from being used for far-field calculations.In order to create a suitable surface for far-field calculations, a second vacuum airbox is used.Not visible in Fig. 4(a) is the second, shorter (along the x-axis) airbox seated tightly inside the first airbox.This airbox is similar in shape to the first, except the walls which are parallel to the y-z plane are at a shorter distance (λ max /2.05) from the cylinders, meaning they do not touch the PML boundaries.None of the walls of the second inner airbox are assigned specific material boundary types, but they are used as Huygen's surface for nearfield to far-field computation.The 0.05 factor used in the length of the second airbox (i.e. the difference in the denominators of λ max /2.05 and λ max /2) is chosen because it is a small number.If another small number is used, such as 0.01, similar results are generated.
One advantage of HFSS is that it allows the reflection and transmission characteristics of the model to be extracted as S-parameters: S 11 and S 21 are the reflection and transmission coefficients respectively.The S-parameters, which are usually expressed in decibels, are defined in [33] either in terms of the electric field E or magnetic field H of a plane wave transmitting from medium 1 into medium 2: where η 1 = (μ 1 /ε 1 ) 1/2 , and η 2 = (μ 2 /ε 2 ) 1/2 .Values of S 11 describe backscattering, whereas S 21 describe transmission; S 21 therefore provides a means for elucidating the behavior of the THz transmission magnitude plots in Fig. 2(a).If S 21 = 1 (or S 21 = 0 dB), then the incident THz signal is 100% transmitted, implying that the transmission magnitude → 0 dB.If S 11 = 1 (or S 11 = 0 dB), then the signal is 100% reflected, implying that the transmission magnitude → ∞.In an ideal condition with no losses, |S 21 | 2 + |S 11 | 2 = 1.The next Section presents the S-parameters from the model in Fig. 4 and compares them with the measured THz transmission magnitudes.

Comparison between experimental and numerical results
Simulations of the HFSS model described in Section 5 are set to cover the frequency range from 0.088 THz to 2.788 THz, with a frequency resolution of 30 GHz.The maximal mesh size on the side of each cylinder is set at one eighth the length of the cylinder; the maximal mesh size at the top and bottom ends is set at a third of the cylinder's diameter.After iterative adaptive mesh refinements, approximately 200,000 tetrahedra are created by HFSS to cover the problem space.
Figure 5(a) presents the S 21 results obtained with the double-cylinder, periodic model shown in Fig. 4. Compared with the measured THz transmission magnitude in Fig. 2(a), the following observations are made with regard to their qualitative and quantitative match.There is a good qualitative match over all frequencies between the measurements and the S 21 results for both the parallel and perpendicular cases.This match is superior to the poor match between the scattering cross section C sca, and C sca, ⊥ plots in Fig. 3(a) with the measured extinction coefficients α in Fig. 2(c).Quantitatively, the match is poor between the measured THz transmission magnitude in Fig. 2(a) and the S 21 results in Fig. 5(a), raising the question whether the single layered HFSS model used is adequate for modeling the actual array, which may contain more than one layer due to imperfections.This question is answered in Section 6.1.
(a) S 21 plots for both parallel and perpendicular orientations, with 1 μm space between adjacent cylinders (b) S 21 plots for both parallel and perpendicular orientations, with 1, 5, 10 and 15 μm spaces between adjacent cylinders Fig. 5. (a) The larger S 21 values for the perpendicular case indicates stronger THz transmission than for the parallel case, agreeing well with the THz transmission magnitude plots in Fig. 2(a) particularly for frequencies below 1 THz.(b) When the space between adjacent strands increases, the S 21 plots for both orientations shift upwards, indicating stronger THz transmission due to "leakage" through the wider spacing.However, the parallel case appears to be more sensitive to the change in strand spacing, shifting more significantly than for the perpendicular case.This resembles the inconsistencies seen in Fig. 2(a) for the two parallel samples.
For the parallel case, the S 21 plot also matches the measured THz transmission magnitudes more accurately at the lower THz frequencies.However, the S 21 plot for the parallel case still does not explain if the upwards-turning curve after 2 THz as seen in Fig. 2(a) for the parallel sample 1 is an artifact.To ascertain if this feature is an artifact, we explore the effect of changing the space between the cylinders in our HFSS simulation.This scenario is realistic, as the fiberglass strands are not secured very tightly to the polyethylene frame in order to prevent them from snapping.The strands may therefore have moved out of place prior to measurements, creating spaces between the strands that are wider than anticipated.
Three values, 5 μm, 10 μm and 15 μm, were arbitrarily chosen to examine the effect of varying the space between adjacent strands.The S 21 plots generated from these new scenarios are presented together with the 1 μm case in Fig. 5(b).Comparing only the perpendicular cases, the S 21 plots shift upwards as the strand spacing increases, with similar profiles except a change in gradient; the changes in gradients are less significant between the 5, 10 and 15 μm cases than between the 1 μm and 5 μm cases.These observations for the perpendicular case suggests that the the fiberglass arrays in the perpendicular orientation are fairly robust to small changes to strand spacing, possibly explaining why no difference in THz transmission magnitude is measured in Fig. 2(a) for the two perpendicular samples.
Like the perpendicular case, the S 21 plots for the parallel case also shift upwards as the strand spacing increases.The extent of shifting is greater than for the perpendicular case, suggesting that the fiberglass arrays in the parallel orientation are more sensitive to small changes to strand spacing, thus the difference in the measured THz transmission magnitudes of the two parallel samples shown in Fig. 2(a).However, there is no evidence in the S 21 plots for the parallel case of the upwards-turning characteristic present in Fig. 2(a).
Observations from the S 21 plots suggests that (i) the fiberglass arrays in the parallel orientation are particularly sensitive to small changes in the strand spacing, manifesting as different THz transmission magnitudes as reported in Section 2 for the two parallel samples; and (ii) the upwards-turn seen in Fig. 2(a) for the parallel sample 1 cannot be reproduced in the HFSS model, hence it can be concluded that the curve beyond ≈ 2.2 THz in Fig. 2(a) is an artifact.However, the lower THz transmission magnitude of the parallel sample 1 when compared to the parallel sample 2 at 2 THz is genuine, and could be due to the aforementioned sensitivity to strand spacing.Alternatively, there could be more than one layer of cylinders in the path of the THz beam.The effect of multiple cylinder layers is investigated in the next Subsection.

Influence of multiple cylinder layers in the array
When Figs. 2(a) and 5(a) are compared quantitatively, the dB values of the plots in Fig. 2(a) are approximately 10 to 15 times larger than those in Fig. 5(a).For example, the magnitude of the perpendicular sample at 1 THz in Fig. 2(a) is ≈ -5 dB, whereas the S 21 value of the perpendicular sample in Fig. 5(a) is ≈ -0.5 dB.This difference may be due to multiple layers of cylinders in the arrays used in the THz measurements.Referring to Fig. 1(b), there is a strong likelihood that multiple layers of fiber threads are present in the path of the THz beam.
To observe the effect of multiple cylinder layers, a periodic model containing only one cylinder as shown in Fig. 6 Figures 6(d) and 6(e) present the S 21 plots of the models containing multiple cylinder layers for both the parallel and perpendicular orientations.The normalized THz transmission magnitudes of the samples as presented previously in Fig. 2(a) are included for comparison.
For frequencies above 0.5 THz, each cylinder layer decreases the S 21 value by ≈ 1 dB for the parallel orientation, and by ≈ 0.5 dB for the perpendicular orientation.This results in a distinct decrease in the S 21 values for the parallel case when compared with the perpendicular case, indicating less transmission for the parallel case that is consistent with the results shown in Fig. 2(a).The S 21 plots of the 24-cylinder model for both orientations are in good quantitative agreement with the normalized transmission magnitudes in Fig. 2(a) (with the exception of parallel sample 1), indicating that the transmission ratio between both orientations (parallel vs. perpendicular) has been very closely found.Therefore, the use of multiple cylinder layers is a more accurate means of modeling the arrays used in this study.
As discussed previously, the lower THz transmission magnitude of parallel sample 1 when compared to parallel sample 2 is likely due to sensitivity to strand spacing.However, based on the observation described above for the parallel orientation that each cylinder layer decreases the S 21 value by ≈ 1 dB, then it is also possible that sample 1 contains more layers than sample 2.  There is good quantitative agreement between the model with 24 layers with the normalized THz transmission magnitude of parallel sample 2. (e) For the perpendicular orientation, each cylinder layer decreases the S 21 value by ≈ 0.5 dB for frequencies above 0.5 THz.There is good quantitative agreement between the model with 24 layers with the normalized THz transmission magnitude of both perpendicular samples up to ≈ 1.8 THz.

Verifying consistency in periodic models
Although the good agreement between experimental and numerical results is encouraging, its success necessitates a check on whether the numerical results are biased by unexpected inconsistencies in the model and/or simulation.One check is to compare the numerical results obtained from the model in Fig. 4(d) with those generated by other HFSS models that contain structural variations.Three models are presented in Fig. 7(a).The first model (model i), which is similar to Fig. 6(a), contains only one cylinder enclosed inside similar periodic and symmetric boundaries as Fig. 4(d).The second model (model ii) contains two cylinders at the chosen spacing, and therefore considers two periods of model i.The third model (model iii), which is similar to Fig. 4(d), introduces a small offset between the cylinders.This offset breaks the periodicity, and introduces some randomness in the cylinder arrangement.
(a) Periodic models containing (i) only 1 cylinder (similar to Fig. 6(a)), (ii) 2 cylinders with no offset between cylinders, and (iii) 2 cylinders with offset between cylinders (similar to Fig. 4(d)).The grid spacing in these three figures is set to 0.5 μm along both the x and y-axes.
(b) S 21 plots from the various periodic models Fig. 7. (a) Alternative periodic models: model i is the simplest periodic model containing only one cylinder; model ii is equivalent to model i but with two cylinders with no offset between them; model iii, which is the same as that shown in Fig. 4(d), comprises of two cylinders with an offset between them.(b) Comparison between the S 21 plots from: model i, model ii, and model iii.The plots from all 3 models mostly overlap with only slight displacement due to meshing differences between the models.The overlap shows consistency in the HFSS simulations, with no unexpected artifacts arising from the presence of multiple cylinders in the problem space.
Model i is structurally simple as it contains only one cylinder, reducing any errors that may arise from the presence of multiple cylinders in the problem space.Model ii is equivalent to model i, and is included to check the consistency of the periodic computational model.Therefore the S-parameters of models i and ii should be similar with only small differences that arise from meshing.If the presence of multiple cylinders causes any inconsistency in the model, then the two sets of S-parameters will differ noticeably.
The S 21 plots of the single-(model i) and the double-cylinder with no offset (model ii) pe-riodic models overlap when juxtaposed in Fig. 7(b), indicating that both models are indeed equivalent.Furthermore, the model in Fig. 4(d) (model iii) also overlaps the model i and model ii plots.This means that artifacts have not been introduced from the presence of multiple cylinders in the problem space, nor by the offset between cylinders.It can therefore be concluded that the numerical results generated by HFSS are not corrupted by meshing inconsistencies.

Conclusion
In this study, we have presented novel THz-TDS measurements of subwavelength fiberglass arrays, and have utilized a full-wave electromagnetic field solver to explain our measurements.The plotting capabilities of the full-wave electromagnetic field solver has permitted us to visualize the near-field characteristics of the fiber structures under investigation, which is invaluable for understanding the mechanisms and interactions between the fibers, providing an extra dimension of information that would otherwise be lacking in transmission-mode experiments.Furthermore, the use of a commercially available full-wave electromagnetic field solver package, such as HFSS, provides for ease of reproduction of the results.We have also included in this study a performance comparison between analytical solutions for scattering, and a mathematical technique for modeling arrays.The analytical approaches do not include the near-field interactions between objects, whereas the full-wave electromagnetic field solver intrinsically takes this into account as it solves Maxwell's equations.One potential extension of this investigation is to replace cylinders with spheres in order to simulate powders and other spherical scatterers commonly encountered in THz experiments.Both cylinders and spheres can also be used to study the influence of surface roughness in the THz regime, which is of interest in the realm of THz communications and medically-inspired THz applications.These investigations will be the focus of our future work.

Fig. 1 .
Fig. 1. (a-b) Strands of fiberglass are taken from woven fiberglass cloth (plain weave E glass, weight 7.5 oz/yd 2 ), and aligned as a near-periodic array.(c) The array is clamped between two polyethylene frames, with a square aperture of size 20 mm×21 mm.One of the polyethylene frames is covered with aluminum foil to prevent transmission of THz through the frames.The frame is placed orthogonal to the incident THz beam, which is vertically polarized.The red circle with a cross illustrates the incident THz beam into the plane of the fibers.The beamwidth is ≈ 5 mm.(d) The parallel and perpendicular orientations of the arrays with respect to the THz beam.(e) Schematic diagram of an infinite circular cylinder in the parallel or transverse magnetic (TM) orientation, illuminated by a plane wave at θ = 90 • , with propagation vector k normal to the long axis of the cylinder (z-axis).The electric field vector E is parallel to the z-axis.After [1].

( a )Fig. 2 .
Fig. 2. Frequency response and measured optical properties of the arrays described in Fig. 1, in the parallel and perpendicular orientations.(a) Terahertz transmission magnitudes (normalized) of two of the four arrays.The perpendicular arrays have similar magnitudes, unlike the parallel arrays where the magnitudes vary from sample to sample.(b) The measured effective refractive indices n of the fiberglass arrays vary slightly between the two orientations, suggesting geometry-induced birefringence.The error bars account for uncertainties in the samples' thicknesses.(c) The measured extinction coefficients α of sample 1 changes directions beyond 2 THz, raising the question of whether this change is an artifact.
(a) Scattering cross sections C sca, and C sca, ⊥ normalized by 1 m 2 , plotted on a logarithmic scale (b) C sca, and C sca, ⊥ multiplied by AF where N = 100, and normalized by 1 m 2 , plotted on a logarithmic scale

Fig. 3 .
Fig. 3. (a) C sca, and C sca, ⊥ normalized by 1 m 2 , and plotted on a logarithmic y-axis.(b) By multiplying C sca with AF, a scaling effect across all frequencies is observed.
Fig. 4. (a) Overview of the HFSS simulation with the cylinders located about the origin, and the surrounding airbox highlighted in purple.The incident wave propagates in the xdirection.(b) Schematic of the radiation boundaries assigned on the walls of the airbox.The airbox walls, which are parallel to the y-z plane, are λ max /2 away from the cylinders, where the minimum frequency ν min = 0.088 THz.The two airbox walls that lie parallel to the x-z plane are assigned master/slave boundaries in order to model the periodic array that extends to infinity along the y-z plane.In order to model each cylinder's infinite length with respect to its diameter, the walls parallel to the x-y plane are assigned symmetrical electric field E boundaries for the parallel case (symmetrical magnetic field H boundaries for the perpendicular case).(c) One of the two PML boundaries with depth of λ max /6; this depth ensures that fields incident on the PML boundaries are nearly completely absorbed and are not reflected back into the problem space.(d-e) Zoomed view showing two cylinders enclosed inside the airbox for both the parallel and perpendicular cases.(f) Plan view of the double cylinder structure with a 2 μm offset to account for misalignments between adjacent strands in the fiberglass arrays used in the THz experiment.To create a small 1 μm gap between adjacent strands, a 0.5 μm wallspace is introduced between each cylinder and the airbox walls which are parallel to the x-z plane.The master/slave boundaries replicates the double cylinder structure along the y axis.The grid spacing in this specific figure is set to 0.5 μm along both the x and y-axes.
(a) is compared with periodic models containing up to 24 cylinder layers.Two examples of models containing multiple layers are presented in Figs.6(b) and 6(c).
(a) Periodic model with only one cylinder layer (b) Periodic model with two cylinder layers (c) Periodic model with four cylinder layers (d) S 21 plots of periodic models containing between 1 and 24 cylinder layers (parallel orientation).The normalized THz transmission magnitudes from Fig. 2(a) of the two parallel samples are included for comparison.(e) S 21 plots of periodic models containing between 2 and 24 cylinder layers (perpendicular orientation).The normalized THz transmission magnitudes from Fig. 2(a) of the two perpendicular samples are included for comparison.

Fig. 6 .
Fig.6.The grid spacing of the geometries in figures a-c is set to 0.5 μm along both the x and y-axes.(a) Simplest periodic model containing only one cylinder.(b) Periodic model containing two cylinder layers.The gap between the cylinders is equal to twice the wallspace to ensure that gaps are uniform in both the x and y directions.(c) Periodic model containing four cylinder layers.(d) For the parallel orientation, each cylinder layer decreases the S 21 value by ≈ 1 dB for frequencies above 0.5 THz.There is good quantitative agreement between the model with 24 layers with the normalized THz transmission magnitude of parallel sample 2. (e) For the perpendicular orientation, each cylinder layer decreases the S 21 value by ≈ 0.5 dB for frequencies above 0.5 THz.There is good quantitative agreement between the model with 24 layers with the normalized THz transmission magnitude of both perpendicular samples up to ≈ 1.8 THz.