1. Introduction

A rapidly growing research interest has been focused lately on the scientific field of terahertz (THz) technology, which deals with the electromagnetic (EM) properties of materials and devices in the so-called THz gap, i.e., the spectral window that links the more extensively investigated microwave and infrared (IR) frequencies. A broad range of novel applications is envisaged, based on material properties that are particular to this spectral region. These infiltrate numerous scientific fields, spanning from chemical detection of explosives or other hazardous materials, safe bio-imaging and detection of diseases, to inherently secure short-range communications and non-destructive testing in industry [1].

For these advances to reach technological maturity, the necessary instrumentation for the generation, manipulation, and detection of THz waves is indispensable. Apart from sources and detectors [2], the extensive functionalities required for the development of high-end systems typically demand for a number of components, such as modulators, switches and filters, whose properties are dynamically tuned by means of an external control signal. In this context, liquid crystalline (LC) materials offer a traditional ultra-low power consumption solution, which has been thus far widely explored and demonstrated in the field of tunable optics, photonics, and plasmonics [3–5]. Their large dielectric anisotropy and their response to applied electric or magnetic fields provide a direct means to dynamically control the refractive and polarization properties of LC-infiltrated elements. These favorable properties are also present in the THz spectrum, stimulating research on both advanced LC materials and technology [6]. Indeed, a whole range of THz tunable components, based on the use of both traditional and novel highly-birefringent and low-loss LC materials [7–11], has been already demonstrated, such as tunable waveplates and polarizers [12,13], filters [14,15], gratings [16], lenses [17], beam-steerers [18], and metamaterials [19, 20].

As it is not always possible to identify analytical solutions, the design and investigation of advanced THz components demands for numerical tools, capable of capturing the EM properties of materials and structures particular to this spectral region. The finite-difference time-domain (FDTD) method is an established and efficient numerical tool for the broadband characterization of EM devices ranging from the microwave to the visible (VIS) frequencies [21] and it has also been successfully implemented for the analysis of advanced THz structures, e.g., metamaterials [22]. When LCs come into play, the FDTD method can be extended in order to include material anisotropy [23] and it can be properly adjusted for the simulation of non-dispersive LC materials [24, 25]. Nevertheless, contrary to the VIS/IR spectrum, where typically a few-μm LC-infiltrated cells and cavities are employed [26, 27], the rigorous investigation LC-THz components needs to take into account both the material absorption losses, which often show strong dichroism, and the dispersion of the LC refractive indices in the broad THz spectrum.

Recently, some of us have introduced an auxiliary differential equation (ADE) FDTD dispersive/anisotropic formulation, which models the LC index dispersion via Lorentzian terms. It was shown that the formulation is capable of describing the LC dispersive properties in the VIS/IR [28]. In this work, we extend the FDTD scheme by modeling the LC dielectric properties via modified Lorentz (m-Lor) terms, a generalized function of dispersion that encompasses as subcases the well-known Debye, Drude, and Lorentz models. The m-Lor model provides extensive degrees of freedom in terms of representing phenomenologically the complex index dispersion of experimentally characterized LC materials, as it was recently demonstrated in the case of isotropic metals and semiconductors [29, 30]. In this work, the complex dispersion of both reference [7, 8] and novel experimentally characterized [9, 31] LC-THz materials is fitted to the m-Lor model, successfully describing their dielectric properties. The proposed formulation is tested and validated in two benchmark examples, transmission through a LC-slab placed between crossed-polarizers and a LC-tunable Fabry-Pérot filter, which have analytical and semi-analytical solutions, respectively. In addition, time-domain spectroscopy (TDS) [1, 32, 33] experimental measurements of THz pulse propagation through LC cells are reproduced numerically, demonstrating the effectiveness of the proposed formulation in the broadband investigation of LC-THz structures.

2. Time-domain modeling of THz liquid-crystals via modified Lorentzian terms

The local molecular orientation of nematic LC materials with an arbitrary spatial alignment profile is described via the nematic director, a unit vector defined via the tilt θ and twist ϕ angles, as shown in Fig. 1(a). The local orientation directly impacts on the LC anisotropic EM properties, which are described via the relative permittivity tensor ε̃ [24, 28]

 

Fig. 1 (a) The local LC orientation is described by the nematic molecular director, which is defined via the spatially dependent tilt (θ) and twist (φ) angles. (b) Unit cell employed in the proposed FDTD formulation.

Download Full Size | PDF

In the following, focusing on the case where the optical axis of the LC lies in the x − y plane (θ = 0), the frequency-dependent non-negative elements of the relative permittivity tensor are given by

In [28] the Lorentz model was chosen for the frequency description of the dielectric tensor elements, which is sufficient for the accurate description of LC index dispersion in the VIS/IR, typically described by the Cauchy model [34]. Nevertheless, in the THz spectrum the LC dispersion depends strongly on the material chemical composition and molecular resonances and, in general, it cannot be described by simple analytical models. Moreover, LC-THz materials show strong dichroism and, contrary to VIS/IR LC-based applications, the absorption losses along both the ordinary and the extraordinary axes must be taken rigorously into account. To this end, the description of ε_o/e(ω) is here modeled via m-Lor terms. It has been shown that the m-Lor analytical function can provide higher quality broadband fitting to experimentally measured permittivity values for many isotropic materials such as silicon [29] and noble metals [35, 36], while it can be proved that it satisfies the Kramers-Kronig relations. Following the proposed formulation, the LC permittivities are expressed as

It can be observed in (6)–(9) that the auxiliary variables W_• and R_• obey the differential equation of a m-Lor medium, involving the dispersive function of either ε_o(ω) or ε_e(ω), which is multiplied with non-dispersive sine/cosine coefficients. The differential equation of W_x1, for instance, is discretized using central finite differences of second-order accuracy

Similar equations can be obtained for the other auxiliary variables. Using the update equations of W_•, a difference equation for $D_x^{n+1}$ can be obtained. The difference equation for $D_y^{n+1}$ derives similarly from the update equations of the R_• variables. Finally, the update expressions of $E_x^{n+1}$ and $E_y^{n+1}$ are calculated by solving the 2×2 system of linear equations that involve D_x and D_y, as described in extent in [28]. It is stressed that the extension of the classical Lorentz model (a₁ = 0) in the proposed implementation comes at no cost in terms of auxiliary variables and, therefore, memory requirements.

3. Numerical results and discussion

In this section, the accuracy of the proposed methodology is validated by investigating LC-THz structures based on three benchmark LC materials, the well known single molecule cyanobiphenyl 5CB [7], the highly birefringent nematic mixture BL037 [8, 37] and the material 1855, an in-house synthetized low-loss mixture of fluorinated biphenyls [9, 38], whose molecular structure is shown in the inset of Fig. 2(b). The dielectric properties of these LCs are described by their complex index of refraction ñ_o/e = n_o/e − jk_o/e, where n_o/e is the refractive index and k_o/e the absorption coefficient. Nematic LC-THz materials are both dispersive and dichroic, which means that the accurate modeling of their properties depends on the fitting of the four functions n_o/e(λ) and k_o/e(λ). Here, using the data provided in [7–9], the experimental values are fitted in the 0.5 – 2 THz window to a m-Lor term via the nonlinear optimization Nelder-Mead technique implemented in the fminsearch function of Matlab^© [39], making use of the known expressions ℜ{ε} = n² − k² and ℑ{ε} = 2nk. The fitted parameter values for the ordinary and extraordinary permittivities of the three LC materials are summarized in Table 1 and a comparison between the experimental data and the fitted functions is shown in Fig. 2. A single m-Lor term is sufficient to capture the dispersive behavior of both the refractive index and the absorption coefficient, despite the fact that their functions show different trends. It is pointed out that the fitted parameters of Table 1 are not necessarily associated with the physical properties of the LC materials, but they provide an analytical description of their permittivity dispersion, in the context of the time-domain studies of LC-THz structures.

Table 1. Fitted m-Lor parameters for the ordinary (o) and extraordinary (e) permittivities of the three LC materials under study in the 0.5 – 2 THz spectral window.

View Table

 

Fig. 2 Experimentally measured and m-Lor fitted values of (a) refractive index and (b) extinction coefficient of the three reference LC materials used in this study, i.e., 5CB [7], BL037 [8], and the low-loss mixture of fluorinated biphenyls 1855 [9,38]. Markers denote experimental values, continuous and dashed lines denote fitted extraordinary and ordinary index values, respectively. The inset in (b) shows the molecular structure of the three constituent molecules of 1855 (n = 2, 3, 5).

Download Full Size | PDF

Figure 3 shows a direct comparison between measured data [9] on the characterization of the mixture 1855 and a FDTD simulation of the experimental setup, using the fitted m-Lor model parameters of Table 1 and the formulation presented in Section 2. The LC was characterized using a standard TDS technique [1]. A 0.5-mm LC layer was sandwiched between two 1.5-mm quartz plates LC cells, modeled with an index n_g = 1.95. Copper wires as separators and electrodes were used, and a homogeneous LC alignment was obtained by applying a high voltage of about 30 kV/m. In this work, the reference signal was transmitted through two directly-adjacent quartz plates, and it was used as the FDTD excitation. The application of a fast Fourier transform (FTT) on the reference and the transmitted signal yields their complex transmission coefficients T(ω), whose amplitude and phase contains the information on the LC complex indices, as described in detail in [31]. The model here proposed accurately reproduces both the retardation, directly linked to the refractive indices n_o/e and the attenuation of the THz pulse, which stems from the material absorption given by k_o/e. The FDTD parameters used were a spatial step Δz = 0.5 μm, and time-step Δt = 0.8Δz/c, where c is the speed of light in vacuum.

 

Fig. 3 Measured THz pulse waveforms for the experimental characterization of the LC mixture 1855 [9] and FDTD simulations using the m-Lor fitted values of Table 1.

Download Full Size | PDF

As a next case example, we study EM wave propagation through a LC layer placed between crossed polarizers as in the inset of Fig. 4. The analytical expression for the transmittance through the layer is given by

 

Fig. 4 Transmittance of a LC slab between cross polarizers calculated by the proposed method for the three reference LC materials. Results are compared with the analytic solution of both the experimentally measured LC index data and their fit to the m-Lor dispersive model. The slab thickness is d = 1.5 mm and the LC twist angle φ = 45°.

Download Full Size | PDF

Finally, the properties of the LC-tunable Fabry-Pérot filter shown in the Fig. 5(a) are investigated via the proposed method and compared to the semi-analytical solution provided by the 4 × 4 Berreman matrix method [41]. The reflective mirrors of the filter are composed of alternating layers of glass and air, whose thickness is matched for a central bandgap frequency of 1 THz, namely d_a = 75 μm and d_g = 38.5 μm, for refractive index values equal to n_a = 1 and n_g = 1.95. The resonant cavity is a 150-μm layer of the LC material 1855, with a uniform alignment in the x − y plane described by the twist angle φ between the director and the x-axis. The simulation parameters used were Δz = 0.5 μm, Δt = 0.8Δz/c, and the computational domain was backed by a 12-cell convolution perfectly matched layer for the efficient absorption of outgoing waves [42]. The FDTD domain was excited by a modulated Gaussian pulse with spectral content in the region from 150 μm to 600 μm.

 

Fig. 5 Transmittance of a LC-THz Fabry-Pérot filter calculated via the Berreman 4 × 4 method and the proposed formulation. The reflective mirror is composed of alternating slabs of glass (d_g = 38.5 μm, n_g = 1.95) and air (d_a = 75 μm, n_a = 1). The resonant cavity of thickness d_LC = 150 μm is infiltrated with the LC mixture 1855, whose director lies in the x − y plane at a fixed angle φ with the x-axis.

Download Full Size | PDF

The filter shows two resonant defect modes within the bandgap at approximately 0.8 and 1.1 THz, whose spectral position depends on the LC twist angle, as in the inset of Fig. 5. When φ = 45° the resonance splits in two defect modes, stemming from the filter’s response to the two orthogonal polarizations of the impinging THz wave, which undergoes polarization rotation while passing through the LC slab. The reference solution is calculated for the experimental index values and it is juxtaposed with the numerical simulation of the fitted LC material. The FDTD model is shown to predict the semi-analytical results with very good accuracy. Again, it has been verified that the small degree of disagreement, found less than 1%, is due to the fitting error of the LC indices. This error can be minimized by using more advanced fitting techniques and/or reducing the noise of the fitted experimental samples. Furthermore, in case the ultra-broadband modeling of LC materials is targeted, where multiple resonances are expected [43], the fitting can be optimized by using additional m-Lor terms, albeit at the cost of increased algorithm complexity and memory requirements to store the additional auxiliary variables.

4. Conclusions

In conclusion, a numerical technique has been presented for the rigorous time-domain analysis of LC-based THz structures, which takes into account the dispersive properties of both the refractive indices and the absorption coefficients of LC materials in the THz spectrum. The fitting parameters for three reference LC materials are provided and numerical studies validate the accuracy of the proposed model in reference LC-tunable THz components. Furthermore, experimental measurements of the time-domain spectroscopy characterization of THz-LC are accurately reproduced, thus demonstrating the efficiency of the presented method as a numerical tool for the broadband investigation of LC-THz devices.