Efficient analysis of deep high-index-contrast gratings under arbitrary illumination

Nojan Motamedi; Amir Shlivinski; Joseph E. Ford; Vitaliy Lomakin

doi:10.1364/OE.23.033472

1- Introduction

Dielectric gratings are used in numerous applications. One example is “imaging fiber bundle”, an array of parallel dielectric fibers, where each fiber is a high index multimode waveguide surrounded by a lower index cladding (Fig. 1). Such fiber bundles (FB) provide a non-imaging transfer of the spatial distribution of light between two surfaces, with applications such as endoscopic imaging [1,2] and image sensing of a curved input surface [3] to name but a few. The ability to efficiently and accurately model such structures is becoming increasingly important. This is especially true as applications emerge with the spatial period of the fiber structure and the image sensor being on the order of a few microns [4–6]. Accurate and efficient modeling of the spatial resolution of the FB is complicated due to the energy transfer between its modes [7–9] as the extended guided modes interfere strongly as the depth of grating increases.

Fig. 1 Schematic of the problem: deep periodic structure illuminated with arbitrary incident beam.

Download Full Size | PDF

Often, the fibers in the bundle can be considered periodic or periodic with multiple fibers per period. The excitation may be a plane wave, but more often the excitation is a beam or a spatially modulated image. When the bundle is considered to be periodic, the computational problem is one having a periodic geometry with an aperiodic excitation. Several methods can be adapted to numerically model such deep gratings [10–13]. One possible approach is to use coupled mode analysis [14,15] but it makes approximations in terms of neglecting the propagating spectrum. Another approach is the Fourier decomposition of the incident field [16], in which the incident field is decomposed into a spectrum of plane waves. For each such plane wave the problem becomes completely periodic. The transmitted diffraction orders can be calculated for each wavenumber using a periodic solver, e.g. Rigorous Coupled Wave Analysis (RCWA) [17,18]. Transmitted output field is then calculated by taking the inverse Fourier transform of the transmitted field for each diffraction order. However, using this approach becomes complicated for very deep gratings. The difficulties arise from the fact that the transmitted diffraction orders have fast phase variations in their wavenumber dependence, which are associated with the large propagation range.

Fourier decomposition of the incident field has been used along with Rigorous Coupled Wave Analysis [16] to analyze the scattering of a confined beam off a periodic structure. In this method, the incident beam is decomposed into a spectrum of plane waves, the transmission coefficients corresponding to these plane waves are calculated using the Floquet theorem via RCWA, and the final transmitted beam is computed from the inverse Fourier integral. However, such an approach results in a large computational cost for very deep gratings. In this paper we propose an alternative method, using the analytic property of Bloch waves propagation in the grating region to compute the transmitted field with a much smaller number of samples in deep gratings. This computationally efficient method requires keeping track of individual Bloch waves created at the input facet of the grating. The method neglects the insignificant multiple internal scattering of Bloch waves in the grating region. The resulting Bloch wave excitation coefficients and the corresponding components in the transmitted field vary slowly as a function of the angle (or wavenumber). Using the slow variation property, we developed a semi-analytical method that computes the transmitted field using FFTs with a much smaller number of samples. This number of samples is nearly independent of the grating depth. The result is a major computational speed up, which enables fast computations of the transmitted fields in very deep gratings. The paper is organized as follows. Section 2 provides problem statement. Section 3 presents the numerical formulation describing all the steps for efficiently computing the transmitted beam. Section 4 presents numerical results supporting the numerical formulation. Section 5 summarizes and discusses the findings.

2. Problem statement

Consider a grating composed of a dielectric material with a complex permittivity that is periodic in the x direction with period Ʌ (Fig. 1). Each period can comprise only a single fiber or multiple fibers, typically with a random position. The grating depth (d) is much greater than the period or the wavelength with typical values of tens of thousand of wavelength. A TE-polarized optical beam with the electric field in the y-direction is incident on the grating from the top. The incident beam is coupled into a set of guided modes of the grating. The set of guided modes propagate through the structure and are coupled into the transmitted beam propagating downwards.

The considered structure is periodic but the exciting beam is aperiodic. The beam, however, can be decomposed into a continuum of plane waves via Fourier transform. For each such plane wave both the excitation and the structure are periodic. For such a periodic problems the Floquet theorem can be used to represent the solution in terms of a discrete set of plane waves. Furthermore, the rigorous coupled wave analysis (RCWA) can be used to find the coefficients of the scattered plane waves, excited (Bloch) waves in the grating, and the transmitted field. This RCWA solution can be further simplified by utilizing the fact the grating is very deep and neglecting the multiple bounces of the grating Bloch waves between the entrance and exit interfaces. As a result, the problem is reduced into three parts: coupling into the grating at the entrance interface, propagating from the entrance to the exit interface, and coupling out of the grating at the exit interface. Once the coefficients of each plane wave corresponding to the incident beam decomposition are found the resulting transmitted beam is calculated via Fourier transform. The following important components need to be addressed to allow efficiently solving this numerical problem: (i) computing the coupling (scattering, reflection, transmission) coefficients, (ii) calculating the Bloch waves, and (iii) computing the Fourier decomposition integral for the transmitted beam with a small number of samples even at very large propagating distances. These aspects are addressed in the next section.

3. Formulation

3.1 Global RCWA formulation

The electric field of the incident beam at the entrance interface can be represented via a continuous Fourier decomposition

E^{i n} (x) = \frac{1}{2 π} \int {\hat{E}}^{i n} (k_{x}) e^{- j k_{x} x} d k_{x},

where

{\hat{E}}^{i n} (k_{x}) = \int E^{i n} (x) e^{j k_{x} x} d k_{x},

represents coefficients of the incident beam’s plane wave spectrum. Once the plane wave spectrum of the incident wave is known the Floquet theorem can be applied to calculate the scattered fields (reflected and transmitted) via diffraction coefficients for the periodic structure. The diffraction coefficients are found by matching the electric and magnetic fields at input and output interfaces using a set of forward and backward Bloch waves for the grating region. We start with a summary of the general RCWA procedure [17,18] and simplify it to account for the decoupling between the entrance and exit interfaces due to the large depth of the structure. The periodic permittivity is first expanded via the Fourier series

ε (x) = \sum_{m} ε_{m} \exp (j \frac{2 π m}{Λ}),

where

ε_{m} = Λ^{- 1} \int_{}^{Λ} ε (x) \exp (- j x 2 π m / Λ) d x

is the Fourier series coefficient of the complex permittivity and Λ is the period of the grating. Using the Floquet theorem the fields in regions I and II are written as

\begin{array}{l} E^{I, y} (x, z) = E^{i n, y} (x, z) + \sum_{i} R_{i} \exp [- j (k_{x i} x - k_{I, z i} z)] \\ E^{I I, y} (x, z) = \sum_{i} T_{i} \exp {- j [k_{x i} x - k_{I I, z i} (z - d)]}, \end{array}

where

R_{i}

is the reflection coefficient in region I and

T_{i}

is the transmission coefficients in region II,

k_{x, i}^{} = k_{x} - i K, K = \frac{2 π}{Λ} .

The wavenumbers

k_{I, z i}

and

k_{I I, z i}

are given by

k_{l, z i} = {\begin{cases} + k_{0} {[n_{l}^{2} - {(k_{x i} / k_{0})}^{2}]}^{1 / 2} k_{0} n_{l} > k_{x i} \\ - j k_{0} {[{(k_{x i} / k_{0})}^{2} - n_{l}^{2}]}^{1 / 2} k_{0} n_{l} < k_{x i} \end{cases}, l = I, I I .

The electric and magnetic fields inside the grating region can be written using a set of forward and backward Bloch waves

\begin{array}{l} E^{g y} = \sum_{i} \sum_{m}^{} w_{i, m} {c_{m}^{+} \exp (- k_{0} q_{m} z) + c_{m}^{-} \exp [k_{0} q_{m} (z - d)]} \exp (- j k_{x i} x) \\ H^{g y} = - j (ε_{0} / μ_{0}) \sum_{i} \sum_{m}^{} v_{i, m} {- c_{m}^{+} \exp (- k_{0} q_{m} z) + c_{m}^{-} \exp [k_{0} q_{m} (z - d)]} \exp (- j k_{x i} x), \end{array}

where

q_{m}

and

w_{i, m}

are the square root of the eigenvalues with a positive real part and the

(i, m)

element of the eigenvector matrix

W

of matrix

A = {(k_{x i} / k_{0})}^{2} - E

for the

T E

polarization.

E

is the Toeplitz permittivity matrix with its

(i, m)

element being

ε_{i - m} .

The quantity

v_{i, m} = q_{m} w_{i, m}

is the

(i, m)

element of the matrix product

V = W Q

, where

Q

is a diagonal matrix with elements

q_{m} .

The coefficients

c_{m}^{+}

and

c_{m}^{-}

are the amplitudes of the forward and backward Bloch waves yet to be determined within the grating region. The summations over the diffraction order index

i

is over an infinite set of positive and negative integers and the summation over Bloch wave index

m

is over an infinite set of positive integers. These infinite summations are truncated to finite numbers for a numerical implementation. Using equation sets (4) and (6) one can match the electric and magnetic fields at input and output interfaces to obtain the following system of equations

\begin{array}{l} δ_{i 0} + R_{i} = \sum_{m}^{} w_{i, m} {c_{m}^{+} + c_{m}^{-} \exp (- k_{0} q_{m} d)} \\ j [δ_{i 0} \cos (θ) n_{I} - (k_{I, z i} / k_{0}) R_{i}] = \sum_{m}^{} v_{i, m} {c_{m}^{+} - c_{m}^{-} \exp (- k_{0} q_{m} d)}, \end{array}

\begin{array}{l} \sum_{m}^{} w_{i, m} {c_{m}^{+} \exp (- k_{0} q_{m} d) + c_{m}^{-}} = T_{i} \\ \sum_{m}^{} v_{i, m} {c_{m}^{+} \exp (- k_{0} q_{m} d) - c_{m}^{-}} = j (k_{I I, z i} / k_{0}) T_{i}, \end{array}

where

δ_{i 0}

is the Kronecker delta function. Equations (7) and (8) can be solved simultaneously to calculate the scattered diffraction orders, which is the complete rigorous solution of the plane wave scattering from the periodic structure known as RCWA. We refer to this complete solution as Global RCWA (G-RCWA) in the rest of the paper.

3.2 Forward RCWA formulation

The coupling between the forward and backward waves can be weak, e.g. due to losses, scattering, or weak reflections at the interfaces, as in the case of deep transmissive gratings or imaging fiber bundles [4–6]. In this case, the equation sets (7) and (8) are decoupled and may be solved separately, which is equivalent to neglecting the coupling between the forward and backward Bloch waves with coefficients $c_{m}^{+}$ and $c_{m}^{-}$ . The result is a two-step solution procedure for calculating the diffraction coefficients, which we refer to as Forward RCWA (F-RCWA). In the first step of F-RCWA, the input single-interface scattering coefficients $c_{m}^{+}$ , are found via the following equation set [Eq. (7) by neglecting $c_{m}^{-}$ ]

\begin{array}{l} δ_{i 0} + R_{i} = \sum_{m = 1}^{} w_{i, m} c_{m}^{+} \\ j [δ_{i 0} \cos (θ) n_{I} - (k_{I I, z i} / k_{0}) R_{i}] = \sum_{m}^{} v_{i, m} c_{m}^{+} . \end{array}

In the second step, the first interface’s transmission coefficients

c_{m}^{+}

are already known and give the amplitudes of the forward Bloch waves, which serve as an excitation for the second interface after propagating through the depth of grating. The scattering problem at the second interface is described by Eq. (8) with the known value of

c_{m}^{+}

, from which the amplitudes of the transmitted diffraction orders

T_{i}

are found for a given

k_{x}

. Once the transmission diffraction coefficients are found the output field may be written as

E^{T} (x) = \sum_{i} \int {\hat{E}}^{i n} (k_{x}) T_{i} (k_{x}) e^{- j k_{x}_{, i} x} d k_{x},

where

T_{i} (k_{x})

is the transmitted diffraction order versus the wavenumber calculated in Eq. (8). The inverse Fourier integral may be evaluated directly via Discrete Fourier Transform (DFT). Fast Fourier Transform (FFT) can be used if the number of samples is chosen as a product of prime numbers. Therefore, one can rewrite Eq. (10) as

E^{T} (x) = \sum_{i} \sum_{n} {\hat{E}}^{i n} (k_{x}) T_{i} (k_{x}) e^{- j k_{x}_{, i} x},

where

k_{x, i} = k_{x} - i K

is the same as before with

k_{x} = n 2 π / L_{x}

and

L_{x}

is the extent of calculations in

x

direction. If approached based on the Nyquist sampling criteria directly, the solution of Eq. (11) requires a very large number of samples for a very deep grating due to the fast phase variations of

T_{i}

as a function of

k_{x}

. These fast variations are due to the propagation factor

\exp (- k_{0} q_{m} d)

in Eq. (8) associated with the imaginary part of

q_{m}

and the long-range propagation of the Bloch waves. The transmission coefficient

T_{i}

is determined by the contribution of all Bloch waves and therefore, it is impossible to smooth its wavenumber dependence, due to its rapid variations and the large number of samples requirement.

3.3 Semi-analytical approach

The effect of the rapid variations associated with the factor $\exp (- k_{0} q_{m} d)$ can be accounted for analytically if $T_{i} (k_{x})$ is written as a sum of partial transmission coefficient contributions corresponding to individual Bloch waves

T_{i} (k_{x}) = \sum_{m} {\tilde{T}}_{i, m} (k_{x}) e^{- k_{0} q_{m} (k_{x}) d} .

Here,

{\tilde{T}}_{i, m}

are the partial transmission coefficients determined by modifying Eq. (8)

\begin{array}{l} w_{i, m} c_{m}^{+} + \sum_{m^{'}} w_{i, m^{'}} {\tilde{c}}_{m, m^{'}}^{-} = {\tilde{T}}_{i, m} (k_{x}) \\ v_{i, m} c_{m}^{+} - \sum_{m^{'}} v_{i, m^{'}} {\tilde{c}}_{m, m^{'}}^{-} = j (k_{I I, z i} / k_{0}) {\tilde{T}}_{i, m} (k_{x}), \end{array}

where

{\tilde{c}}_{m, m^{'}}^{-}

are the modified Bloch wave reflection coefficients of the second interface as a response to the

c_{m}^{+}

Bloch wave excitation. An important feature of Eq. (13) as compared to Eq. (8) is that the set of Eq. (13) does not have the factor

\exp (- k_{0} q_{m} d)

, i.e.

{\tilde{T}}_{i, m}

and

q_{m}

are independent of the grating depth. As a result, the variations of

{\tilde{T}}_{i, m}

and

q_{m}

with respect to the wavenumber

k_{x}

(angle) only depend on the wavelength and the permittivities, and therefore, are on a much slower scale than the variations of

T_{i} (k_{x})

. However propagation of the Bloch waves form the input to the output interface in a deep grating still introduces fast phase variations in the summation terms of Eq. (12) explicitly via the factor

\exp (- k_{0} q_{m} d)

. For this reason the accurate calculation of Eq. (11) would still require dense sampling if evaluated directly and purely numerically. However, one can use the fact that the factor

\exp (- k_{0} q_{m} d)

is known analytically to reduce the integral sampling rate and make it independent of the grating depth. To this end, one can substitute Eq. (12) into Eq. (10) to rewrite the transmitted field as

E^{T} (x) = \sum_{i} \sum_{m} \int {\hat{E}}^{i n} (k_{x}) {\tilde{T}}_{i, m} (k_{x}) e^{- k_{0} q_{m} (k_{x}) d} e^{- j k_{x, i} x} d k_{x} .

The summation in Eq. (14) can be taken only over the propagating Bloch waves while the integrand is negligible for cutoff (decaying) Bloch waves. For propagating modes the integrand has two components, including the amplitude

{\hat{E}}^{i n} (k_{x}) {\tilde{T}}_{i, m} (k_{x})

that is slowly varying function of

k_{x}

and the exponential

\exp (- k_{0} q_{m} d)

that is rapidly varying for large propagation distances. However, in comparison to

\exp (- j k_{x, i} x)

, the exponential power related to

q_{m}

is a much slowly varying function of

k_{x}

. Therefore, one can tabulate the amplitude

{\hat{E}}^{i n} (k_{x}) {\tilde{T}}_{i, m} (k_{x})

and the eigenvalues

q_{m}

at a coarse set of

k_{x}

samples, sufficient to resolve their behavior, e.g. by local interpolation. The tabulated values can further be used to calculate the Fourier integral in Eq. (14). To this end, the Fourier integral in Eq. (14) can be represented via a discrete sum of integrals

E^{T} (x) = \frac{1}{Δ k_{x}} \sum_{i} \sum_{m} \sum_{n} \int_{k_{x_{n}}}^{k_{x_{n + 1}}} {\hat{E}}^{i n} (k_{x}) {\tilde{T}}_{i, m} (k_{x}) e^{ϕ_{i, m, n} (k_{x})} e^{- j k_{x}_{, i} x} d k_{x_{n}},

with

ϕ_{i, m, n} (k_{x}) = - k_{0} q_{m} (k_{x}) d .

First order Taylor expansion of the phase term can be used to evaluate the fast integrand

e^{ϕ_{i, m, n} (k_{x})}

in Eq. (15) analytically.

I_{i, m, n} = \frac{1}{Δ k_{x}} \int_{k_{x n}}^{k_{x_{n + 1}}} e^{ϕ_{i, m, n} (k_{x})} d k_{x} = \frac{1}{Δ k_{x}} \int_{k_{x_{n}}}^{k_{x_{n + 1}}} e^{[ϕ_{i, m, n} (k_{x n}) + {ϕ^{'}}_{i, m, n} (k_{x n}) . (k_{x} - k_{x n})]} d k_{x} = e^{ϕ_{i, m, n} (k_{x n})} \frac{e^{{ϕ^{'}}_{i, m, n} (k_{x n}) Δ k_{x}} - 1}{Δ k_{x} {ϕ^{'}}_{i, m, n} (k_{x n})} .

Here

I_{i, m, n}

is the result of piecewise analytical inverse Fourier integral. This way the fast integrand is computed accurately to the first order between the sampling points of the coarse table. The analytical evaluation of the fast integral also makes the calculations independent of the grating depth or propagation distance d, such that there is no need for taking larger number of samples for deeper gratings. Substituting Eq. (16) into Eq. (15) yields

E^{T} (x) = \sum_{i} \sum_{n} {\hat{E}}^{i n} (k_{x}) e^{- j k_{x}_{, i} x} \sum_{m} {\tilde{T}}_{i, m, n} I_{i, m, n} .

The DFT summation over n in Eq. (17) is the most time consuming part of the calculations. Comparing Eq. (17) with Eq. (11) one should note that the summation over the Bloch waves (m) in Eq. (17) is already taken in Eq. (11). However, the inverse DFT summation in Eq. (11) is taken over a much larger number of samples compared to the same calculation in Eq. (17). The novelty in using Eq. (17) relies in taking the piecewise inverse Fourier transform analytically before summing up the contribution of the Bloch waves as opposed to Eq. (11). In other words, the scattering coefficients

{\tilde{T}}_{i, m}

are weighted by

I_{i, m, n}

before calculating the inverse DFT, which was performed by using the FFT algorithm. Since

I_{i, m, n}

coefficients are calculated analytically from a simple function, the summation

(m)

is computed rapidly as compared to overall computation time.

4- Numerical analysis and results

We start with verifying that the F-RCWA approach compares well to the G-RCWA one. To that end, a straight grating (waveguide array) of 1.5 μm period and 1 cm depth was analyzed at a wavelength of 550 nm with the F-RCWA and G-RCWA solutions. The grating depth corresponded to 18,200 wavelengths and depth-period ratio was 12,133. Schott’s imaging fiber bundle’s (FB) material [19,20] was used to make the comparisons more realistic. The FB consisted of 30% cladding with refractive index of 1.48 and 70% core with refractive index of 1.81.

Figure 2(a) shows the total transmitted power calculated by G-RCWA and F-RCWA versus incident angle when the medium in region II is vacuum. Total transmitted power was calculated by summing up the transmitted diffraction orders. The rapid fluctuations of the transmitted power can be observed in both calculated methods due to the fast variations of the exponential terms describing the Bloch waves propagation through the grating. Figure 2(b) shows the transmitted power when an index matching oil with refractive index of 1.6 was used in region II similar to what was used in the experiment. The index matching oil reduces the amplitude of multiple reflections and the scattering of the field at the output interface. The transmitted beam was also calculated in both non-matched and matched cases as shown in Figs. 2(c) and 2(d), respectively. Although the variations in the transmitted diffraction orders and power are high in the non-matched case, the transmitted intensity profiles are very similar when the sampling is dense enough for both cases. This agreement shows that the F-RCWA approach is sufficiently accurate.

Fig. 2 Transmitted power calculated through G-RCWA and F-RWA when coupling out to (a) vacuum and (b) index matched material with refractive index 1.6; (c),(d) normalized transmitted intensity corresponding to (a) and (b) respectively.

Download Full Size | PDF

We proceed with the analysis of validity of the semi-analytical approach of Sec. 3.3. The same waveguide parameters were used to analyze the F-RCWA solution using both the proposed semi-analytical and direct DFT approaches. Medium II was chosen to be vacuum. Although the negligible multiple reflections are disregarded in the forward solution, the problem of fast phase variations still exists and requires a large number of samples for the direct solution.

Figure 3(a) shows the transmission efficiency of the zeroth diffraction order calculated with F-RCWA. The rapid variations of the individual diffraction orders versus the incident angle are due to the fast phase variations of the propagating Bloch waves. Accurate calculation of the output field would then require a large number of samples in the Fourier domain to be able to follow these rapid variations. However, the phase term $- k_{0} q_{m} d$ and the amplitude term ${\tilde{T}}_{i, m}$ in Eq. (15) are slowly varying functions of $k_{x}$ as seen in Fig. 3(b) and Fig. 3(c), which show the phase and amplitude of the first two Bloch waves. Yet, for a large propagation distance, the exponential summation of the involving phases yields the rapidly varying transmission coefficients as a function of wavenumber (or angle). One can take advantage of the knowledge of these smooth functions by using Eqs. (15) and (16) to analytically calculate the rapidly varying integrand with much smaller number of samples.

Fig. 3 (a) Fast variation of 0_th diffraction order with incident angle calculated with F-RCWA, (b), (c) slowly varying phase and amplitude as functions of angle.

Download Full Size | PDF

Figure 4 compares the semi-analytical and direct FFT approaches using F-RCWA solution as well as the G-RCWA solution for the transmitted output intensity. The same 1 cm long waveguide array is considered and the incident wave is a Gaussian beam with a $1 / e$ beam radius of $0.8 μ m$ and normal angle of incidence. Medium II was index matched to the refractive index of the fiber core. For this case 2,389 Fourier domain $(M k)$ samples are sufficient to represent both the incident field’s spectrum and the variations of the Bloch eigenmodes of the waveguide array. This number of samples is enough to calculate the output field properly using the proposed semi-analytical approach of Eq. (17). However, if the conventional direct inverse (FFT) approach is used one has to increase the number of samples in the Fourier domain from 2,389 to 286,515 (a factor of 120 for the given examples) in order to calculate the transmitted field properly. Using the proposed method one can calculate the output field efficiently and accurately with a much smaller (sparser) number of samples compared to the conventional direct inverse FFT method. One can clearly note the excellent matching between the efficient semi-analytic approach with 2,389 points compared to the direct FFT approach with 286,515 points.

Fig. 4 Transmitted field after 1cm of propagation through the waveguide array calculated with F-RCWA.

Download Full Size | PDF

In order to quantify the efficiency of the proposed method, Normalized Root Mean Square Error (NRMSE) of the output intensity was calculated versus the depth of the same waveguide array through F-RCWA. The error is defined as

N R M S E = \sqrt{\frac{\int_{L} {(I_{1} - I_{2})}^{2} d x}{\int_{L} I_{1}^{2} d x}}

where

I_{1}

is the accurate solution calculated with G-RCWA,

I_{2}

is the solution with semi-analytic and direct F-RCWA methods for smaller number of samples, and L is the integration range. Figure 5 shows the calculated normalized RMSE through direct IFFT and the proposed semi-analytic approach. The error is then calculated with respect to the accurate field in the ± 10 μm range (i.e. L = 20 μm) for all depths of the waveguide array. As the depth of grating is increased a larger number of samples is required to achieve the same level of accuracy. The reason for that, as described before, is the faster variation of the transmitted diffraction orders and propagating Bloch waves with spatial frequency (angle), which requires a larger number of samples to calculate the field properly. For the same number of samples, the proposed semi-analytic approach yields a more accurate result as it calculates the inverse Fourier transform peace-wise analytically.

Fig. 5 Normalized RMSE for F-RCWA calculated with direct IFFT and proposed semi analytic integration.

Download Full Size | PDF

Importantly, the accuracy of the proposed semi-analytical method is to an acceptable level independent of the depth of grating. Therefore, the electric or magnetic field can be calculated at any depth without the need to increase the number of samples as long as the initial number of samples is high enough to smoothly follow the Bloch waves variations with incident angle. The computation time of the proposed semi-analytic method, direct F-RCWA method and direct G-RCWA method on a single core of a 2.6 GHz Intel Core i7-3720QM CPU is presented in Fig. 6. The number of samples was chosen such that the normalized RMSE was fixed at the value of 2% for all grating depths by using Eq. (18) for the three calculation methods. One can note that the computation time for the proposed method is nearly independent of the grating depth. The slight increase of the computational time for deeper gratings in the proposed method is due to the fact that in this case higher order Bloch waves contribute more to the output beam and hence a slightly greater number of samples is required. The proposed method is much faster than both direct G-RCWA and direct F-RCWA methods.

Fig. 6 Computation time for 2% NRMSE error using this work, F-RCWA and G-RCWA methods.

Download Full Size | PDF

Table 1 summarizes the improvement in computation time and number of samples for the given problem. It is evident that the proposed method is up to 33 faster than the G-RCWA and F-RCWA methods for the given grating depths.

Table 1. Improvement in computation parameters compared with G-RCWA

View Table

Finally, Fig. 7 shows the normalized intensity of a 1D periodic waveguide array calculated with the proposed method. The waveguide core and cladding refractive indices are 1.554 and 1.550 and the widths are 3.5 μm and 5 μm respectively. The incident wave is a Gaussian beam with $1 / e$ beam width of 3 μm at a wavelength of 633 nm and normal angle of incidence. We note that that the intensity profile at each height was calculated using the same Bloch eigenmodes and eigenvectors. As the Gaussian beam propagates through the waveguide, it couples into the adjacent waveguides and the diffraction pattern of the waveguide array is observed.

Fig. 7 Intensity pattern of a 1D periodic waveguide array.

Download Full Size | PDF

5- Summary

An efficient numerical method was proposed to calculate the transmission of an electromagnetic beam through a deep periodic dielectric grating. The incident beam was decomposed into its Fourier spectrum of plane waves and the propagating Bloch waves for the periodic grating region were calculated for each plane wave component using RCWA. The RCWA solution was simplified by considering the forward propagating contributions and it was shown that there is an insignificant difference between the F-RCWA and G-RCWA solutions for deep gratings. As the depth of grating increases the variation of the transmission coefficients becomes faster as a function of spatial frequency, and therefore a larger number of samples is required for accuracy. Individual treatment of the propagating Bloch waves enabled us to calculate the inverse Fourier transform semi-analytically using both analytical integration of individual Bloch waves and FFT. The advantage was using the FFT to maintain the speed while accounting for the fast phase variations through analytical integration to maintain the accuracy of calculations with a smaller number of samples. It was shown that the presented formulations lead to accurate and efficient calculation of the output field.

Acknowledgments

This research was supported by the DARPA SCENICC program under contract W911NF-11-C-0210. The research was conducted while A. Shlivinski was on sabbatical leave at the University of California San Diego. The authors would like to thank Dr. Salman Karbasi for invaluable discussions.

References and links

1. O. Coquoz, C. D. Depeursinge, R. Conde, and E. B. de Haller, “Microendoscopic holography with flexible fiber bundle,” Proc. SPIE 2132, 466–474 (1994). [CrossRef]

2. G. Oh, E. Chung, and S. H. Yun, “Optical fibers for high-resolution in vivo microendoscopic fluorescence imaging,” Opt. Fiber Technol. 19(6), 760–771 (2013). [CrossRef]

3. T. S. Axelrod, N. J. Colella, and A. G. Ledebuhr, “The wide-field-of-view camera,” in Energy and Technology Review (Lawrence Livermore National Laboratory, 1988).

4. J. Ford, I. Stamenov, S. J. Olivas, G. Schuster, N. Motamedi, I. P. Agurok, R. Stack, A. Johnson, and R. Morrison, “Fiber-coupled monocentric lens imaging,” in Imaging and Applied Optics Technical Papers, OSA Technical Digest (Optical Society of America, 2013), paper CW4C.2.

5. S. Karbasi, A. Arianpour, N. Motamedi, W. M. Mellette, and J. E. Ford, “Quantitative analysis and temperature-induced variations of moiré pattern in fiber-coupled imaging sensors,” Appl. Opt. 54(17), 5444–5452 (2015). [CrossRef] [PubMed]

6. A. Arianpour, N. Motamedi, I. P. Agurok, and J. E. Ford, “Enhanced signal coupling in wide-field fiber-coupled imagers,” Opt. Express 23(4), 5285–5299 (2015). [CrossRef] [PubMed]

7. K. L. Reichenbach and C. Xu, “Numerical analysis of light propagation in image fibers or coherent fiber bundles,” Opt. Express 15(5), 2151–2165 (2007). [CrossRef] [PubMed]

8. X. Chen, K. L. Reichenbach, and C. Xu, “Experimental and theoretical analysis of core-to-core coupling on fiber bundle imaging,” Opt. Express 16(26), 21598–21607 (2008). [CrossRef] [PubMed]

9. J.-H. Han and J. U. Kang, “Effect of multimodal coupling in imaging micro-endoscopic fiber bundle on optical coherence tomography,” Appl. Phys. B 106(3), 635–643 (2012). [CrossRef] [PubMed]

10. F. Montiel and M. Neviere, “Differential theory of gratings: extension to deep gratings of arbitrary profile and permittivity through the R-matrix propagation algorithm,” J. Opt. Soc. Am. A 11(12), 3241–3250 (1994). [CrossRef]

11. K. Knop, “Rigorous diffraction theory for transmission phase gratings with deep rectangular grooves,” J. Opt. Soc. Am. 68(9), 1206–1210 (1978). [CrossRef]

12. T. Delort and D. Maystre, “Finite-element method for gratings,” J. Opt. Soc. Am. A 10(12), 2592–2601 (1993). [CrossRef]

13. T. Iizuka and C. M. de Sterke, “Corrections to coupled mode theory for deep gratings,” Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 61(44 Pt B), 4491–4499 (2000). [CrossRef] [PubMed]

14. D. Marcuse, Light Transmission Optics (Van Nostrand Reinhold, 1972).

15. A. Yariv, “Coupled-mode theory for guided-wave optics,” IEEE J. Quantum Electron. 9(9), 919–933 (1973). [CrossRef]

16. S.-D. Wu, T. K. Gaylord, E. N. Glytsis, and Y.-M. Wu, “Three-dimensional converging-diverging Gaussian beam diffraction by a volume grating,” J. Opt. Soc. Am. A 22(7), 1293–1303 (2005). [CrossRef] [PubMed]

17. M. G. Moharam, E. B. Grann, D. A. Pommet, and T. K. Gaylord, “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings,” J. Opt. Soc. Am. A 12(5), 1068–1076 (1995). [CrossRef]

18. M. G. Moharam, D. A. Pommet, E. B. Grann, and T. K. Gaylord, “Stable implementation of the rigorous coupled-wave analysis for surface-relief gratings: enhanced transmittance matrix approach,” J. Opt. Soc. Am. A 12(5), 1077–1086 (1995). [CrossRef]

19. “An Introduction to Fiber Optic Imaging” download available at http://www.us.schott.com/lightingimaging/english/download/fo.book.pdf.

20. N. Motamedi, S. Karbasi, J. E. Ford, and V. Lomakin, “Analysis and characterization of high-resolution and high-aspect-ratio imaging fiber bundles,” Appl. Opt. 54(32), 9422–9431 (2015). [CrossRef] [PubMed]

Depth (um)	100	550	1000	5500	10000
Improvement in Computation Time	5x	8.5x	10.7x	30.2x	33x
Improvement in # of Samples	4.1x	7.5x	8.5x	22x	25.1x

Efficient analysis of deep high-index-contrast gratings under arbitrary illumination

Abstract

1- Introduction

2. Problem statement

3. Formulation

3.1 Global RCWA formulation

3.2 Forward RCWA formulation

3.3 Semi-analytical approach

4- Numerical analysis and results

5- Summary

Acknowledgments

References and links

Cited By

Figures (7)

Tables (1)

Equations (18)

Optics Express