Open Access
Add to BibSonomy
Abstract
We use femtosecond laser-based two-photon polymerization (TPP) to fabricate a 2.5D micropillar array. Using an angular detection setup, we characterize the structure’s scattering properties and compare the results against simulation results obtained from a novel electrodynamics simulation method. The algorithm employs a modified formulation of the Lorentz Oscillator Model and a leapfrog time differentiation to define a 2D coupled Oscillator Finite-Difference Time-Domain (O-FDTD). We validate the model by presenting several simulation examples that cover a wide range of photonic components, such as multi-mode interference splitters, photonic crystals, ring resonators, and Mach-Zehnder interferometers.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Polymeric optical components have become increasingly popular for their low-cost, ease of processing and tunable optical properties such as refractive index thermos-optic effects and birefringence [1]. Most remarkably, when exposed to femtosecond pulsed laser light, certain photosensitive polymers can undergo Two-Photon Absorption (TPA) and a subsequent photopolymerization process. The nonlinearity of the effect confines effect in three dimensions creating an excitation voxel [2]. This is known as Two-Photon Polymerization (TPP) and it has enabled exploring Direct Laser Writing (DLW) of 3D microstructures with applications ranging from biology [3] to photonics [4]. Recent employments of TPP for photonic applications include optical components [4], free-standing optical interconnects for optoelectronic circuits [5–7] and liquid crystal photonic devices [8].
Compared to conventional microfabrication techniques, TPP enables the fabrication of high-resolution low roughness free-standing complex 3D microstructures in a single writing step, with no need for support structures [9]. Its main drawbacks relate to femtosecond laser costs, requirement of suited transparent materials and low large-area fabrication speed compared to shadow-mask or master-based techniques such as conventional photolithography and nanoimprint lithography, respectively. Several parallel writing strategies have been investigated, such as hologram-based multi-focus [10], spatial light-modulator-based multi-focus [11] multimode fiber-base writing [12], ultrafast random-access digital micromirror scanning binary holography [9] and diffraction optical element-based multi-focus [3].
Towards the establishment of TPP as a photonic structure fabrication technique, we study simple yet highly versatile micropillar array structures. In photonics, micro and nanopillar structures have been recently developed for subwavelength Terahertz micropillar lasers [13], waveguide-coupled nanopillar LEDs [14] and quantum-dot micropillar-based neuromorphic computing [15]. Furthermore, polymeric micropillar structures have been developed for photonic crystal-based sensors [16], tunable optical properties of textured surfaces [17] such as structural color or transparency, among others [18].
More general implementations of polymeric pillar arrays include doped polymer arrays for analysis of cellular mechanical responses [19,20], functionalized micropillars for fluid wicking-based DNA detectors [21] and self-cleaning smart surfaces [22].
Photonic simulation algorithms are crucial tools not only for the understanding of the optical properties of such structures but also for application-driven design optimization. Discrete numerical time-domain electrodynamic simulations can be performed using a variety of models such as Finite-Difference Time-Domain (FDTD), Finite-Volume Time-Domain (FVTD), Finite-Element Time-Domain (FETD), Discontinuous Galerkin Time-Domain (DGTD), Finite Difference Frequency-Domain (FDFD) and Pseudospectral Time-Domain (PSTD) [23]. Most of these methods either solve or apply Maxwell’s field equations over time. FDTD is most famous for its accuracy, stability, versatility, and reasonable simulation times. However, there is a constant search for novel ways of modeling electric field-matter interactions optimized for specific systems [24–27]. Here we strive to develop a new modeling approach optimized for dielectric material systems, with reduced calculation complexity, simulation time, and CPU requirements.
In this work, we use femtosecond-based TPP to fabricate a highly ordered micropillar array and experimentally characterize its transmission scattering properties. Also, we develop a novel 2D computational electrodynamics model derived from a modification of the Lorentz Oscillator model, which uses the motion of abstract coupled-oscillators to model the propagation of the electric field waves. We call it Oscillator Finite-Difference Time-Domain (O-FDTD), and it can be visualized as the propagation of a mechanical wave on the surface of a fluid.
In the O-FDTD model, we use a regular square meshing and a leapfrog method to perform spatial and temporal differentiation, respectively. Doing so results in a single global electric field updating equation, which excludes the need for alternating orthogonal electromagnetic field components found in FDTD. Figure 1 illustrates the approach, and the derivation of the model and algorithm implementation can be found in Section 2.
Fig. 1. 3D illustration of an oscillating electric dipole (centered at the origin, $\lambda$=500 nm). The electric field is modelled as a planar network of coupled oscillators whose vertical displacement is proportional to the electric field amplitude.
We apply the method to simulate light interactions with a dielectric micropillar array and compare the results against theory and the experimental characterization of a TPP-fabricated sample. We believe this novel approach holds great potential for the understanding, design and optimization of diverse photonic components and circuits, which we demonstrate by the simulation of well-known and highly relevant structures: The Multimode Interference (MMI) splitter, with important applications in light splitting and coupling [28,29]; Dielectric-based photonic crystals, which can be engineered for light control [30], but can also be found in natural systems [31–33]; The ring resonator [34], which finds applications as varied as sensors e.g. of refractive index [35] and temperature [36], in integrated photonic circuits such as micro-lasers [37], photonic decoders [38], and optical Digital-to-Analog Converters (DAC) for Pulse Amplitude Modulation (PAM) transmitters [39]; And finally, the Mach-Zehnder interferometer (MZI), with applications ranging from gas-refractive index sensors [40], label-free biosensors [41], mechanical strain sensors [42], plasmonic optical modulators [43], to characterization of quantum material properties [44,45].
2. Theoretical model
The description of electric field oscillations through oscillating electrons has been previously addressed by the well-known Lorentz Oscillator Model (LOM) [46,47]. However the LOM can only be applied when electrons are present and thus is not suited for vacuum modelling. In this work a new model is proposed which, similarly to LOM’s approach, considers a set of oscillators whose displacement is proportional to an electric field amplitude. However, such oscillators are not electrons but rather abstract field carriers, also present in vacuum. The oscillators are coupled in a regular rectangular mesh, so that the field propagates similarly to mechanical waves propagating on a water pond surface, see illustration in Fig. 1. The “force” between adjacent oscillators determines the field transfer rate from one oscillator to its neighbors, and consequently, the wave propagation speed. We regard this “force” as the coupling strength between oscillators. Hence, the speed of light in different propagation media (vacuum included), as well as interface effects between different media are governed by a refractive index-dependent coupling strength.
In this section, the motion of isolated oscillators is introduced according to the LOM, followed by the modified formulation for the coupled oscillators.
2.1 Isolated oscillators
The motion of isolated oscillators can be described in an analogy to Hooke’s law for the case of a driven, damped harmonic oscillator, oscillating in the $z$ direction
Equation (3) can be recognized as part of the well-known LOM.
It is important to notice that in vacuum, the density of electrons $N=$ 0, and thus the right side of Eq. (3) becomes zero, rendering it unsuited for electric field quantification in vacuum. In fact, the LOM can be used to describe the dielectric function of certain materials, by defining the expression of $E(t)$, typically as a sine function, and solving Eq. (3) with respect to $\epsilon (\omega )$. However, that is not the scope of this work. Instead, we are interested in the LOM description of the electric field as a damped harmonic oscillator. In the following section, we use an analogy of the LOM to propose a new electron density-independent formulation of the second time-derivative of the electric field, to describe the propagation of electric field waves across a mesh of coupled oscillators.
2.2 Coupled oscillators
In the coupled oscillator O-FDTD approach, the spatial discretization is performed using a square mesh with resolution $\delta x = \lambda /R$, while the temporal discretization is performed by iterations of $\delta t = T/R$ time intervals. $\lambda$ is the vacuum wavelength, $T$ is the oscillation period and $R$ is the number of mesh points per wavelength, which we shall refer to as the resolution coefficient. Using a leapfrog time backward differentiation scheme, the electric field in a given cell at the iteration $i$ can be obtained as follows
Applying the same backward finite difference to the LOM (Eq. (3)), the expression of $\mathrm {d}^2E_{i-1}/\mathrm {d}t^2$ for isolated oscillators can be written as:
The first term of Eq. (6) relates to an external or incident electric field, the second term is equivalent to the restoring spring motion and the third term relates to the energy damping or absorption. Using a similar structure, we propose the following expression for the electric field in a network of coupled oscillators (electric field carriers)
While the LOM describes oscillators with a natural frequency $\omega _0$ excited by an external oscillating electric field $E^{ext}$, the O-FDTD model Eq. (7) describes a mesh of oscillators with no explicit resonance frequency, though which an electric field wave flows at a rate defined by $\Gamma$. Light-absorption is accounted by the damping coefficient $\gamma$. The neighbourhood electric field, $E^{neigh}$, is the average field of the eight surrounding oscillators, which can be efficiently calculated using kernel matrix convolution.
using the kernel matrixWe define the coupling strength as
where $n$ is the real part of the refractive index and $C(\Omega )$ is a constant that depends on the size of the kernel matrix $\Omega$. The term $n^2$ modulates the light speed in the medium and for the $\Omega$ specified in Eq. (9), $C(\Omega )=8/3$. Hence, $\Gamma$ defines the proportionality between the temporal oscillation frequency and the wavelength and thus the speed of light and wavelength in different media. An empirical verification of this formulation can be found in the Supplement 1 (SD). In Fig. S1 it is shown that the correct wavelength and propagation speed is achieved if $C(\Omega )=8/3$, while any other value will lead to incorrect results. In Fig. S2, the normalized simulation of an oscillating electric dipole is compared against the FDTD simulation and the analytical solution. In Fig. S3, the normalized simulation of a plane wave incident on an absorbent material is compared against the FDTD simulation and the analytical solution. Furthermore, Fig. S10 shows the simulation of a plane wave incident on a slab of material, as well as the interference-induced reflectance spectra of dielectric-coated Si substrates.We define the damping coefficient $\gamma$ as
where $c$ is the speed of light, $\kappa$ is the imaginary part of the refractive index, $\lambda$ is the vacuum wavelength and $\gamma$ is the absorption coefficient that accounts for the exponential decrease of the light intensity inside a light-absorbent material. Again, $n$ modulates the speed of light and the wavelength inside the material, while $\kappa$ is quantifies the light absorption.The motion Eq. (7) can thus be re-formulated as
Together, Eqs. (4), (5) and (12) describe the complete O-FDTD model.
3. Materials and methods
3.1 Custom O-FDTD simulation algorithm
An implementation of the O-FDTD model into a simulation algorithm has been developed to facilitate structure design and real-time result visualization. Light sources are implemented by introducing a sinusoidal electric field oscillation in a defined region of the simulation space. For the boundary conditions, we match the effect of Perfectly Matched Layers (PML) by assigning a linear gradient of absorbent layers around the simulation space to prevent reflections at the boundary interface [48,49]. periodic boundary conditions (PBC) are implemented by a geometrical inversion of the boundaries [50]. Light sources, boundary conditions and parameter-varying simulation times are described in detail in Section 1 of the SD. Further, the algorithm was implemented into a standalone MATLAB toolbox which is provided as an executable. A typical screenshot of the simulation interface is shown in Fig. S7 of the SD. The O-FDTD simulation results presented in this manuscript were obtained using a resolution coefficient of R = 40 points per wavelength.
3.2 3D microfabrication setup
A Newport $\mathrm {\mu }$FAB writing station with a fixed focused laser beam is used for fabrication. Sample translation stages in X-Y have a traveling range of 100 mm with 1 nm resolution and minimum step size of 10 nm while in the Z-direction a traveling range of 4.8 mm with 20 nm resolution and minimum step size of 60 nm. The presented nano-pillar array is fabricated with a 40x Nikon dry objective microscope (NA 0.75). Sample focusing and fabrication process was monitored via a mounted CCD camera. As light source for TPP, a femtosecond pulsed laser (Tsunami, Spectra Physics) was used and tuned to a wavelength of 795 nm and a pulse length of approximately 90 fs, after a prism compressor. For power control a half-wave plate and polarizer are inserted into the beam path. The fabricated structure was designed using $\mu$FAB, Newport.
3.3 Sample preparation
Cleaned 170 $\mathrm {\mu }$m-thick glass microscope coverslips are rinsed first in ethanol, then water and heated at 100 °C on a hotplate for 30 minutes to evaporate solvent residues. 20 $\mathrm {\mu }$L of a low-shrinkage hybrid organic-inorganic zirconium containing sol-gel polymer SZ2080 is drop-casted on the glass substrate and subsequently baked for 30 minutes at 100 °C. After baking the sample is mounted in the $\mathrm {\mu }$FAB writing station and the designs are traced in the polymer layer. Development of the fabricated structures is performed for 45 min in a 1:2 solution of 4-methyl-2-pentanone and 2-propanol with subsequently air drying.
3.4 Structure design and TPP writing parameters
The fabricated structure consists of a 120x120 pillar square, with a pitch of 5 $\mathrm {\mu }$m, resulting in a 600x600 $\mathrm {\mu }$m array. The fabrication speed was set to 25 $\mathrm {\mu }$m/s and the laser power 12 mW before the back aperture of the microscope objective. A more detailed description of the TPP setup and experimental parameters used for TPP microfabrication in SZ2080 polymer can be found in [3].
3.5 Spectral-angular characterization setup
The spectral-angular diffraction profiles are characterized using a fiber-coupled Xe lamp white light source, focused on the sample with a 4x microscope objective (PLN4X, Olympus). The sample was fixed at the center of a 300 mm rotating breadboard (RBB300A/M, Thorlabs) and a fiber-coupled spectrometer (AvaSpec-ULS2048x64-EVO, Avantes) was used to measure the transmission spectrum as function of angle (angular resolution $\delta \theta =$ 0.5 °).
4. Results
4.1 Experimental results
A uniform periodic 2.5D array of micropillars is fabricated using femtosecond laser-based TPP DLW. Figure 2(a) shows a 3D illustration of the TPP-fabricated structure with optical microscope images as insets. The laying pillar (right inset) was obtained by pushing it with a contact profilometer tip and from the side-profile of the pillar, a height and width of about 5 $\mathrm {\mu }$m and 1 $\mathrm {\mu }$m, respectively, could be extracted. Optical microscope images of the complete fabricated array can be found in Fig. S8 of the SD.
Fig. 2. Light scattering of TPP-fabricated micropillar array. a) 3D schematic of the fabricated sample (120x120 pillars, 5 $\mathrm {\mu }$m spacing). Insets show microscopic images of the array unit cell and a knocked over pillar. b) Experimental spectral-angular light scattering profile. The inset shows a photograph of the scattering pattern displayed on a white screen 20 cm away from the sample. c) Simulated time-integrated squared electric field $E^2$ map around a micropillar. d) Experimental, simulated and calculated diffraction intensity profile for a wavelength of 650 nm.
The scattering properties of the structures are experimentally characterized. Figure 2(b) shows the complete spectral-angular profile in a 2D logarithmic color-scaled map. The black lines trace the position of the diffraction maxima calculated from Fraunhofer diffraction theory [51]. The inset image shows a photograph of the diffraction pattern produced on a white screen positioned 20 cm away from the sample. In the following, we compare the experimental results with O-FDTD simulations.
4.2 Simulation results
The light-scattering properties of TPP-fabricated micropillars are simulated using the O-FDTD model via the implemented 2020 WaveBox software for a section of the periodic arrays. In order to do that, the 2.5D pillar system was reduced to 2D, by considering a simulation plane cutting vertically through the center of the pillars. A planar wavefront (line source) emitting with a vacuum wavelength of $\lambda =650$ nm is placed inside the glass substrate. A refractive index of 1.46 + i0.017 ($\lambda$ = 650 nm) is considered for the SZ2080 polymer. The boundary conditions are PBC in x and PML in y. Figure 2(c) shows the resulting time-integrated $E^2$ map, where the white and red dashed lines indicate the edge of the pillars and light source, respectively. Figure 2(d) show the far-field angular projection profile. The simulation results are compared with Fraunhofer diffraction calculations and a selected line from the experimental results corresponding to 650 nm [red line from Fig. 2(b)]. Details regarding the Fraunhofer diffraction calculations and comparison with simulations for different slit conditions can be found in Section 2B of the SD.
Agreement between simulations, experimental and simulated results is observed in terms of diffraction order angle and peak Full Width at Half Maximum (FWHM) up to the 7th diffraction order, extending over 3 orders of magnitude of light intensity.
4.3 Applications of O-FDTD for other key photonic structures and devices
The agreement obtained between experimental, theoretical and simulation results has prompted us to explore the applicability/advantages of the developed O-FDTD method towards some of the most relevant and widely used photonic building blocks. Here, the applicability and versatility of O-FDTD is showcased via simulations of an MMI splitter, a photonic crystal waveguide, a ring resonator and a Mach-Zehnder interferometer. These structures were selected to highlight the possibility to simulate micro-waveguiding, photonic bandgap and light resonance effects using O-FDTD, thus creating opportunities for many other relevant photonic structures and devices. Supplementary results on light-matter interactions at material interfaces are presented in Fig. S10 of the SD. The supplemental video, Visualization 1 shows a simulation example using the 2020 WaveBox graphical user interface. The executable MATLAB toolbox files are available in supplemental Code 1 [52] for download.
We present a simple 1 $\times$ 4 rectangular MMI design composed of a Si waveguide core surrounded by a SiO2 cladding. A splitter length of 9 $\mathrm {\mu }$m was used to split an input signal into four equal outputs. The time-integrated squared electric field $E^2$ and real electric field $E$ simulation results can be found in Fig. 3(a). The dashed lines indicate the Si/SiO2 interface. By comparing the input power with the sum of output power, we calculate a coupling efficiency $\eta _{coupling}=$ 90$\%$. Optimized splitter cavity geometries and output/output tapering can be used to reduce reflections and improve the splitting efficiency [28].
Fig. 3. Photonic building blocks. a) Time-integrated $E^2$ and time-frame $E$ of a Si/SiO2 core/shell 1x4 Multimode Interference (MMI). The dashed lines indicate the Si/SiO2 interface. $\lambda$ = 1550 nm. b) $E$ map of a photonic crystal waveguide. $\lambda$ = 1550 nm. c-e) Ring resonator. Light is pumped through the in waveguide at wavelength c) 1945 nm and d) 1970 nm. Geometry parameters: $w=$ 200 nm, $S=$ 130 nm, $R=$ 3 $\mathrm {\mu }$m. e) Resonance peak at out depending on the spacing $S$. f) Broadband spectrum showing resonance periodicity.
Photonic bandgap effects in periodic subwavelength structures is demonstrated using a photonic crystal waveguide. The designed photonic crystal consists of a 2D array of Si pillars (circles) with radius of $r=$ 0.158 $\times a$, surrounded by air. The separation between closest neighbors is given by the lattice constant $a$, with an angle of 60$^o$ between layers. The photonic bandgap was calculated using MIT Photonic-Bands [30,53], see bandgap diagram in Fig. S11 of the SD. The chosen frequency for the bandgap was $\nu =$ 0.41 c/a. For a wavelength of 1550 $\mathrm {\mu }$m, the lattice constant of the crystal is $a=\nu \times \lambda =$ 0.6355 $\mathrm {\mu }$m and the pillar radius becomes 0.1004 $\mathrm {\mu }$m. Figure 3(b) shows the O-FDTD-simulated $E$ map of a light beam propagating in the designed “z-shaped” photonic crystal waveguide. Two sharp waveguide bends of 120$^o$, which would not be possible in a “conventional” waveguide.
The effect of light resonance is illustrated by the simulation of a double-waveguide ring resonator, see Figs. 3(c), 3(d). We use Si/SiO2 core/cladding waveguides with width $w=$ 350 nm and ring radius $R=$ 3 $\mathrm {\mu }$m. For most wavelengths, the phase shift is expected to generate destructive interference inside the ring, causing most of the light to be transmitted to the through (top-right) waveguide. For specific resonance wavelengths, constructive interference inside the ring leads to the out-coupling of light via the out (bottom-left) waveguide. Figures 3(c) and 3(d) show the time-integrated $E^2$ and $E$ amplitude maps for input vacuum wavelengths of 1945 nm and 1970 nm, corresponding to off-resonance and on-resonance, respectively. Pin, Pthrough and Pout are the measured powers at the in, through and out waveguides, respectively. Figure 3(e) shows the intensity spectra at the through and out waveguides, $T_{through}=P_{through}/P_{in}$ and $T_{out}=P_{out}/P_{in}$, respectively, for a resonance around 1970 nm and spacing $S$ ranging from 50 to 170 $\mathrm {\mu }$m. Shorter waveguide spacing results in stronger coupling to the ring but broadened resonance peaks, as reported elsewhere [34]. We choose a spacing of $S=$ 120 nm and simulate the transmittance spectrum from 1600 nm to 2600 nm - Fig. 3(d). Due to rather small selected ring radius, the resonance peaks are rather spectrally far, since the periodicity of the resonances scales inversely with the ring radius. Thus, due to the wide wavelength range selected, it is possible to observe the effect of wavelength-varying coupling efficiency into the ring, which changes the spectral response of the resonances.
The applicability of resonance effects to the sensing, for instance, of refractive-index-based biological samples is illustrated by the simulation of an MZI. Figure 4(a) shows the design for a Si/SiO2 MZI. In this design, a semi-triangular taper is used to couple an incoming light beam (Source waveguide) into a 20 $\times$ 4 $\mathrm {\mu }$m$^2$, 1 $\times$ 2 MMI beam-splitter. In the upper waveguide, a sensing arm is designed with a length of 3.35 $\mathrm {\mu }$m. A symmetric structure is used to re-couple the two beams into the Output waveguide.
Fig. 4. Mach-Zehnder Interferometer. a) Schematic of the Si/SiO2 design. b) Time-integrated $E^2$ maps for effective refractive indices of $n=$ 3.47 and $n=$ 3.25. c) Transmittance as function of the sensing arm effective refractive index $n$ for $L=$ 3.34 and 6.7 $\mathrm {\mu }$m. The black lines for simulated, color for calculated results.
An effective refractive variation at the sensing arm will induce a phase shift, which in turn creates destructive/constructive interference response at the deice output. Figure 4(b) shows the time-integrated $E^2$ maps of a 300 fs light pulse (vacuum wavelength $\lambda =$ 1.5 $\mathrm {\mu }$m) for effective refractive indices of 3.47 and 3.25 at the sensing arm. The refractive index change induces a relative light intensity at the out waveguide variation from 0.777 to 0.028. The phase delay depends on the sensing arm’s length and refractive index according to $T_{Out}\propto \mathrm {cos}(\frac {2\pi L}{\lambda }(n - n_{eff}) )$ where $T_{Out}=P_{Out}/P_{in}$ is the transmittance at the Output waveguide, $n$ is the refractive index of the waveguide, $n_{eff}$ is the effective refractive index of the sensing arm and $L$ is its length. Figure 4(c) shows the simulated transmittance (black markers) as function of $n_{eff}$ in comparison with the theoretical curve for sensing arm lengths of 3.35 and 6.7 $\mathrm {\mu }$m. We find a very good agreement between the simulated transmittance and the theoretical curve, with a transmittance $T_{out}$ oscillation between 2.8% and 78%.
5. Discussion
To test the O-FDTD method, we perform comparative simulation and calculation studies using FDTD (Lumerical) simulations and the Fraunhofer diffraction theory. An empirical demonstration of the O-FDTD formulation can be found in the SD, Figs. S1, S2, S3 and S10. The agreement obtained in the comparison of the O-FDTD results against the FDTD simulations and analytical results, can only be achieved by a correct reproduction of the propagation speed, wavelength and surface interactions. In particular, any other value of $C(\Omega )$ other than 8/3 leads to an incorrect speed of light, in which case the propagation results stray from the analytical solution, and phase-dependent interference effects become distorted. Figure S3 demonstrates that the imaginary part of the refractive index is correctly being taken into account in terms of material absorption, since it perfectly matches the Beer Lambert law. Furthermore, the reflectance spectra shown in Fig. S10 demonstrate that both the refractive index contrast at material interfaces, as defined by the coupling strength $\Gamma (\Omega )$, and the interference effects occurring inside the dielectric material are being modelled correctly.
As an experimental test of the proposed O-FDTD model, we use TPP-fabricated 2.5D array uniformly distributed in the xy plane. This symmetry allows reducing the 2D pillar array to a 1D grating problem, since the the angular response of the 2D grating [inset of Fig. 2(b)] can be obtained by the superposition of the 1D grating along the perpendicular directions. However, this approximation may not be possible for more complex, less symmetric 3D architectures.
We present the versatility of O-FDTD via the simulation of a variety of photonic components, which cover a wide range of potential applications. The design parameters were chosen arbitrarily as a proof-of-concept and are not optimized for any specific application. All simulations were performed with a mesh size of 20 points per wavelength unless mentioned otherwise. A convergence coefficient is used to enforce a time step shorter than the analogous spatial resolution by a factor of $\sqrt {2}$, which improves the simulation accurateness and stability. Most time-domain electrodynamics simulation algorithms, such as FDTD, operate by solving Maxwell’s equations over time and space or, like FETD, by applying the solutions of those equations to isolated mesh elements [23]. While algorithms such as FDTD, FDFD, and PSTD use rectangular meshes, which lead to stair-stepping edges in round geometries, approaches like FETD, FVTD, and DGTD make use of unstructured polygon meshes, which adapt much better to the round, irregular, and sharp-angled geometries [54]. Irregular or unstructured meshes can be very appealing for representation and solution accurateness in advanced modeling applications [55]. However, defining such meshes can be very challenging, especially outside of commercial software, and recent works have expended extensive effort in algorithms for mesh optimization [56,57]. Also, FETD methods are often hindered by subdomain nonlocality issues that lead to longer time and higher CPU requirements to enforce continuity between subdomains [23].
Other works focus on hybrid solutions that combine different methods and mesh systems [27,48]. Alternative meshing schemes include the ones used in Nyquist theorem-limited PSTD [58] and meshless FDTD [27] algorithms. The former promise reduced computational requirements but suffer from long-time instability, while the latter ones allow Lagrangian simulations at the cost of longer simulation times and computation complexity. Frequency-domain algorithms such as FDFD provide steady-state solutions rather more efficiently than time-domain methods but are limited to narrowband analysis and require the individual computation of each frequency [23]. A novel and more exotic approach uses electronic circuit elements to model electromagnetic field interactions [26].
The O-FDTD algorithm developed in this work stands out from previous approaches because it is derived from an analogy with a solution of the Lorentz Oscillator Model. Hence, a single global calculation provides the electric field update for the entire simulation space, via leapfrog time backward differentiation, independently from the magnetic field or polarization constraints. Several proof-of-concept examples demonstrate the agreement between our simulations and both theoretical and experimental results. Also, the developed algorithm is compatible with parallel multi-wavelength solving. Thus, an efficient implementation of the algorithm holds the potential to reduce the memory and CPU requirements considerably, compared to normal FDTD. The remarkable algorithmic simplicity and ease of implementation of O-FDTD in any programming platform makes it useful for many scientists and engineers that may lack the resources for commercial tools and the expertise to implement highly specialized approaches. Our implementation of O-FDTD in MATLAB includes a real-time monitoring feature for instant result inspection. The algorithmic simplicity of O-FDTD owes partly to the use of a first order finite difference approximation, the use of a single material parameter (refractive index), the limitation to TM polarization, and the use of a regular square mesh, which, like conventional FDTD, makes it vulnerable to stair-stepping edge issues. It should also be noted that the proposed formulation does not take the material conductivity, i.e. free electrons, into account, which limits its application when simulating metals. Without the conductivity, the proposed O-FDTD method can still reproduce certain properties of metals, such as the absorbance and the refractive index contrast-dependent reflectance, but it fails to describe electric charge effects, such as electric current induction or plasmonic effects. Nevertheless, as demonstrated by several examples in this article, many dielectric and semiconductor simulation applications can benefit from its use. Future developments may include the use of higher-order time finite differences, and a modified mesh size-dependent coupling strength to allow variable-size rectangular meshing for multi-scale problems.
Well-established FDTD and other computational electrodynamic modeling software tools consider other material properties such as electric conductivity. On the other hand, the numerous materials parameters, simulation space settings, light sources, and monitors used by commercial solutions may lead to nonphysical or unrealistic results if not correctly set, for instance, by inexperienced users. Also, the possibility to obtain instantaneous feedback over the simulation results, as opposed to the usual long waiting time required to obtain the final calculation only, holds a considerable practical advantage of our implementation compared to most commercial software. Future work on an improved O-FDTD model may include electric field polarization and magnetic field calculation. Also, in our implementation of the O-FDTD, we define the PML and PBC boundary conditions using geometrical operations. While this can be advantageous in terms of simulation convergence and homogeneity, faster solutions may be achievable using customized boundary condition formulations.
6. Conclusions
We use femtosecond laser-based TPP to fabricate a highly ordered and homogeneous 2.5D array of polymeric micropillars and characterize its light-scattering properties. This type of structure finds many applications, such as in microlensing, light outcoupling, surface reflectivity and transmission tuning, to mention a few examples.
We develop a novel O-FDTD electrodynamics model for the propagation of electric field waves from a modified formulation of the Lorentz Oscillator Mode. We use leapfrog time differentiation and obtain a single global field-equation formulation suitable for parallel wavelength calculations. An O-FDTD simulation algorithm is described and made available to the readers in the form of a software package 2020 WaveBox. The simple modus operandi of our method can be of advantage for less experienced users, including students or researchers from diversified scientific backgrounds.
As examples, we present the simulations of several photonic structures that deploy interference, waveguide, and resonance effects, using the 2020 WaveBox software. We obtain excellent agreement between simulation, theoretical, and experimental results. We envision a wide range of scientific and technological applications of O-FDTD spanning from integrated photonic circuits to applications in biophotonics. The model is scalable to other wavelengths and thus holds the potential for applications, for example, in terahertz and radiofrequency technologies.
Due to the single-equation nature of the model, we foresee the suitability for machine learning-assisted optimization strategies, for example, for photonic component and circuit optimization.
Funding
H2020 Future and Emerging Technologies (828841); Interreg (0181_NANOEATERS_1_EP); H2020 Marie Skłodowska-Curie Actions (713640); European Regional Development Fund (NORTE-01-0145-FEDER-000019).
Acknowledgment
Ricardo Adão acknowledges the PhD programme "Laser, Photonics and Vision" of the University of Vigo, Spain. This research contributes to the goals associated to the PhD project.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
Supplemental document
See Supplement 1 for supporting content.
References
1. M. C. Oh, W. S. Chu, J. S. Shin, J. W. Kim, K. J. Kim, J. K. Seo, H. K. Lee, Y. O. Noh, and H. J. Lee, “Polymeric optical waveguide devices exploiting special properties of polymer materials,” Opt. Commun. 362, 3–12 (2016). [CrossRef]
2. S. Maruo, O. Nakamura, and S. Kawaia, “Three-dimensional microfabrication with two-photon-absorbed photopolymerization Shoji,” Opt. Lett. 22(2), 132–134 (1997). [CrossRef]
3. C. Maibohm, O. F. Silvestre, J. Borme, M. Sinou, K. Heggarty, and J. B. Nieder, “Multi-beam two-photon polymerization for fast large area 3D periodic structure fabrication for bioapplications,” Sci. Rep. 10(1), 8740 (2020). [CrossRef]
4. T. Gissibl, S. Thiele, A. Herkommer, and H. Giessen, “Two-photon direct laser writing of ultracompact multi-lens objectives,” Nat. Photonics 10(8), 554–560 (2016). [CrossRef]
5. N. Lindenmann, G. Balthasar, D. Hillerkuss, R. Schmogrow, M. Jordan, J. Leuthold, W. Freude, and C. Koos, “Photonic wire bonding: a novel concept for chip-scale interconnects,” Opt. Express 20(16), 17667 (2012). [CrossRef]
6. R. Woods, S. Feldbacher, D. Zidar, G. Langer, V. Satzinger, G. Schmid, W. Leeb, and W. Kern, “Development and characterization of optoelectronic circuit boards produced by two-photon polymerization using a polysiloxane containing acrylate functional groups,” Appl. Opt. 52(3), 388–393 (2013). [CrossRef]
7. Z. Gu, T. Amemiya, A. Ishikawa, T. Hiratani, J. Suzuki, N. Nishiyama, T. Tanaka, and S. Arai, “Optical transmission between III-V chips on Si using photonic wire bonding,” Opt. Express 23(17), 22394 (2015). [CrossRef]
8. Z. He, G. Tan, D. Chanda, and S.-T. Wu, “Novel liquid crystal photonic devices enabled by two-photon polymerization,” Opt. Express 27(8), 11472 (2019). [CrossRef]
9. Q. Geng, D. Wang, P. Chen, and S. C. Chen, “Ultrafast multi-focus 3-D nano-fabrication based on two-photon polymerization,” Nat. Commun. 10(1), 2179 (2019). [CrossRef]
10. G. Vizsnyiczai, L. Kelemen, and P. Ormos, “Holographic multi-focus 3D two-photon polymerization with real-time calculated holograms,” Opt. Express 22(20), 24217 (2014). [CrossRef]
11. L. Yang, A. El-Tamer, U. Hinze, J. Li, Y. Hu, W. Huang, J. Chu, and B. N. Chichkov, “Parallel direct laser writing of micro-optical and photonic structures using spatial light modulator,” Optics and Lasers in Engineering 70, 26–32 (2015). [CrossRef]
12. E. E. Morales-Delgado, L. Urio, D. B. Conkey, N. Stasio, D. Psaltis, and C. Moser, “Three-dimensional microfabrication through a multimode optical fiber,” Opt. Express 25(6), 7031 (2017). [CrossRef]
13. M. Krall, M. Brandstetter, C. Deutsch, H. Detz, A. M. Andrews, W. Schrenk, G. Strasser, and K. Unterrainer, “Subwavelength micropillar array terahertz lasers,” Opt. Express 22(1), 274 (2014). [CrossRef]
14. V. Dolores-Calzadilla, B. Romeira, F. Pagliano, S. Birindelli, A. Higuera-Rodriguez, P. J. Van Veldhoven, M. K. Smit, A. Fiore, and D. Heiss, “Waveguide-coupled nanopillar metal-cavity light-emitting diodes on silicon,” Nat. Commun. 8(1), 14323 (2017). [CrossRef]
15. T. Heuser, J. Grose, S. Holzinger, M. M. Sommer, and S. Reitzenstein, “Development of Highly Homogenous Quantum Dot Micropillar Arrays for Optical Reservoir Computing,” IEEE J. Sel. Top. Quantum Electron. 26(1), 1–9 (2020). [CrossRef]
16. X. Gan, H. Clevenson, and D. Englund, “Polymer Photonic Crystal Nanocavity for Precision Strain Sensing,” ACS Photonics 4(7), 1591–1594 (2017). [CrossRef]
17. F. H. Chen, S. Pathreeker, S. Biria, and I. D. Hosein, “Synthesis of Micropillar Arrays via Photopolymerization: An in Situ Study of Light-Induced Formation, Growth Kinetics, and the Influence of Oxygen Inhibition,” Macromolecules 50(15), 5767–5778 (2017). [CrossRef]
18. E. Lee and S. Yang, “Bio-inspired responsive polymer pillar arrays,” MRS Commun. 5(2), 97–114 (2015). [CrossRef]
19. K. Nagayama, T. Inoue, Y. Hamada, and T. Matsumoto, “A novel patterned magnetic micropillar array substrate for analysis of cellular mechanical responses,” Journal of Biomechanics 65, 194–202 (2017). [CrossRef]
20. M. L. Rathod, N. Pareek, S. Agrawal, S. Jaddivada, D. W. Lee, and N. Gundiah, “Engineered ridge and micropillar array detectors to quantify the directional migration of fibroblasts,” RSC Adv. 7(81), 51436–51443 (2017). [CrossRef]
21. M. Geissler, L. Malic, K. J. Morton, L. Clime, J. Daoud, J. A. Hernández-Castro, N. Corneau, B. W. Blais, and T. Veres, “Polymer Micropillar Arrays for Colorimetric DNA Detection,” Anal. Chem. 92(11), 7738–7745 (2020). [CrossRef]
22. H. Ko, Z. Zhang, K. Takei, and A. Javey, “Hierarchical polymer micropillar arrays decorated with ZnO nanowires,” Nanotechnology 21(29), 295305 (2010). [CrossRef]
23. U. S. Inan and R. A. Marshal, Numerical Electromagnetics (Cambridge University Press, New York, 2011), 1st ed.
24. K. Luo, Y. Yi, Y. Duan, B. Xu, and B. Chen, “2D Efficient Unconditionally Stable Meshless FDTD Algorithm,” Int. J. Antennas Propag. 2016, 1–6 (2016). [CrossRef]
25. G. S. Cheng, Z. H. Fan, D. Z. Ding, and R. S. Chen, “An efficient high order plane wave time domain algorithm for transient electromagnetic scattering analysis,” Eng. Analysis with Boundary Elem. 85, 13–19 (2017). [CrossRef]
26. T. M. Minteer, “Mdeling the Interaction of Electrodynamic Fieldsto Circuit Elements,” IETE J. Educ. 59(2), 70–84 (2018). [CrossRef]
27. Q. Sun, R. Zhang, Q. Zhan, and Q. H. Liu, “3D Implicit-Explicit Hybrid Finite Difference/Spectral Element/Finite Element Time Domain Method Without A Buffer Zone,” IEEE Trans. Antennas Propag. 67(8), 5469–5476 (2019). [CrossRef]
28. L. Nikolaevsky, T. Shchori, and D. Malka, “Modeling a 1 x 8 MMI green light power splitter based on gallium-nitride slot waveguide structure,” IEEE Photonics Technol. Lett. 30(8), 720–723 (2018). [CrossRef]
29. R. Reese, H. Tesmer, M. Jost, E. Polat, M. Nickel, R. Jakoby, and H. Maune, “A Compact Two-dimensional Power Divider for a Dielectric Rod Antenna Array Based on Multimode Interference,” J. Infrared, Millimeter, Terahertz Waves 39(12), 1185–1202 (2018). [CrossRef]
30. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic Crystals: Molding the flow of light (Princeton University Press, Woodstock, Oxfordshire, 2011), 2nd ed.
31. L. H. Wu and X. Hu, “Scheme for achieving a topological photonic crystal by using dielectric material,” Phys. Rev. Lett. 114(22), 223901 (2015). [CrossRef]
32. J. Teyssier, S. V. Saenko, D. Van Der Marel, and M. C. Milinkovitch, “Photonic crystals cause active colour change in chameleons,” Nat. Commun. 6(1), 6368 (2015). [CrossRef]
33. M. Lopez-Garcia, N. Masters, H. E. O’Brien, J. Lennon, G. Atkinson, M. J. Cryan, R. Oulton, and H. M. Whitney, “Light-induced dynamic structural color by intracellular 3D photonic crystals in brown algae,” Sci. Adv. 4(4), eaan8917 (2018). [CrossRef]
34. W. Bogaerts, P. de Heyn, T. van Vaerenbergh, K. de Vos, S. Kumar Selvaraja, T. Claes, P. Dumon, P. Bienstman, D. van Thourhout, and R. Baets, “Silicon microring resonators,” Laser Photonics Rev. 6(1), 47–73 (2012). [CrossRef]
35. V. Donzella, A. Sherwali, J. Flueckiger, S. M. Grist, S. T. Fard, and L. Chrostowski, “Design and fabrication of SOI micro-ring resonators based on sub-wavelength grating waveguides,” Opt. Express 23(4), 4791–4803 (2015). [CrossRef]
36. H.-T. Kim and M. Yu, “Cascaded ring resonator-based temperature sensor with simultaneously enhanced sensitivity and range,” Opt. Express 24(9), 9501–9510 (2016). [CrossRef]
37. B. Stern, X. Ji, A. Dutt, and M. Lipson, “Compact narrow-linewidth integrated laser based on a low-loss silicon nitride ring resonator,” Opt. Lett. 42(21), 4541–4544 (2017). [CrossRef]
38. F. Mehdizadeh, M. Soroosh, and H. Alipour-Banaei, “A novel proposal for optical decoder switch based on photonic crystal ring resonators,” Opt. Quantum Electron. 48(1), 20 (2016). [CrossRef]
39. S. Moazeni, S. Lin, M. Wade, L. Alloatti, R. J. Ram, M. Popovic, and V. Stojanovic, “A 40-Gb/s PAM-4 Transmitter Based on a Ring-Resonator Optical DAC in 45-nm SOI CMOS,” IEEE J. Solid-State Circuits 52(12), 3503–3516 (2017). [CrossRef]
40. N. L. P. Andrews, R. Ross, D. Munzke, C. van Hoorn, A. Brzezinski, J. A. Barnes, O. Reich, and H.-P. Loock, “In-fiber Mach-Zehnder interferometer for gas refractive index measurements based on a hollow-core photonic crystal fiber,” Opt. Express 24(13), 14086–1499 (2016). [CrossRef]
41. K. Qin, S. Hu, S. T. Retterer, I. I. Kravchenko, and S. M. Weiss, “Slow light Mach–Zehnder interferometer as label-free biosensor with scalable sensitivity,” Opt. Lett. 41(4), 753–756 (2016). [CrossRef]
42. X. Dong, H. Du, X. Sun, Z. Luo, and J. Duan, “A novel strain sensor with large measurement range based on all fiber mach-zehnder interferometer,” Sensors 18(5), 1549 (2018). [CrossRef]
43. C. Haffner, W. Heni, Y. Fedoryshyn, J. Niegemann, A. Melikyan, D. L. Elder, B. Baeuerle, Y. Salamin, A. Josten, U. Koch, C. Hoessbacher, F. Ducry, L. Juchli, A. Emboras, D. Hillerkuss, M. Kohl, L. R. Dalton, C. Hafner, and J. Leuthold, “All-plasmonic Mach-Zehnder modulator enabling optical high-speed communication at the microscale,” Nat. Photonics 9(8), 525–528 (2015). [CrossRef]
44. C. M. Wilkes, X. Qiang, J. Wang, R. Santagati, S. Paesani, X. Zhou, D. A. B. Miller, G. D. Marshall, M. G. Thompson, and J. L. O. Brien, “Mach – Zehnder interferometer,” Opt. Lett. 41(22), 5318–5321 (2016). [CrossRef]
45. E. G. Idrisov, I. P. Levkivskyi, and E. V. Sukhorukov, “Dephasing in a Mach-Zehnder Interferometer by an Ohmic Contact,” Phys. Rev. Lett. 121(2), 026802 (2018). [CrossRef]
46. J. Li, “Finite element study of the Lorentz model in metamaterials,” Comput. Methods Appl. Mech. Eng. 200(5-8), 626–637 (2011). [CrossRef]
47. A. F. J. Levi, “The Lorentz oscillator model,” in Essential classical mechanics for device physics (Morgan and Claypoll Publishers, San Rafael, 2016), chap. 5, p. 101, 1st ed.
48. A. Fathi, B. Poursartip, and L. F. Kallivokas, “Time-domain hybrid formulations for wave simulations in three-dimensional PML-truncated heterogeneous media,” Int. J. for Numer. Methods Eng. 101(3), 165–198 (2015). [CrossRef]
49. Z. Y. Huang, L. H. Shi, Z. Sun, Q. Lei, and J. B. Wang, “Implementation of Associated Hermite FDTD Method to Periodic Structures,” IEEE Antennas Wirel. Propag. Lett. 16, 2696–2699 (2017). [CrossRef]
50. A. A. Jamali, M. R. Amirzada, M. R. Anjum, and K. Kanwar, “Implementation of Periodic Boundary Conditions in FDTD Algorithm,” J. Inf. Commun. Technol. Robotic Appl. 9, 1–3 (2018). [CrossRef]
51. E. Hecht, Optics (Addison Wesley, San Francisco, 2002), 4th ed.
52. R. Adão, M. Caño-Garciá, J. B. Nieder, and C. Maibohm, “WaveBox Software for Electric Field Propagation Simulations,” figshare (2021), https://doi.org/10.6084/m9.figshare.13359494.
53. S. G. Johnson and J. D. Joannopoulos, “MIT Photonic-Bands,” (2015).
54. A. Gansen, M. E. Hachemi, S. Belouettar, O. Hassan, and K. Morgan, “EM modelling of arbitrary shaped anisotropic dielectric objects using an efficient 3D leapfrog scheme on unstructured meshes,” Comput. Mech. 58(3), 441–455 (2016). [CrossRef]
55. D. A. Fadeev, “Numerical simulation of 2D electrodynamic problems with unstructured triangular meshes,” Comput. Opt. 43(3), 385–390 (2019). [CrossRef]
56. H. Fahs, “High-order discontinuous Galerkin method for time-domain electromagnetics on non-conforming hybrid meshes,” Math. Comput. Simul. 107, 134–156 (2015). [CrossRef]
57. A. Gansen, M. El Hachemi, S. Belouettar, O. Hassan, and K. Morgan, A 3D Unstructured Mesh FDTD Scheme for EM Modelling, 0123456789 (Springer Netherlands, 2020).
58. M. Yeung and E. Barouch, “Long-Time Instability of Pseudospectral Time-Domain Method in Curvilinear Coordinates,” IEEE Trans. Antennas Propag. 68(4), 2993–3001 (2020). [CrossRef]
