1. Introduction

Dielectric and metallic artificial structured media have allowed the tailoring of light propagation with unprecedented freedom in the conception of photonic devices. The engineering of photonic band structures, in analogy with the electronic ones, allows to express a plethora of different behaviors and functions. Advanced control of the wavefront and orbital properties of light, transformation optics, photonic crystal membranes [1,2], cavities with ultra-high Q [3,4], and slow-propagating light waveguides [5,6] that enable strong light-matter interaction for sensing and nonlinearities excitation are just a few examples of the possibilities that this paradigm enables.

Traditionally, waveguiding in dielectric photonic crystals (PhCs) can be achieved through the introduction of defects [7], which localize states within the photonic band gap. These defects will be embedded immutably within the periodic structure and will result in guiding occurring in a very specific region of the artificial crystal. A different approach is based on the engineering of spatial dispersion which can enable self-collimation (SC) of light and which can lead to self-guided beams that can cross each other without giving rise to energy exchange. This mechanism enables light guiding that relies neither on bandgap, nor on the presence of defects in the periodic structure, nor on total internal reflection, but on the exploitation of flat Iso-Frequency-Curves (IFCs) that allows beam-propagation without deformation. Indeed, when a beam with an arbitrary shape impinges on a given structure, as long as its angular spectral content lies within the flat portion of the IFC of interest, it propagates undistorted, by conserving its shape and size. Conversely, any angular component of the beam falling outside of this flat region will propagate according to the local curvature of the IFC. Ultimately, the flat region of an IFC can act as a filter for the angular components of the incident beam: those within the flat region will propagate unperturbed, the rest being filtered out as they propagate. Another interesting property of these structures is that they allow guiding for multiple light beams within the same structure without cross-coupling between wavepaths.

Self-collimation (SC) has been widely investigated in long PhCs [8]. More recently, the existence of an analogous propagation regime based on pure spatial dispersion, termed mesoscopic self-collimation (MSC), has been demonstrated also in mesoscopic photonic crystals (MPhC) [9]. An MPhC is an artificial crystal whose elementary cell is organized in such a way as to allow compensation of spatial dispersion in certain directions and frequencies. In its simplest configuration, an MPhCs appears as alternating slabs of natural and/or artificial material with opposite spatial dispersion. For example, as presented in [9], it may consist of alternating layers of homogeneous material with layers consisting of PhCs (see Fig. 1 (a)). In particular, in [9], it was demonstrated that MSC can be achieved along the direction of high symmetry for the MPhC, i.e., the direction normal to the interfaces between PhC and bulk slabs. This peculiar propagation regime occurs when the curvature of the dispersion band averages to zero across the overall mesoscopic period.

 

Fig. 1. Mesoscopic structure: (a) Principle and periodicity; (b) phase and (c) energy propagation in a single mesoscopic period.

Download Full Size | PDF

In [10], we proposed an hybrid numerical-analytical design method that conjugates MSC condition to the control of impedance matching between the alternating layers composing the MPhC. In particular, this method allows to find the optimal geometrical parameters needed to achieve MSC with tailored reflectivity. Although it was limited to high symmetry directions for the crystal, this simple method allows to easily design MPhC in which light confinement and total transmission or total reflection (both close to 100 %) are simultaneously obtained. This method relied on perfectly satisfying the MSC condition and on an approximate control of the interference conditions at the interface to obtain a given overall reflectivity. In particular, mesoscopic mirrors, which are conceived by combing the requirements of MSC and high reflectivity, were shown to offer much larger angular acceptance than traditional Bragg mirrors, which is especially true along the direction of high symmetry for the MPhC. In fact, along such direction, it is possible to find rather large regions, in the space of the wavevector, where the curvatures of the different layers compensate each other. Thanks to this particular and important feature, it was possible to numerically [11] and then experimentally [12] demonstrate complete confinement of a stable optical mode within a 1D Fabry-Pérot cavity bounded by a pair of mesoscopic mirrors. Such work revealed that, barring radiation losses in the remaining spatial dimension, the proposed 2D models were sufficiently robust and accurate to allow prediction of the behavior of a real, three-dimensional device. However, radiation losses, due to the leaky nature of the modes along the orthogonal direction to the design plane, have been shown to be significant enough to limit the design of certain real devices. In particular, whilst the propagation length in mirrors was low enough that the diffraction losses did not prevent their experimental demonstration [6], they are strong enough that any demonstration of propagation over a significant length could not be achieved.

Differently from [10], in [13] we proposed a numerical method dedicated to the analysis of arbitrary MPhCs that relaxes the high symmetry constraint. In particular, an extensive numerical study was performed using this method. This allowed us to reach a twofold conclusion: (i) it is indeed possible to collimate light even along some directions that are not of high symmetry for the crystal. These directions can be found by searching the whole band diagram of the elementary cell of a MPhC and looking for locally flat isofrequency curves (IFCs); (ii) some of these solutions lie below the light-cone, i.e., in the region of purely guided modes, which do not suffer from the aforementioned radiation leakage. However, whilst the proposed technique allows to find MSC directions in a complex band-structure by an a-posteriori numerical analysis, it is not a model that allows design and parametric study of MSC below the light cone, at oblique incidence.

In this manuscript we propose a novel hybrid numerical-analytical method dedicated to the design of MPhCs at oblique incidence and under the light cone. This has the following threefold goal: (i) ensuring the simultaneous coexistence of MSC and total reflectivity control (and thus limiting or eliminating insertion losses due to impedance mismatch); (ii) including among the solutions those involving a direction of propagation at an oblique angle with respect to the direction of high symmetry for the crystal; (iii) limiting the solutions to those that exclude all diffraction losses, either out-of-plane losses, previously discussed, and in-plane scattering losses at the PhC interfaces.

This design method has the advantage of an extreme analytical simplicity and, differently from the analysis method proposed in [13], allows direct design of oblique incidence MPhC structures rather than discovering them in existing structures through a-posteriori numerical estimation of high-order derivatives of the dispersion surface, which were affected by computational noise [13]. Such a simple yet powerful design tool may pave the way to the conception of advanced photonic interconnections and circuits based on MSC, which are expected to be immune to crosstalk and out-of-plane losses.

The paper is organised as follow. In the following section, the proposed design method will be presented. The quantities required to illustrate the problem and method will be introduced and referred to a generic MPhC. The necessary conditions to ensure the MSC regime, those to control the impedance matching of the multilayer, and those necessary to avoid in-plane diffraction and out-of-plane radiation losses will be detailed. Then, in Section 3, through a practical example, it will be explained how to search for solutions that satisfy all the previously introduced constraints. Both a closed-form design equation and a figure of merit for the assessment of self-collimation efficiency will also be presented. Finally, in Section 4, the design obtained through the proposed method is tested by FDTD simulations as the incident beam width varies. The results are compared with the broadening of a beam of the same lateral size propagating in a homogeneous medium. The impact of the residual mean curvature is also be discussed. To conclude, we investigate the impact of varying the length of the homogeneous layers around the optimal value.

2. Design algorithm description and methodology

The structure supporting mesoscopic self-collimation under oblique incidence is schematically depicted on Fig. 1 (a). A narrow beam of light is impinging under oblique incidence on a periodic structure of period $D = d_1 + d_2$ along the $x$ direction, alternating slabs of bulk material of length $d_1$ and PhC slabs of length $d_2$. Each PhC slab consists of $N$ rows of holes and their length extends beyond the centre of the first and last row (see Fig. 1 (b)). The total length of each PhC slab is therefore $d_2=N$, in units of the lattice parameter $a$. At the frequency of design, the natural lateral energy spreading of the beam in the bulk slabs (beam divergence) is compensated for in the PhC slabs by a self-focusing effect and the beam propagates over many periods without net spreading, providing so called MSC, mesocopic equivalent of SC in a PhC.

The incident beam comes from bulk material with an incident angle $\theta _{1}$ and then propagates inside the structure, at an angle $\Theta$ (see Fig. 1 (a)). Inside the structure the wavevectors are described by their angle $\theta _i$ with the $x$ direction, with $i=1$ in the bulk slabs and $i=2$ in the PhC slabs (see, Fig. 1 (b)). The k-vector components are $k_{xi}$ along $x$, and, since all interfaces are along the $y$ axis, $k_y$ is a natural invariant of the problem, and $k_{y1}=k_{y2}=k_{y}$ (see Fig. 1 (b)). As for the energy propagation, it occurs over the distances $L_i$ along the Poynting vector directions described by the angle $\theta _{\pi i}$ (see Fig. 1 (c)). Note that in the bulk material, $\theta _1 = \theta _{\pi 1}$ as these slabs are isotropic. For the following computations, we use dimensionless units (such as those used in MEEP [14]): lengths are normalized to $a$, representing the period of the square lattice PhC, and the reciprocal vectors are normalized to $2\pi /a$.

Designing an effective AR-MSC structure (i.e. an impedance matched structure in which the overall beam spreading is negligible) involves simultaneously satisfying both anti-reflection (AR) and MSC conditions. Those two conditions can be expressed with simple mathematical formulas, also in the case of arbitrarily oblique incidence by generalizing the equations reported in [10].

The Anti-Reflection conditions in both slabs can be written in dimensionless units as:

The Mesoscopic Self-Collimation condition can be written (in dimensionless units) as:

Equation (3) generalizes the formula reported in [10] to the case of arbitrary angle of incidence. In fact, in [10], the energy propagation was only considered under normal incidence ($\theta _1 = 0$) with respect to the interfaces of the MPhC alternating slabs and the propagation lengths $L_1$ and $L_2$ were replaced by the geometrical lengths $d_1$ and $d_2$. This normal incidence corresponds to the high-symmetry condition for the PhC slab and thus to a propagation regime that does not give rise to anisotropy. Indeed, in this case, the energy and the phase are propagating along the same direction (the $x$-direction). Conversely, in the current case of arbitrarily oblique incidence, anisotropy must be taken into account since energy (Poynting vector) and phase (incident $k-$vector) are not propagating in the same direction. More specifically, in the PhC slab, the energy flows along the normal to the IFCs.

Equation (3) stems from the generalization, to oblique incidence, of the propagator described in [9], which integrates the phases of the multiple plane-waves composing a beam along its propagation path. In our case, the phases are integrated over a length $L_i$ defined by the Poynting vector direction $\theta _{\pi i}$ and by the slab length $d_i$ :

 

Fig. 2. (a) Iso-frequency curves of the infinite square lattice PhC similar to the PhC slabs in the MPhC ($n_1=2.85$, $r/a=0.28$). (b) local curvature index $n_{c2}$ limited to a $[-5$:$+5]$ range. The white line depicts perfect self collimation ($1/n_{c2} = 0$). The coloured vertical segments correspond to the solutions of Eq. (1) in the $d_2=7[a]$ case studied in section 3. Dashed lines correspond to general solutions, continuous lines to the ones with a negative curvature index.

Download Full Size | PDF

Following the notations in Eqs. (1) and (3), our algorithm is rather straightforward. For a given PhC slab length $d_2$ (containing an integer number N of row of holes), we will :

In the following section, we discuss in more details the implementation of this algorithm on a practical example.

3. Practical implementation

For a practical implementation, we focus on an MSC structure composed of a square lattice of air holes (Fig. 1) with a radius $r=0.28\,[a]$ in a medium of optical index $n_1 = 2.85$ which corresponds to the effective index of the guided mode in a $270$-nm-thick GaAs membrane suspended in air at the wavelength $1.55\,\mu$m. In this section, we choose to consider a PhC slab constituted of $N=7$ rows of holes (and thus $d_2 = 7\,[a]$). The first step of our implementation is the calculation of the band diagram for the PhC used in the PhC-slab. The band diagram is calculated for an infinite PhC, modelled by simulating a single period surrounded by periodic boundary conditions.

3.1 PhC slab characteristics

PhC band diagram is calculated with MPB open source software [17]. As the most computationally expensive part of our algorithm, the whole first irreducible Brillouin zone was computed once and exploited in the following calculations. It is worth stressing that, since we are considering an infinitely periodic PhC, the unit cell is small (one single hole in a $a \times a$ square area) and, as we restrict ourselves to the first band of the PhC, the computation time is really short, less than 15 minutes using a standard desktop computer. For the calculations, we use a two dimensional grid of $500\times 500$ k-points and a spatial resolution of 32 pixels per unit $a$. Figure 2 (a) shows the calculated IFCs on the first band for the perfect PhC (solid gray curves).

As explained above, to maximize the transmission through the structure, we restrict our observations to $k_{x2}=p/(2d_2)=p/(2N)$ values satisfying the AR condition in the PhC slab (Eq. (1)). For $N=7$, this condition is satisfied along the 6 vertical dashed lines numbered $p=1$ to $p=6$ superimposed in Fig. 2. Furthermore, we also restrict the computations to points where the curvature index is negative, which is a necessary condition to fulfil the MSC condition (Eq. (3)). Thus, in Fig. 2 (b), the points satisfying both conditions simultaneously are those contained within the short solid segments. Only these points are therefore under investigation.

Each of these 1D segments of potential solutions corresponds to a monotonic relation between $k_y$ and $u$. Plotting the points $k_y(u)$ for each solution results in a set of parametric curves as shown in Fig. 3 (a). The gray-shaded areas in Fig. 3 identify the regions where the interface diffraction condition (hatched gray area) and the light-cone condition (solid gray area) are not satisfied, thus giving, respectively, in-plane diffraction losses and out-of-plane radiation losses.

 

Fig. 3. Parametric curves for the set of solutions: (a) transverse wave-vector $k_y$, (b) curvature index $1/n_{c2}$ and (c) angle of propagation of energy in the PhC slabs $\theta _{\pi 2}$ as a function of $u$. Color curves correspond to specific values for $p$ AR order, as in Fig. 2. The solid gray area corresponds to the light cone (air losses) and the hatched gray area to diffraction by the PhC interface (diffraction losses). The solid thick curves are AR solutions without light-cone or diffraction losses.

Download Full Size | PDF

In order to ensure MSC, we need to calculate the accumulated curvature in the PhC slab (i.e. $L_2/n_{c2}$) which compensates the one in the bulk slab (i.e. $L_1/n_{c1}$), according to Eq. (3). For this purpose, Fig. 3 (b) shows the curvature index $n_{c2}$ calculated from the band diagram of the PhC as a function of the frequency $u$, for the six sets of solutions highlighted in Fig. 2 (b). These curves are obtained from the local curvature of the considered IFC [9].

In addition, Fig. 3 (c) shows, as a function of the frequency $u$, the angle $\theta _{\pi 2}$ between the $x$-direction and the direction of the energy flow (along the normal to the IFC). $\theta _{\pi 2}$ is needed to calculate the propagation length $L_2$ and the MSC condition (Eq. (3)).

As discussed above, we can further reduce our solution set by discarding $k_y$ values corresponding to a propagation above the light cone (solid gray area on Fig. 3 (a)) or that can diffract at the bulk-PhC interface (hatched gray area on Fig. 3 (a)).

The remaining parameter space (thick color curves in Fig. 3) corresponds to potential MSC solutions where energy should flow forward in the structure without any in-plane or out-of-plane diffraction losses.

3.2 Derivation of the bulk slab length

For the sake of clarity, we now focus on the solution $p=6$, but similar analysis can be carried out for the other $p$ solutions. For a given $p$, once the suitable solution space for $k_y(u)$ is defined as described above, we can determine the corresponding values of both $k_{x1}(u)$ and $\theta _1(u)=\theta _{\pi 1}(u)$ (being the energy/phase propagation direction in the bulk slab, defined with respect to the $x$-direction, as a function of the reduced frequency).

The last parameter of the meso-structure that needs to be fixed is the length $d_1$ of the bulk slab. According to the proposed algorithm, we obtain $d_1$ as a solution of Eq. (1) (AR-condition) and then we find the best approximate solution for the MSC-condition Eq. (3). We thus first express Eq. (3) as a function of $d_1$ and $d_2$ using Eqs. (4) and (5):

Figure 4 (a) shows the values of $\theta _1$ as a function of the reduced frequency $u$ calculated according to Eq. (6) for the $p=6$ branch. Moreover, Fig. 4 (b) shows the length of the bulk $d_1$ as a function of $u$ calculated according to Eq. (13). We can see that $d_1$ exhibits a quasi staircase variation: each plateau corresponds to a given integer value for $m$ (the order of the AR condition in the bulk slabs) as shown in Fig. 4 (b). Within each plateau a small variation with $u$ can be observed: it is due to the term $1/[2k_{x1}(u)]$ in Eq. (10) or Eq. (13). Of course the total accumulated curvature over one meso-period $D$ is not strictly null for all the solutions as we choose to prioritize the AR condition which introduces the rounding in Eq. (12) and results in approximated MSC. In order to evaluate the quality of the self-collimation in our structures, we define a figure of merit (FOM) as :

 

Fig. 4. Parametric curves of $\theta _1$, $d_1$ and FOM (see text for definition) as a function of $u$ for the $p=6$ branch of solutions. In (b) each almost constant plateau corresponds to a given order $m$ for the AR condition in the bulk slabs. In (c) three different solutions, labelled $S_1$, $S_2$, and $S_3$, are indicated with black dot markers. The gray regions in (a), (b) and (c) corresponds to the AR order $m=8$.

Download Full Size | PDF

Moreover, we can clearly see the effect of the rounding in Eq. (12): for each constant value of $m$ (and thus each plateau in Fig. 4 (b)), the FOM goes from positive to negative values. As an example, we highlighted the AR order $m=8$ in Fig. 4 (gray area). For this AR order, we can identify three different approximated MSC regimes: $S_1$ near the centre of the plateau, where the $FOM=-2\times 10^{-4}\simeq 0$ and the MSC is almost perfect; $S_2$ towards then end of the plateau where the FOM is negative ($\simeq -2\%$) and the overall MPhC slightly focussing and $S_3$ towards the beginning of the plateau where the FOM is positive ($\simeq +2\%$) and the MPhC slightly diverging.

A similar analysis can be carried out for every value of $m$. As a consequence, despite our initial choice to focus on perfect AR and only consider MSC as an afterthought, for each pair of $(p,m)$, that is each AR solution order for PhC and bulk slabs, we can find a specific reduced frequency, near the centre of each plateau on $d_1(u)$, where the MSC will be perfectly respected.

However, as will be detailed in section 4.2, in practice for devices of finite size, even a FOM of a few percent results in satisfactory collimation.

It is worth pointing out that prioritizing AR-condition is a further element of novelty with respect to the approach in our previous work [10], where, with the propagation being restricted along the x-direction, we first derived $d_1$ from the MSC-condition (Eq. (3)) and then we searched for the reduced frequency $u$ that simultaneously approximate the AR-condition (Eq. (1)), with a desired degree of precision. This change is aimed at ensuring the best energy propagation by minimizing the reflections. Indeed, MSC Eq. (3) only holds for a forward propagating beam that goes only once through each slab. However, when the AR is approximated like in the previous method, due the high number of interfaces in the structure, a significant part of the energy is propagating multiple times back and forth between the interfaces, and the MSC Eq. (3) no longer holds. This leads to unwanted and difficult to predict beam spreading that directly depends on the quality of the AR approximation.

With our new method, we limit these reflections and accept a small residual curvature that will lead to a small and predictable beam spreading that can be negligible for finite-size devices.

4. FDTD validations of the designs

4.1 MSC efficiency of the structures

We arbitrarily decided to study extensively a structure chosen in the middle of the $p=6$ branch (Design $D1$ on Fig. 4). The parameters from our model are $\theta _1 = 27.16^\circ$, $d_1 = 7.542 \,[a]$ (corresponding to a $m=8$ order AR structure in the bulk slabs). The predicted MSC occurs at $u=0.2091\left [ a/\lambda \right ]$ with a quasi-null $\textrm {FOM} = -2 \times 10^{-4}$.

All the simulations presented in this section are 2D-FDTD simulations carried out using MEEP [14]. We have chosen to validate the design method through 2D-FDTD simulations since the 2D model are considerably less expensive in terms of computational resources (i.e. computation time, RAM, etc.) than 3D-FDTD ones, but they are yet as accurate in predicting the actual behavior of fabricated devices. As an example, in [13,18] the authors have shown that 2D-FDTD simulations predict well the spectral features, such as the transmittance behavior and the resonant peaks, of strongly resonant devices. It is worth pointing out that, when present, out-of-plane losses can not be directly evaluated by the 2D model, but require 3D simulations [12]. However, the design method proposed in this paper is intended to perform a systematic search for optimal MSC structures, which are all below the light-line. This condition minimizes the losses by vertical outcoupling in the air superstrate, thus making 2D simulations even more accurate.

We first perform simulations on a $200 \,[a]$-long MSC structure ($\simeq 14$ mesoperiods) with a spatial resolution of 16 pixels$/a$. The source is launched along the computational cell horizontal direction $X$, and the entire MSC structure is tilted so that the angle of incidence on the first interface is $\theta _1$. We use a pulsed source centered at the predicted frequency and with a spectral width of $\delta u =0.1$. The beam profile is Gaussian along the $Y$-direction, with various waist sizes. The field over the simulation cell is saved every $dt = 1/(6 \delta u) = 1.65625$ rounded to the nearest FDTD time-step interval, and afterward frequency-resolved with a sub-sampled time-Fourier-analysis. The method used for this analysis is the same as in [10], and allows to compute with a single simulation the MPhC properties (waist evolution, reflectivity and transmission) and the beam propagation trace at several frequencies. From these beam traces, a Gaussian fit was performed to extract the beam-waist after propagation in the MSC structure.

Figures 5 (a) and (b) show the spatial distributions of the magnetic field energy for source beam-waists of $W_0=5 \,[a]$ and $W_0=40 \,[a]$ at the predicted MSC reduced frequency (dashed line in Fig. 5 (c) and (d)). In both figures, to facilitate the identification of the structures, the layout of the simulation cell is superimposed, together with the angle of propagation inside the structure predicted by the model $\Theta = \arctan [(d_1 \tan (\theta _1)+d_2 \tan (\theta _{\pi 2}))/D]$.

 

Fig. 5. (a) Magnetic field intensity distribution at the predicted frequency for $W_0=5\,[a]$; (b) Magnetic field intensity distribution at the predicted frequency for $W_0=40\,[a]$. In (a) and (b) the layout of the simulation cell is superimposed. (c) Waist evolution as a function of the reduced frequency for sources with different waist sizes. (d) Transmission spectrum for sources with different waist sizes. In (c) and (d) the vertical dashed line indicates the frequency corresponding to the predicted MSC. In (c) the color dashed lines represents the waist evolution for the propagation in bulk only.

Download Full Size | PDF

Concerning this angle of propagation, our model predicts an angular deviation of $\Theta \simeq 34.75$° with respect to the $x$-axis. For both Figs. 5 (a) and (b), the FDTD simulation gives $\Theta \simeq 35.86$° which is a fair agreement (relative error of $3.2\%$).

In the case of narrow beam size, we can see in Fig. 5 (a) the apparition of side beams below the main one. In first approximation, these mean that within each layer, the initial beams and those that have been reflected twice do not spatially overlap in a perfect way, which degrades the interference effect on which the AR property arises. This effect is more pronounced for smaller beams (which are wider in the $k$-space) and for higher angles. Conversely, Fig. 5 (b) shows that, for broader beams, no side beams are visible. The initial and twice reflected beams can spatially overlap, allowing a better interference and thus better AR.

From the beam trace, we can extract the beam waist size at the output of the MSC, and compare it to the waist of the source. From this, we can study the impact of the spatial extent of the incident beam on self-collimation. Figure 5 (c) shows the transmitted waist spectra. We can see that for $\Delta u \simeq 0.05$ around $u=0.208$ the transmitted beam waist is almost equal to the initial waist, which is the signature of self-collimation. This range is smaller than for standard SC in PhC, and is limited by the multiple internal reflections in the MPhC at oblique incidence. This point justifies a posteriori our strategy of privileging the AR efficiency rather than concentrating on self-collimation. One can also note that the best MSC frequency is well predicted by our simple model, as the minimum transmitted waist occurs at a frequency almost equal to the predicted one (vertical dashed line on the figure).

To further quantify the dependence between the AR quality and the beam waist, we computed the MPhC overall trasmittance $T$ for different source beam-waists. Figure 5 (d) shows that $T$ increases with the waist size, with values ranging from $70\%$ up to $95\%$ when the waist increases from $5 \,[a]$ to $40 \,[a]$. Indeed, a signature of the mesoscopic periodicity of the MPhC is clearly visible in the transmission spectra. In particular, the alternating slabs with period $D$ in the MPhC results in the opening of small band gaps in which the transmittance is reduced (see the dip in Fig. 5(d) at u=0.2091). This behaviour is entirely similar to that of a traditional multilayer structure. Besides, when the waist is large enough (i.e. $W_0>20 \,[a]$), a dip in the transmittance curve can be observed, appearing at the design frequency. Both these effects can be compared to the transmittance observed at normal incidence in [10], where values of $100\%$ were reached and the spectral signature of a small band-gap was observed slightly above the design frequency. We attribute the slightly lower values of transmittance and slightly blurred band-gap opening in our design to the large angle of incidence of the beam on the MPhC. Indeed, both these properties depend on the angular acceptance of the multilayer structure, which is known to decrease with the angle of incidence for simple quarter wavelength structures as used here. To further establish that the lower AR properties arise from the large angle of incidence, we tested our new algorithm at normal incidence using the same PhC structures as in [10], and found structures with very similar dimensions and performances (not shown here).

In order to further demonstrate the validity of our structure, we performed different 2D-FDTD calculations over a large simulation cell of $1200 \,[a]$ (83 mesoperiods). Figure 6 shows the steady state intensity of the magnetic field over the simulation cell for source frequency of $u=0.2091$ (predicted frequency for this MSC according to our algorithm) on Fig. 6 (a); $u=0.2083$ (best MSC frequency according to our previous simulations) on Fig. 6 (b) and a reference simulation in a bulk medium Fig. 6 (c). For all simulations the beam waist was kept equal to $20 \,[a]$.

 

Fig. 6. (a) Intensity map ($\lvert H_z \lvert ^2$) at frequency $u = 0.2091 \left [ a/\lambda \right ]$, which corresponds to the predicted frequency; (b) Intensity map ($\lvert H_z \lvert ^2$) at frequency $u = 0.0.2083 \left [ a/\lambda \right ]$, which corresponds to the best frequency for collimation through a 83 mesoperiods structure. (c) Intensity map ($\lvert H_z \lvert ^2$) at frequency $u = 0.2083 \left [ a/\lambda \right ]$ through bulk material.

Download Full Size | PDF

As seen, over such a length, the MSC is obvious when comparing the beam traces of the MPhC to the reference simulation. One can also see that at the predicted frequency (Fig. 6 (a)) narrow side-beams appear during the simulation. At the frequency offering the best MSC (Fig. 6 (b)) they have a much lower intensity. More precisely at the design frequency $u=0.2091$ the intensity transmission is $74.3\%$ and the transmitted beam has a waist of $19.9 \,[a]$. At $u=0.2083$, the transmission is higher, $82.5\%$, with a transmitted beam waist of $19.4 \,[a]$. For the reference run, the waist at the same position is already $39 \,[a]$.

From these results we can conclude that our methodology is sound and permits to study and design structures offering mesoscopic self-collimation over a non-trivial distance, under oblique incidence and under the light-cone.

At first glance, the overall transmission of the designed structure looks significantly lower than our previous results at normal incidence [10]. However, these previous 2D studies were not accounting for vertical losses and were limited to structure working way above the light cone, with extremely large vertical losses observed both on 3D simulations and experimental demonstration [12]. Here for the first time, the structures are designed to be simultaneously AR and below the light cone and vertical losses should remain small.

These results also highlight the main limitation of structures ensuring mesoscopic self-collimation at oblique incidence: as any multilayer structure under oblique incidence, they offer limited angular acceptance and when working with small beams, the AR degrades and side beams appear as clearly seen in Figs. 6 (a) and (b).

4.2 Impact of the residual curvature

The first studied structure presented a very low residual curvature ($FOM = -2 \times 10^{-4}$). We studied two other solutions, $S_2$ and $S_3$ in Fig. 4, with different FOM to investigate its impact on MSC performances. Both solutions are relatively similar but for the sign of the FOM and its magnitude (two orders of magnitude larger than that of $S_1$). Solution $S_2$ has a bulk length of $d_1=7.53 \,[a]$ and is designed for a beam at $u=0.2101$ incident at $27.46^{\circ }$, and exhibits a negative $FOM = -0.019$. Solution $S_3$ has a bulk length of $d_1=7.56 \,[a]$ and is designed for a beam at $u=0.2079$ incident at $26.78^{\circ }$, and exhibits a positive $FOM = 0.025$.

FDTD simulations for a 14-mesoperiod-long MSC structure were carried out for both $S_2$ and $S_3$ solutions (with similar resolution and settings than for $S_1$). Figure 7 shows the retrieved transmission (Fig. 7 (a)) and waist evolution (Fig. 7 (b)) for the three solutions. We choose to plot only the transmission spectrum, for the sake of figure simplicity. As seen the three structures exhibit characteristics that are fairly similar (if not exactly equal) to each other.

 

Fig. 7. For Designs 1, 2 and 3 and a propagation distance of about $200\,[a]$: (a) Transmission spectrum and (b) waist amplitude evolution, as a function of the reduced frequency. The dashed lines represent, for each design, the reduced frequency where the best self-collimation occurs. The blue line represents the waist after propagation in an homogeneous medium.

Download Full Size | PDF

Our interpretation of these comparisons is that one does not need a perfect MSC for a device of finite length. Provided that the residual curvature times the device length stays small enough as compared to the Rayleigh range of the incident beam, the output beam will be virtually unchanged.

4.3 Versatility of the algorithm

To further check the validity of our algorithm, we have also studied structures with varying sizes for the PhC slabs. We found many different structures achieving MSC under oblique incidence, and for illustration purpose, we show here one structure with $N=8$ rows of holes in the PhC slabs. Figure 8 shows the algorithm results for the $p=7$ AR branch for the PhC that supports angles of incidence similar to that of $S_1$. Once more, a large set of solutions can be predicted swiftly over a large range of frequencies. As expected from the larger PhC length, longer bulk sections (ie higher values of $d_1(u)$) are required to balance out the curvature accumulated in the PhC slab as compared to the $N=7$ structures studied previously. To limit the FDTD simulation cell whilst keeping a relatively high number of MSC periods, we choose to validate our calculations using structures that exhibit low $d_1$ values (see Fig. 8 (b), the blue stars identify the selected structures).

 

Fig. 8. Algorithm results (a) $\theta _1$, (b) $d_1$ and (c) $FOM$ for $N_{rows} =8$, we represent only the order $p=7$, the blue star shows the selected structure parameters.

Download Full Size | PDF

The 2D simulations are carried out using a MPhC about $200 \,[a]$ long, with still roughly 14 mesoperiods, like for $S_1$. The source beam waist is $20 \,[a]$. Figure 9 (a) shows the evolution of the magnetic field intensity at the reduced frequency predicted by the model (dashed line in Figs. 9 (b) and (c)). Figure 9 (c) shows the transmitted (green) and reflected (red) spectra and Fig. 9 (c) the transmitted waist as a function of the reduced frequency (solid green curve), together with the transmitted waist when propagating in a bulk material of same length (dashed green curve).

 

Fig. 9. (a) The magnetic field intensity at the predicted frequency. (b) Transmission and reflection spectra, and (c) the waist evolution of the transmitted beam, as a function of the reduced frequency. In (c) the green dashed curve shows the waist size after propagation in bulk material of same length. In (b) and (c) the vertical dashed line indicates the predicted frequency for MSC.

Download Full Size | PDF

From Fig. 9 (a), we see that the beam stays well collimated, with an propagation direction forming an angle of $\Theta = 30.3$° with respect to the $x$-direction, while the model predicted an angle of $\Theta = 29.3$°.

As seen in Figs. 9 (b) and (c), the structure exhibits good self-collimation performances: the transmitted waists is only $25 \,[a]$ wide with an energy transmission larger than $80\%$. Similar results were also found for $N=9$ structures. Increasing the size of the PhC increases the curvature to be compensated, and thus increases the size and $m$ order of the bulk section. The MSC quality is not affected by this change, but the angular and spectral acceptance of the AR condition is decreased, leading to an overall decrease of the transmission performances of the AR-MSC designs.

5. Conclusion

In this paper, we proposed a novel design method to achieve MSC under oblique incidence and under the light cone. The method allows to define the optimal geometrical parameters of MPhC through a simple hybrid numerical-analytical procedure. In particular, our algorithm targets solutions that simultaneously allow the propagation of self collimated beams along arbitrary directions (which are not of high symmetry for the crystal) and are under the light cone. This latter condition is particularly important for an efficient light confinement in the vertical plane and thus for reducing out-of-plane losses. Moreover, we prioritize the AR condition and then we find the solutions that best approximate the MSC condition. This choice allows to minimize the effect of multiple reflection between the bulk-PhC interfaces, which tend to degrade the self-collimation performance and can lead to difficult-to-predict beam spreading.

The self-collimation performance was quantified through a FOM defined as the ratio between the residual curvature per unit length in the MPhC structure and the curvature per unit length in the bulk material. The FOM indicates that the beam spreading of the designed MPhCs is less than $2~\%$ of the spreading expected in the bulk material.

The proposed design method was validated through a series of FDTD simulation examples showing a good agreement between the predicted overall angular direction $\Theta$ (defined with respect to the x-axis) and the FDTD simulated one (e.g. design beam direction equal to $\Theta \simeq 34.75^\circ$ and simulated beam direction equal to $\Theta \simeq 35.86^\circ$). Moreover, we evaluated the waist of the beam transmitted through the simulated MPhCs, showing that is is almost equal to the initial one in a wide spectral range around the nominal frequency $u$, thus confirming the occurrence of self-collimation.

We also put in evidence the dependence of the MPhC overall transmittance T, showing that higher transmittance (up to $95~\%$ in the proposed examples) is achieved for larger source beam-waists.

We demonstrated that parametric study of MPhC under arbitrary incidence is possible when reflectivity control becomes a priority. These parametric studies allow for fast MPhC design as we are able to have a global view of the geometry parameter space, thus opening the possibility to the conception of advanced MSC-based devices and interconnections. This procedure allows for parametric extensive design studies of MPhCs and provides meaningful starting points for FDTD validation analysis. As an extension, the numerical method developed in [13] could also be applied around the identified solutions to assess their angular acceptance. This would certainly provide a insight on both the angular acceptance of the structures and on their beam-filtering (as in k-space filtering) properties.