Directional spontaneous emission in photonic crystal slabs

: Spontaneous emission is one of the most fundamental out-of-equilibrium processes in which an excited quantum emitter relaxes to the ground state due to quantum fluctuations. In this process, a photon is emitted that can interact with other nearby emitters and establish quantum correlations between them, e.g., via super and subradiance effects. One way to modify these photon-mediated interactions is to alter the dipole radiation patterns of the emitter, e.g., by placing photonic crystals near them. One recent example is the generation of strong directional emission patterns—key to enhancing super and subradiance ef-fects—in two dimensions by employing photonic crystals with band structures characterized by linear iso-frequency contours and saddle-points. However, these studies have predominantly used oversimplified toy-models, overlooking the electromagnetic field’s intricacies in actual materials, including aspects like geometrical dependencies, emitter positions, and polarization. Our study delves into the interaction between these directional emission patterns and the aforementioned variables, revealing the untapped potential to fine-tune


Introduction
Spontaneous emission is a paradigmatic out-ofequilibrium phenomenon that depends on both the emitter and its interaction with the surrounding media, as proposed by Purcell [1].It consists of a quantum system transitioning from an excited energy state to a lower energy state by emitting a photon with quantized energy equal to the difference between the two energy states.This process depends on the intrinsic properties of the emitter, such as the emitter's transition frequency and dipole moment, as well as the surrounding electromagnetic environment, which can be modified by placing nearby materials.This opens the exciting possibility of tuning such fundamental phenomena by positioning the emitter close to macroscopic materials, like cavities [2], placing nearby emitters [3,4], or combinations of both.
One of the most paradigmatic examples of modifying individual and collective spontaneous emission of quantum emitters are photonic crystals (PhC) [5][6][7], which are now routinely interfaced with many different types of emitters [8][9][10][11][12][13][14][15][16][17].In these systems, the propagation of light becomes modulated by the periodic variation of the refractive index, resulting in nontrivial band structures, e.g., featuring photonic bandgaps [18], Dirac-like energy dispersions [19] or saddlepoints [20][21][22].These can substantially modify the individual and collective radiative patterns of emitters placed nearby.For example, when the emitter's optical transition lies within photonic band-gaps, spontaneous emission is largely suppressed, leading to the formation of atom-photon bound states [23][24][25][26].These states can induce perfect coherent exchanges between emitters with different shapes depending on the lattice geometry [27][28][29][30][31][32][33].In the opposite regime, when the emitter lies within the band, spontaneous emission occurs but with renormalized lifetimes and radiation patterns.A particularly intriguing situation is the spontaneous emission modification that occurs in twodimensional photonic crystals due to the presence of saddle-points and straight iso-frequencies in the band structures [20,21].While it has long been known that such energy dispersions can lead to unconventional light propagation, such as supercollimation [20,21], it has been recently pointed out that when emitters are resonant to these frequency regions, spontaneous emission can be accelerated to the point where it becomes non-perturbative [34], and their individual radi-ation patterns become highly directional [34][35][36][37][38] leading to strong super and subradiant effects when several emitters are present [34,35].However, most of these predictions were made for simplified coupled-mode descriptions [34][35][36][37], in which the complexity of the electromagnetic field, such as its geometrical dependence or polarization, were neglected, or for particular emitter positions [21,22,38].The interplay of these factors with such unconventional individual and collective phenomena remains an open question.
In this work, we undertake the challenge of studying the interplay of such directional emission with the emitter's dipole polarization and the emitter's position in a particular two-dimensional PhC slab, i.e., the honeycomb structure studied in Refs.[30,31].To achieve this, we theoretically investigate the electromagnetic wave propagation of dipole emitters with different positions and polarizations using the Guided Mode Expansion (GME) [39] and the finite difference time domain (FDTD) [40] algorithms.From this study, we observe the expected directional emission into six rays predicted by toy-models [41] but also note that some of these directions can be enhanced or suppressed depending on the emitter's polarization and position.The manuscript is organized as follows: In Sec. 2, we introduce the PhC structure we are considering and the theoretical methods and foundations used in this work.Sec. 3 discusses the relationship between the classical Green's function and the individual and collective spontaneous emission phenomena.In Sec. 4, we describe the emission of a dipole emitter at the center of the unit cell with circular polarization.Sec. 5 examines different emitter positions and dipole moment polarizations to analyze their interaction with directional emission.Finally, we summarize our findings and outline potential applications in Sec. 6.

Photonic structure
We consider the photonic-crystal (PhC) slab depicted in Fig. 1, introduced for the first time in Ref. [30] as a way of frequency-isolating Dirac points.The slab consists of a hexagonal lattice of GaP ( GaP = 10.5625) with a thickness  and primitive vectors a 1/2 = ( √ 3/2, ±1/2), where  is the lattice constant.The slab has six bigger and six smaller air holes of radii   and   , respectively, as shown in Fig. 1(b).The bigger holes have a distance to the unit cell center, denoted by ..(e)and (f) Electric field amplitude for the first band mode for  1 (e) and  2 (f), the vectors  1 and  2 in the irreducible Brillouin Zone are showed as green arrows in (d).
To compute the eigenfrequencies and eigenmodes of this PhC structure, we use the GME proposed in Ref. [31] for this structure.The GME method uses the solutions of an effective homogeneous slab denoted by ℎ g, (r, ) to expand the modes of the PhC slab.Here, we use r to denote the position in the plane of periodicity ( plane) and put explicitly the  coordinate.In principle, the homogeneous slab has two types of solutions: radiative modes oscillating outside the slabs and guided modes that decay exponentially outside the slabs.Initially, the method considers only the guided modes that are characterized by ℎ g, (r, ) =  g•r h g, (r, ) [31,39], where the vector g = k + G is the wave vector in the  plane, with G a reciprocal lattice vector that contains information of the system periodicity and k a wave vector in the first Brillouin zone; and  is an index for the quantization of the guided modes.By using these modes, the magnetic field for the PhC slabs is expanded in the frequency domain (harmonic fields) as: where   (k + G, ) are the coefficients of expansion and  denotes the different modes of the PhC slab that appear with frequency  k, .A matrix eigenvalue problem is obtained by putting these fields in the timeindependent wave equation for the magnetic field where the sum over g = k + G is equivalent to the sum over the lattice vector G because k is fixed for each eigenvalue calculation.The matrix elements ℋ g,;g ′ , read where   is the integration volume defined by the unit cell in the -plane and (−∞, ∞) in the -direction; details about the analytic expression for matrix elements can be found in Ref. [39].In principle, equation (2) considers infinite guided modes in the expansion of fields, so the eigenproblem is infinite.However, the guided modes basis must be truncated to a finite size to compute it numerically.Solving the finite eigenproblem, we obtain the eigenvalues 2 and the eigenvectors containing the coefficients   (g, ) that expand the magnetic field mode H k, (r, ) in the truncated basis.The band diagram of our photonic structure is displayed in Fig. 1(c) for the TE-like (even) modes.We are interested in describing the emission of the dipole source in the middle of the slab ( = 0) with components in the  plane, so the TM-like (odd) modes are neglected because they have only -component at  = 0. Here, the first band has a frequency region near the symmetry point  = 2 3 ( √ 3, 0) that presents a saddle-point, which will result in an associated Van-Hove singularity (VHS).To analyze in detail the directional emission in this frequency region, it is necessary to consider the contribution of all the k vectors in the first Brillouin, which is why in Fig. 1(d), we plot a density plot of the first band.The iso-frequency curve at   is displayed by the black-dashed line showing a quasi-straight hexagonal curve, which suggests that the contribution of this k vector propagate in six direction perpendicular to each side of the hexagon.Additionally, for all the k vectors in this curve, we calculate the group velocity through the following relation [42]: where S corresponds to the Poynting vector, and  E(H) denotes electric(magnetic) energy density integrated into the volume   .To compute the Poynting vector, we use the magnetic field previously calculated using the GME method.The group velocity for the k vectors in the iso-frequency curve point almost preferentially in six directions, each one perpendicular to a side of the iso-frequency, and the maximum contribution of the group velocity is at specific k vectors in these specific directions, that are angles of /3 + /6 rad, with  = 0, 1, . . ., 5.These results explain the self-collimation phenomena that cause the directional emission [21,22].Finally, in Fig. 1(e,f), we present the electric field distribution for two k vectors in the isofrequency curve that are represented with green arrows in the panel (d) of the same figure.The field distribution shows that for every single k, the maximum field in the unit cell changes.This ultimately means that the coupling of the emitter with the PhC structure will strongly depend on which position in the unit cell the emitter has.This is what we will explore in the next sections.

Classical Green's function and quantum spontaneous emission
To calculate the effect of the emitter positions and the field polarization in the directional spontaneous emission near PhC, we use the macroscopic QED formalism [43].This formalism gives an effective description for the quantum emitter dynamics after tracing the photonic degrees of freedom of the PhC.Assuming we have  identical emitters with transition frequency   , the description is given in terms of a Born-Markov master equation which reads: where  is the density matrix describing the emitter degrees of freedom,  0 = ∑︀  ℏ     corresponds to the independent emitters' Hamiltonian, and    = |⟩  ⟨| is the notation introduced for the emitter operators.Besides, as a result of the interaction with the photons, the emitters' dynamics contain two additional terms,  eff that represents the coherent part of the photonmediated interactions, which read: and, ℒ eff () that describes the incoherent dynamics induced by the photonic bath that reads: The key point of this formalism is that both the coherent and incoherent photon-mediated interactions,  , and Γ , , respectively, are expressed in terms of the classical Green's function as follows: with d  being the dipole moment of the -th atom.
Here, Γ , is known as the cooperative or collective decay rate for  ̸ = , and corresponds to the spontaneous emission rate of a single emitter for  = .When several emitters are present, both the individual and collective decay terms strongly renormalize their spontaneous emission.To evidence it, we plot schematically the effect of Γ , and  , in the Hilbert space of two twolevel emitters composed by the states: 2. There, we observe how   shifts the energies of the |±⟩ states, whereas Γ  renormalize their lifetime to be Γ  ± Γ  .
Summing up, calculating the real and imaginary parts of classical Green's function, we indirectly have access to the collective dynamics of emitters.This is why, in the next sections, we focus precisely on calculating the classical Green's tensor defined by the equation: where, G(r, r ′ , ) denotes the Green's tensor with components  , (r, r ′ , ).Its classical meaning is that a point source oscillating with frequency  and dipole moment d at the position r ′ produces an electric field at the position r determined by the following expression: To find Green's function, we use an FDTD algorithm that numerically calculates the field in the time domain produced by a dipole source.Finally, to obtain Green's tensor in the frequency domain, we perform a fast Fourier transform considering a time window long enough not to include the effects of turning the source on/off.By simulating the field produced by a dipole with only -component, we can calculate the Green tensor components   and   and do the same for a dipole with -component to find   and   .In what follows, we are interested in putting the source in the middle of the slab where the -component of the TE-like modes is zero, so   =   = 0.

Directional spontaneous emission
As mentioned in the introduction, recent works in the literature with toy-model honeycomb geometries [41] predict a highly directional radiation pattern when the emitter is tuned to frequency regions with straight isofrequencies and saddle-points.In this section, we intend to explore what happens in the PhC slab structure described in Section 2. For that, we use Ansys Lumerical -FDTD software to numerically simulate a PhC slab in a finite hexagonal region with 19 unit cells from side to side.After that, we consider a dielectric homogeneous slab with permittivity   .For concreteness, we consider a lattice parameter  = 0.28044 × 10 −6 m, corresponding with a frequency   =   /2 = 399, 595 THz.The rest of the geo-metrical parameters are fixed in terms of  equal to the ones described in Fig. 1.
We start considering a dipole source at the center of the unit cell with a frequency   and a left-handed circular polarized ( + ) dipole moment, i.e., d = dσ + , where σ± = ∓(x ± ŷ)/ √ 2. In Fig. 3(a), we plot the corresponding electric field intensity radiated into the structure, showing that in this case, the emitter radiates into the structure along in six directions at ±/6 rad, ±/2 rad, and ±5/6 rad.Thus, for this emitter's position and polarization, the behavior is qualitatively similar to the one predicted by toy-models [41].
To make a more detailed analysis of the radiation, in Fig. 3(b), we make a polar plot of the electrical field intensity at a fixed distance  = 4 from the source in terms of the angular coordinate .Like this, it is more clear that field intensity, represented by the radial coordinate, is maximal slightly shifted from the /3 + /6 rad.The underlying reason is that the geometrical structure of the unit cell favors the field to accumulate near the air holes.Apart from that intensity information, we also include a color map in the polar plot, which represents the  1 Stokes parameter defined as: where,  () represents the ()-component of the electric field.With this definition,  1 = 1 (red color) indicates horizontal (H) polarization, while  1 = −1 (blue color) indicates the vertical (V) one [44,45].With this color code, one can immediately notice a strong relation between the direction of propagation and the polarization, reminiscent of what happens in one-dimensional chiral quantum optical systems [46], where certain propagation directions are locked to certain polarizations.For example, along the vertical propagation around ±/2 rad, the field is clearly Hpolarized, while in oblique directions, the V polarization dominates.These results suggest that selecting the emission direction using different linear polarization should be possible, as we explore in the next section.
Finally, to ensure that the peaks in the polar plot are not an artifact of the distance chosen (4), we have analyzed the field intensity as a function of the distance for the four directions depicted with dashed lines in Fig. 3(a).The field intensity in these directions is plotted in Fig. 3(c), showing that directions at /6 rad and /2 rad dominate over the other one at /12 and /3 in almost all distances.

Interplay of directional emission with field polarization and emitter position
From the results obtained in the previous section, it became clear that there is a distinct interplay between the direction of propagation of the field and its polarization.In this section, we aim to explore how one can use this degree of freedom, or others like the emitter's position, to tune the directional emission patterns.In particular, we will first study the dependence of linearly polarized dipoles at the center of the unit cell, and then consider the effect of the emitters' positions by moving them to locations outside of the unit cell center.

Linear polarization
Let us here consider a linear dipole source at the center of the unit cell.As discussed above, the paths in the vertical direction ( = /2 rad) have H polarization, so we start with a linear dipole in this polarization at the center of the unit cell, as shown in the inset of Fig. 4. In this case, the dipole emitted light preferentially in the vertical direction, as the intensity profile shown in Fig 4(a).This figure shows the maximum field intensity at angles around ±/2 and no fields at the other four main paths (±/6 and ±5/6 rad).To observe more clearly the directional emission, we appeal to the polar plot previously explained in Sec. 4. In the case of H polarization, the polar plot in Fig. 4(b) shows only peaks at ±/2 rad, which means that only one of the three main directions was selected by this linear polarization, while the other directions are turned off.Additionally, we analyze other dipole moment orientations, shown in panels (c-d) of Fig. 4. In panel (c), we consider a dipole moment d = √ 3/2x + 1/2ŷ, which excites the paths at −/6, ±/2 and 5/6 rad.In this oblique case, two directions of propagation are activated because the perpendicular direction of the dipole orientation corresponds to 4/6 rad, which is in the middle of the two excited directions.This is because, in this case, the dipole couples with the same amplitude to the modes of the PhC slab that propa-gate in these two main directions and does not couple to the modes propagating in the third direction (/6 rad).Similar to this case, one can consider any arbitrary orientation of the dipole to excite two of the three directions, except when considering orientations of 0, ±/3, and ±2/3 rad.Finally, the particular cases of dipole orientations 0, ±/3, and ±2/3 rad select only one of the three main directions corresponding with the perpendicular direction of the dipole orientations, as shown in Fig. 4(d) for a dipole orientation of /3 rad.In the Supplementary Material, we provide an animated video displaying the field intensity distribution produced by different dipole orientations (see Supplementary Material, Video 1), showing continuously how the radiation patterns are modified as one moves the polarization of the dipole continuously.
Let us now discuss the interplay between the field polarization and the three directions in more detail.In the vertical directions, we have H polarization equal to the dipole orientation that excites this direction.However, in the oblique directions, the polarization becomes a linear combination of both H and V polarization equal to the combination of dipole orientation that activated these oblique directions at /6 and 5/6 rad, where the dipole orientation and field polarization are defined by d = ± √ 3/2x ± 1/2ŷ.This possibility of choosing the direction by using linear polarized dipoles is allowed by the fact that the electric field at the unit cell center always has a non-zero value for all modes in the iso-frequency, and thus, the dipole source can always couple to all momenta irrespective of the polarization.As we will see in the next section, this is not true for all positions of the unit cell.

Linear polarization and emitter out of center
When the emitter position is the unit cell center, the PhC modes at all momenta contribute to the propagation, including all possible linear polarization contributions.However, there are other positions outside the center where some modes have a non-zero electric field while others do not have an electric field, see panels (e-f) of Fig. 1.There, we observe modes for two wave vectors in the iso-frequency curve, one pointing in the horizontal direction (panel e) and the other in the vertical (panel f).The latter has a null electric field in the dielectric region between the two nearest bigger air holes, for example, around r = (0, 0.35).In contrast, the former has a significant field contribution in this region.This difference leads to substantial differences in the control of the directional emission.To evidence that, let us consider a linear polarized dipole at position r = (0, 0.35).In this situation, the intensity profile of the electric field displays a vertical wave propagation as shown for an H-polarized dipole in Fig. 5(a), similar to what happens at the center of the unit cell.However, if the dipole orientation changes to an oblique angle, i.e., d = (x + ŷ)/ √ 2, the intensity profile is still vertical, as shown in Fig. 5(b).Let us, however, note the maximum intensities in the color scale, where we can observe that the intensity decreases by approximately a factor of two from the H-polarized dipole to the oblique one.Similar results for other orientations and circularly polarized dipoles are found, as shown in the videos we provide in the Supplementary Material (see Videos 2-3).Additionally, we must mention that the vertical propagation has a well-defined H polarization independent of the dipole moment.
To understand the origin of this directionality that is independent of the polarization, we analyze the contribution of each mode from the iso-frequency curve to the field propagation by use of the field intensity [Fig.5(c)] and Stokes parameter [Fig.5(d)] of each mode at this position r = (0, 0.35).First, we discuss it in terms of the intensity that gives us information about the coupling of the dipole with each mode.We find that the modes with more intensity correspond to the wave vectors near the vertical directions, which explay why they dominate the radiative decay of an emitter in this position.This explains the directional emission, but it does not explain neither its polarization nature nor the intensity decrease from panel (a) to (b) of Fig. 5.To understand the invariant polarization of the emitted field, we examine the  1 parameter at r = (0, 0.35) for all the modes in the iso-frequency, as shown in Fig. 5(d).Here, we observe that almost all the modes in the iso-frequency curve have only H polarization, which means that only the horizontal component of the dipole couples to the PhC structure and couples with half intensity due to the smaller overlap with the predominant polarization of the modes.

Circular polarization and emitter out of center
Finally, let us analyze the case where the emitter is in the center of the smallest air hole at r = (0, / √ 3).In contrast to the center of the unit cell, here the six principal directions are not equal; for example, if we consider the vertical direction, the photonic structure is different upward than downward from the center of air holes.We expect that this anisotropy causes an upward and downward difference in the emission patterns that can be controlled by the polarization of light, for example, considering circular polarizations  ± .To test this intuition, we consider first a left-handed circular polarized ( + ) dipole source, which selects three of the six emission paths, as shown in Fig. 6 where emitted electric field intensity is displayed.In this case, the  + polarization can be used to filter the paths at /2, 7/6, and 11/6 rad, turning off the other main paths at /6, 5/6 and 3/2.The choice of paths can be inverted if we consider right-handed circularly polarized ( − ) dipoles, where the paths activated are at /6, 5/6 and 3/2 rad, while the inhibited paths are at /2, 7/6 and 11/6 rad, as show the intensity profile in Fig. 6(b).The change in the electric field pattern passing from left-to right-handed circularly polarized dipole is shown in the Supplementary Material (see Video 4).
Additionally, for these two cases,  + and  − dipoles, we quantify the circular polarization of the radiative modes by using the Stokes parameter  3 [44, 45], defined as: where  () are the Cartesian's components of the field and  +(−) are the components in the circular polarization basis (σ ± ).The  3 = ±1 indicates left and righthanded circular polarization, while intermediate values indicate both contributions that can be elliptical or linear polarization.The activated paths contain several centers of the smallest air holes connected by a lattice vector, specifically by a lattice parameter in the radial direction.For the case of the  + emitter, each center in the three activated paths has a dominated  + polarization.In contrast, the field has an almost equal contribution of both left and right polarization for other distances throughout the path, as we can observe from the  3 parameter in Fig. 6(c).For the case of  − , the  3 parameters are plotted in Fig. 6(d) as a function of the distance in some of the three activated paths (the rest show qualitatively similar physics), where one can observe that at each center the  − dominate over the  + .In contrast, both contributions are similar for other distances, suggesting a possible approximately linear polarization.This interplay between directional emission and circular polarization can be used in several applications in chiral quantum optics, as anisotropic dis-persions or anisotropic coupled between emitters [46][47][48].

Conclusions
Summing up, we study the spontaneous emission patterns near a two-dimensional photonic crystal slab featuring band structures with straight iso-frequencies and saddle-points.As expected from simplified descriptions in the literature, we are able to find regimes of highly directional emission, a requisite to observe strong super and subradiant effects in two dimensions.However, thanks to our systematic study, we uncover an interesting interplay between such directional emission patterns and the emitter's polarization and position, which has passed unnoticed in other works.In particular, when the emitter is at the center of the unit cell, we find that some of the radiative directions can be canceled by aligning linearly polarized dipoles properly.On the contrary, for out-of-center positions, we find that one can activate three out of the six propagation directions by judiciously choosing the left/right-polarized nature of the emitter's dipole.This can be considered a two-dimensional generalization of the spin-orbit coupling of light appearing in chiral quantum optical setups [46].Our results highlight the importance of carefully considering the electric field's structure in the study of individual and collective spontaneous emission near photonic crystals and pave the way for observing unconventional collective phenomena in such systems.

Fig. 1 :
Fig. 1: Photonic crystal structure.(a) Schematic representation of an atom-like emitter embedded in a photonic crystal slab of a hexagonal lattice with periodic parameter  and thickness  = 0.25, the frequency emission source is   , this photonic structure was proposed in previous work [30].(b) Top view of the photonic crystal slab, the yellow region represents a dielectric media with   = 10.5625, and the white regions show air holes.The unit cell has six smaller holes of radii  = 0.0833 and six bigger holes of radii   = 0.15 with center at a distance  = 0.4 of the unit cell center.(c) Band diagram of the PhC slab in the inverse Brillouin Zone edges (red dashed line in the inset), the gray shaded region shows the frequency region of the van Hove singularity about   = 0.3849 (︀ 2  )︀.(e) and (f) Electric field amplitude for the first band mode for  1 (e) and  2 (f), the vectors  1 and  2 in the irreducible Brillouin Zone are showed as green arrows in (d).

Fig. 2 :
Fig. 2: Schematic representations of the states and the renormalization produced in the energies and emission rates by the mediated interactions.

Fig. 3 :
Fig. 3: Directional emission in six paths.(a) A density plot of the electric field intensity (V 2 /m 2 ) of a point dipole source with  + polarization at the center of the unit cell as depicted in the inset.The dashed lines show four directions to analyze field propagation, as shown in (c).(b) The polar plot of the electric field intensity at a fixed distance of dipole  = 4, the radius of the polar plot represents the field intensity, and the color represents the Stokes' parameter  1 , indicating the linear polarization of the field.(c) Electric field intensity as a function of distance to the source in four different directions ( 1 = /12,  2 = /6,  3 = /3 and  4 = /2 displayed en (a)).It is noticed that the field intensity for the  2 and  4 is bigger than the other in almost all distances .

Fig. 4 :
Fig. 4: Choosing paths with linear polarization.(a) Intensity of the electric field from a linear polarized dipole source with H polarization at the center of the unit cell, as depicted in the inset.(b-d) Polar plot of the electric field intensity at a distance of 4 from the source with dipole orientation at 0 rad [(b)], /6 rad [(c)], and /3 rad [(d)] as indicates the black arrow in each panel; the radial coordinate display the intensity at each angular direction.The color scale of the line represents the  1 parameter, indicating the dominant linear polarization of the field, H polarization (red color), V polarization (blue color), and similar contribution (gray color).

Fig. 5 : 2 .
Fig. 5: Linear polarized dipoles outside the center at r = (0, 0.35).(a) Electric field intensity of an H-polarized dipole in the whole PhC slab.(b) Electric field intensity in the whole PhC slab produced by a dipole moment d = (x + ŷ)/ √ 2. To compare (a) and (b), one has to take into account the color scale in each panel.(c) Electric field intensity at r = (0, 0.35) for all the modes inside the iso-frequency curve.The wave vector of each mode is represented by its angular coordinate   , as depicted in the inset of (d).(d) Stokes parameter  1 at the position r = for all the modes in the iso-frequency curve.

Fig. 6 :
Fig. 6: Circularly polarized emitter at the center of the smallest air holes.(a-b) Electric field intensity is produced for a point source with  + and  − dipole moment, respectively.The dashed line shows the paths that are activated in each case.The inset displays the emitter position inside the unit cell and its polarization.(c-d) Stokes parameter  3 throughout the dashed line as a function of the distance from the source, for  + and  − dipole moment, respectively.