Mesoscopic self-collimation along arbitrary directions and below the light line

: Self-collimation (SC) and mesoscopic self-collimation (MSC) have been successfully demonstrated along the directions of high symmetry of photonic crystals. Indeed, wide angular acceptances are obtained only in these directions which oﬀer extremely ﬂat isofrequencies. In this article, we numerically demonstrate that mesoscopic self-collimation with large angular acceptance can be achieved along arbitrary directions that are not of high symmetry. In particular, we propose a simple method that allows to easily ﬁnd all the non-trivial collimation directions and corresponding frequencies. Thanks to the double periodicity of the mesoscopic crystal, these solutions can be eﬀectively tailored in terms of direction and frequency. Moreover, non-trivial MSC solutions can be found well below the light cone. These MSC features open up the possibility of designing complex systems by combining diﬀerent conﬁgurations, such as high reﬂection (HR) or anti reﬂection (AR) ones, or active materials.


Introduction
Self-collimation (SC) allows self-guiding propagation in photonic crystals (PhCs) without any typical guiding mechanism, such as index or band gap guiding. As explained for the first time in 1999 by Kosaka et al. [1], this phenomenon arises when, for a certain frequency, the corresponding isofrequency curve (IFC) of the PhC band diagram shows a flat region. In this region, for a whole range of wavevectors, the corresponding group velocity vectors are all parallel to each others as they stay perpendicular to the flat IFC. These wavevectors can thus form a beam that propagates inside the PhC without any lateral spreading. The angular spectrum of this beam is directly determined by the extent of the flat region of the IFC: a wider flat zone leads to a wider angular spectrum and thus to a smaller beam. In traditional PhCs, these flat regions usually occur at inflexion points of the IFC [1] and are trivially angularly narrow and only allow the propagation of broad beams. However, self-collimation can also occur along directions of high symmetry of the PhC where angularly wider flat regions can be achieved [2] and smaller beams can propagate under self-collimation. All of the proposed [3][4][5][6][7][8][9][10][11] or demonstrated [2,[12][13][14][15][16][17][18] devices rely on self-collimation along a direction of high symmetry of the supporting PhC. Moreover, broad-band and broad-angle operation have been achieved in bi-dimensional PhCs [19][20][21][22].
By combining alternating slabs of PhC and bulk materials, one can form a mesoscopic photonic crystal (MPhC) [23]. A MPhC can be seen as a 1D supercrystal constituted by slabs of a 2D PhC and bulk material. This kind of structures support mesoscopic self-collimation (MSC), introduced in [23], which consists in a succession of small defocusing, in the bulk slabs, and refocusing, in the PhC slabs, that compensate each other, resulting in a perfectly collimated beam. This occurs when the PhC used in the slabs is designed to work in its focusing regime, beyond the SC, showing an abnormal lateral dispersion. As previously demonstrated along the ΓX-direction of the supercrystal lattice, MSC only requires that the total spreading over one mesoscopic period averages to zero. Therefore, since it depends neither on the overall mesoscopic periodicity nor on the phase index in the PhC, one can easily design a mesoscopic structure that also exhibits many interesting optical properties as tailored overall reflectivity, slow-light and stable full optical confinement by 1D Fabry-Pérot-like microcavities [23][24][25][26][27][28].
However, so far MSC has been restricted to the directions of high symmetry for the MPhC, i.e. the direction normal to the slabs interfaces that form the 1D mesostructure.
In this contribution, we introduce a simple and generic method to find all the non-trivial self-collimation solutions in a mesoscopic photonic crystal. These solutions exhibit large flat IFCs that can be determined using a simple analysis based on plane wave expansion method (PWEM). Based on this method, we demonstrate that MSC can be achieved along arbitrary directions, and not necessarily along direction of high symmetry. We also demonstrate that, for certain conditions, MSC can be achieved below the light cone, paving the way towards low-loss MSC propagation.

Structure and methods
As depicted in Fig. 1, we will consider a MPhC consisting of a periodic arrangement of slabs of bulk material (length d b , refractive index n b = 2.9, i.e. effective index of a GaAs membranes at 1.55 µm wavelength [28]) interleaved with slabs of a 45 • -tilted square-lattice PhC (length d c , lattice constant a, radius of the air holes r = 0.28 a, etched into the same material that constitutes the bulk slabs). The x-direction is assumed normal to the interfaces of the slabs. The length of each row of holes is a/ √ 2 and, accordingly, we have d c = aN row / √ 2, where N row is the number of rows constituting the PhC slab. The overall periodicity is D = d c + d b . Inside the PhC slabs, the x-direction corresponds to the ΓM c direction (see Fig. 2 In this manuscript the length is taken in unit of a and reciprocal lattice vector in the k-space is assumed to be in normalised unit of 2π/a. Along the x-direction, the periodicity of the MPhC arises from two contributions: the rotated PhC, with a 2D periodicity, and the bulk medium, aperiodic. Along y-direction, its periodicity is fixed only by the rotated PhC and, therefore, is a √ 2 periodic. Consistently, the axis of the first Brillouin zone (FBZ) of the MPhC extends as follows: the k x -axis from −a/2D to a/2D, whilst the k y -axis from − √ 2/4 to √ 2/4. It is worth pointing out that, in the k x -direction, the length of the ΓX 1 segment of the MPhC is √ 2D/a times smaller than the ΓM c segment of the PhC. Furthermore, since the proportion between the FBZs of the two structures is fixed by N r ow , d b and a, in order to have an integer number of bands N bands to cover the entire ΓM c segment, the d b has to be chosen according to the following equation: In order to illustrate our approach, we consider a first structure with N row = 3 and d b = a/ √ 2. For such a structure we have exactly ΓX 1 = 1 4 × ΓM c , as seen in Fig. 2(a) that depicts the FBZ of the MPhC superposed to that of the PhC (tilted by 45 • ). This means that exactly 4 bands of the MPhC are necessary to cover the same k-range than the first band of the PhC. We thus call this structure F 4 as in "4 Foldings". These 4 foldings can be seen comparing Fig. 2(b), which presents the first band for the PhC, and Fig. 2(c), which presents the first 4 bands of the MPhC, properly unfolded.
The band structure involved in this analysis was calculated using a well established PWEM algorithm [29]. The 2D unit cell shown in Fig. 1, surrounded by periodic boundary conditions, has been discretised with a resolution of 32 points per a. The band structure is obtained by calculating the reduced eigenfrequencies f k x , k y , in unit of a/λ, across the whole FBZ for the four first bands, using a resolution for the ì k-plane of 128 point per 2π/a. The resolution was chosen in the simulations after a convergence analysis and it gives a good compromise between computation time and result stability.
The IFCs-diagram is then retrieved by projecting each band on the ì k-plane. However, as the unit cell of the MPhC is four times larger than the plain PhC one, its FBZ is four times smaller, leading to 4 foldings. To ease the understanding of the tangled birefringence of the MPhC, an unfolded representation of its IFCs-diagram is provided ( Fig. 2(b)). To further provide a direct comparison with the IFCs-diagram of the 45-degree-tilted plain PhC, a portion of the same size of the corresponding MPhC FBZ has been depicted in Fig. 2(c). Comparing Figs 2(b) and 2(c), one can see that the overall structure of the IFCs of the PhC are still present in the IFCs of the MPhC. The meso-periodicity mainly introduces new IFCs distortions, most notably around each of the 4 foldings. These distortions will provide new locations for non-trivial MSC solutions, as will be discussed later on. As it has been shown in [23], the mesoscopic self-collimation, along the high symmetry ΓX 1 -direction, requires a null curvature locus on a given IFC. The curvature κ of a certain IFC, at a give point k x , k y , can be estimated from the unfolded band structure, according to [30], as follows: where the x and y subscripts indicate partial derivatives with respect to k x and k y and where ì v g is the group velocity. However, only asking for a null curvature does not suffice to ensure wide MSC in directions that are not of high symmetry. Indeed, every inflexion point, where the flatness of the IFC is trivially narrow, will show a null curvature. In order to find wide MSC zones we have to ask for IFC to be maximally flat in the range of interest.
To determine the conditions ensuring maximally flat IFC, let us consider a Gaussian beam S I (u, v) at a frequency f 0 , that propagates energy along an arbitrary direction u, being v the local transverse direction. By definition of the IFC, the local coordinate system is such that u is normal to the IFC while v is tangential to the IFC. In the [k u , k v ]-space it can be described -plane, the f 0 -IFC, given by the dispersion relation, can be locally expanded as a Taylor series: where a n = 1 and a 0 = 0. Thus, the generic propagator in the u-direction, for an arbitrary medium, can be expressed as it follows: The source, after being propagated for a distance L u , can be written as it follows: Therefore, in order to have a maximally flat IFC and to obtain a collimated beam, all the coefficients a n have to be zero when the u-direction corresponds to the propagation direction. This is a stricter set of constraints than what is usually asked for when searching for SC or MSC solutions [23,31]. Hereinafter, we will only consider the first three coefficients a 1 , a 2 , a 3 to retrieve the position P 0 in the k−space of the maximally flat MSC ranges. These can be written as: To calculate numerically these coefficients, we rely on the very definition of the IFC under study as f (k u , k v ) = f 0 , it implies that its total differential is null : This allows to express a 1 , a 2 and a 3 as a function of the partial derivatives of the band (e.g. we can write a 1 = − f v / f u and so on).
where the u and the v subscripts indicate partial derivatives with respect to k u and k v . However, these coefficients are still expressed in the local coordinates (u, v). In order to express them in the (x, y) coordinates and to calculate them numerically, we have to perform a rotation of the reference axis by θ = arctan(v g,y /v g, x ), where v g, x and v g, x are the components of the group velocity. As the group velocity is always normal to the IFC, after performing this rotation we find that the condition a 1 = 0 is always satisfied. Therefore, in order to find the position of the maximally flat MSC ranges we only need to look for a 2 = a 3 = 0. Graphically, requiring the IFC to be maximally flat in P 0 implies that the IFC and the curve κ(k x , k y ) = 0 must have the same tangent at this point: this provides a geometrical insight for the MSC condition. Moreover, this also provides an alternative way to search for it. Since the tangent angle of a given IFC and of the curve κ(k x , k y ) = 0 can be written as φ I FC = arctan(− f x / f y ) and φ κ=0 = arctan(−κ x /κ y ) κ=0 , respectively, we can also ask for |φ κ=0 − φ I FC | to be zero along the zero-κ locus.
Finally, it is worth observing that the a 2 can be related to the curvature κ (which is invariant with respect to axis rotations) as it follows: In order to validate this method we apply it to the F 4 structure introduced previously (N row = 3 and d b = a/ √ 2, having N bands = 4). We find the MSC points by numerically computing the respective unfolded band structure and by looking for the crossing points of the zero a 2 and a 3 loci. Then, some of these working points are simulated using 2D-FDTD [32] with a computational cell surrounded by scalar absorber boundary condition (electric and magnetic conductivity losses that slowly increase within the absorbing layer). The simulation region is excited by a continuous wave (CW) source. Two different profiles are used: either an omnidirectional point source or a tilted Gaussian spatial profile, with a full width at half maximum of the intensity ∆ x = 10 × a. Moreover, by computing the Fourier transform of the z-component of the magnetic field (h z ), it is possible to directly observe the IFCs responsible for the propagation of the self-collimated beam.

Results
Figure 3(a) shows the loci satisfying the conditions a 2 = 0 (black lines) and a 3 = 0 (gray lines), superimposed on the IFC of the F 4 configuration. As it can be seen, several points satisfy the MSC conditions and provide non-trivial flat zones. The length of these flat zones strongly depends on the particular position of the MSC point. Wider flat zones will allow self-collimation of narrower beams. It is worth observing that, owing to the higher order derivatives involved, the numerical computation of the a 3 = 0 curve is noisier than the a 2 = 0 one. Numerical noise is the reason why, in practice, it is not meaningful to try canceling higher order terms of the expansion given in equation (3).
In particular, the thick segments in Fig. 3(b) highlight the flat zone for each of them. Moreover, the same figure displays both the phase (black thin arrows) and the group velocity (grey thick arrows) vectors for three specific solutions labelled P1, P2 and P3. It is worth recalling that the group velocity vector determines the direction of propagation of the self-collimated beam.  (Fig. 4(a)), at f = 0.1609[a/λ] (Fig. 4(b)) and at f = 0.2043[a/λ] (Fig.4(c)) that correspond to the reduced frequency of the IFCs onto which the points P1, P2 and P3 (see Fig. 3(b)) are located. Figs. 4(d)-4(f) show the corresponding FDTD calculated IFCs in the FBZ, obtained by the spatial Fourier transform of the propagated magnetic fields. As it can be inspected by comparing Fig. 4(d)-4(f) and Fig. 3(b), the shape of the FDTD calculated IFCs is in complete agreement with that obtained by means of PWEM calculations. To ease the figure reading, the beam formation has been The omnidirectional point source is capable of exciting all the self-collimated beams at once as shown in Fig. 4. However, we can excite a specific non-trivial solution independently by imposing a directive source. This can be numerically achieved by a Gaussian beam, provided that the phase components of the excitation lay on a IFC flat region. A case of study that gives an example of this property is a MPhC forming a flat vertical interface (parallel to the y-axis) with the hemispace (x > 0) filled by the background bulk material. We performed the simulation of this structure excited by a tilted source having a Gaussian spatial distribution (∆ = 10 × a) placed near the lower left corner of the computational cell. The source axis (the direction normal to  its waist corresponding to the Gaussian beam axial direction) forms an angle Θ s with respect to the x-axis. Figure 5 shows the Poynting vector magnitude calculated by FDTD for the F 4 MPhC configuration forming a flat interface with the bulk medium. The structure is excited by a Gaussian beam with f = 0.1368[a/λ] and Θ s = 70 • , corresponding to the point P1 in Fig. 3 (b).
It is worth pointing out that the solution P1 is below the light cone and, therefore, the excited beam propagates with no intrinsic out-of-plane scattering losses. By inspecting Fig. 5, we can see that, once the beam has reached the flat interface, it is reflected and refracted according to the generalised Snell Law [33]. According to the power conservation law, the reflected and refracted beams exhibit different beam widths and their intensities depend on the incidence angle on the interface. Moreover, the reflected beam remains collimated, as it still propagates inside the MPhC, whereas the refracted beam spreads as expected for the propagation in the bulk medium. In order to better show the versatility of the proposed design method, we consider another MPhC configuration, already studied by the authors [24]. This configuration, named AR1, was designed to offer low reflectivity at the interfaces between the bulk material and the MPhC for a beam propagating along the direction of high symmetry. Using our method, we found alternative solutions along arbitrary directions. In particular, here we analyze the excitation of a tilted self-collimated beam in a MPhC layer delimited by two semi-infinite layers of bulk medium (the interfaces are at x < −90a and x > 100a). The aim of this example is to show that we can excite a self-collimated beam at a well defined SC direction, from outside the MPhC, by a Gaussian source placed in the bulk medium. Figure 6(a) shows the calculated Poynting vector distribution for the AR1 MPhC when it is excited by a Gaussian beam in the bulk medium with f = 0.1868[a/λ] and incidence angle in the bulk Θ s = 18.7 • . Considering the refraction phenomenon at the bulk-MPhC interface, the input beam is capable of exciting a self-collimated beam in the MPhC propagating at angle Θ s = 13.2 • . In conclusion, in order to excite the SC beam at a given angle (i.e. Θ s = 13.2 • ) from the bulk medium, the source incidence angle must be properly chosen to satisfy the conservation of the wavevector component parallel to the interface. To identify the required incidence angle, we can inspect Fig. 6(b) that shows the logarithm of the self-normalised absolute value of the spatially Fourier transformed z-component of the magnetic field. The superimposed dashed black curve represents the IFC pertaining to the AR1 MPhC calculated by means of PWEM, whereas the dot-dashed thin black curve are the corresponding bulk isofrequency circle. The blue and the violet arrows denote the wavevectors of the MPhC and of the bulk, respectively, whereas the pale and the dark green arrows are the MPhC and the bulk group velocity vectors, respectively. The wavevector component k // parallel to the interface is highlited by the red arrow. The conservation of the wavevector parallel component occurs when the beam impinges on the interface with an angle Θ i n = 18.7 • . This is in very good agreement with the full-wave FDTD results of Fig. 6(a). Similar considerations can be made for the output interface at x = 100a.

Conclusion
In this paper, we have introduced a simple numerical method for determining all non-trivial solutions (with wide angular acceptance) of mesoscopic self-collimation. Using this method we have numerically demonstrated MSC along non-trivial directions, while MSC was so far limited to highly symmetrical directions for the structure. Albeit not shown in this contribution, this method is not restricted to MSC and works equally well for normal self-collimation. It is worth pointing out that by changing the composition of MPhC, it is actually possible to adapt the sets of angles and frequencies in which MSC can be achieved. In addition, we have found MSC solutions even below the light cone, while MSC had so far been proposed and observed only above the light cone, resulting in large propagation losses. This paves the way for low-loss propagation of MSC. Finally, MSC solutions can be combined with other optical properties, such as minimal reflectivity at the PhC interfaces for efficient coupling with the surrounding environment.