Photon-mediated interactions near a Dirac photonic crystal slab

Dirac energy-dispersions are responsible of the extraordinary transport properties of graphene. This motivated the quest for engineering such energy dispersions also in photonics, where they have been predicted to lead to many exciting phenomena. One paradigmatic example is the possibility of obtaining power-law, decoherence-free, photon-mediated interactions between quantum emitters when they interact with such photonic baths. This prediction, however, has been obtained either by using toy-model baths, which neglect polarization effects, or by restricting the emitter position to high-symmetry points of the unit cell in the case of realistic structures. Here, we develop a semi-analytical theory of dipole radiation near photonic Dirac points in realistic structures that allows us to compute the effective photon-mediated interactions along the whole unit cell. Using this theory, we are able to find the positions that maximize the emitter interactions and their range, finding a trade-off between them. Besides, using the polarization degree of freedom, we also find positions where the nature of the collective interactions change from being coherent to dissipative ones. Thus, our results significantly improve the knowledge of Dirac light-matter interfaces, and can serve as a guidance for future experimental designs.


INTRODUCTION
The Dirac energy spectrum of graphene is the source of many of its extraordinary electronic properties [1]. This has triggered the quest to translate these energy dispersions to other systems, like photonic crystals (PhC) [2], by exploiting the analogy between the electronic and electromagnetic wave propagation [3][4][5][6][7][8]. In this photonic context, these energy dispersions have already been predicted to lead to many non-trivial phenomena such as realizing topologically [3] or pseudo-diffusive [4] photon transport, observing Klein-tunneling [6], or enhancing the Purcell factor over large areas [8]. One of the latest additions to the exciting features of Dirac photonics has been the possibility of obtaining decoherence-free, longrange (power-law) interactions between emitters when many of them couple to these type of structures [9,10].
Long-range interactions are instrumental for many quantum information and simulation applications. For example, they can be harnessed to induce long-distance entanglement [11], to improve the speed of quantum state transfer protocols [12][13][14][15], or to explore (non-) equilibrium phenomena in frustrated spin models [16][17][18][19][20][21][22]. In free-space, photon-mediated interactions are naturally long-ranged (scaling with 1/r 3(1) in the near (far) field), but they are unavoidable accompanied by dissipation [23,24]. Photonic band-gaps can be used to cancel this dissipation [25,26], but at expense of exponentially attenuating the decay of the resulting interactions. Dirac energy dispersions in two [9,10] and threedimensions [27][28][29][30] have been recently pointed out as a way of avoiding this trade-off, showing how they lead to power-law decaying interactions without any associated dissipation thanks to the singular nature of the density of states around the Dirac points.
In the two-dimensional case [9,10], these interactions have been so far predicted to have an overall scaling with a fixed 1/r-decay with the distance, r, between emitters. However, these calculations were done either using toymodel coupled-resonator baths [9], or fixing the emitter positions at high-symmetry points of the unit cell for the case of realistic structures [10], thus limiting the generality of the results. In this work, we go beyond these analyses and derive a general theory for realistic structures that enables us to characterize the emergent photon-mediated interactions for emitters placed at any position of the unit cell. To illustrate its power, we apply this theory to the same Dirac photonic structure analyzed in Ref. 10 and find several remarkable results. First, we study the positions that optimize the interaction strength beyond the one considered in Ref. 10. We characterize both the regions within the dielectric and the air holes, that can be used to optimize the coupling of solid-state [31,32] or natural [33][34][35][36] atomic systems, respectively. Second, we find that certain positions display an effective longer-ranged decay exponent, 1/r γ , with γ < 1, and study the trade-off between the strength of interactions and their range. This tunability is important because it opens the exploration of other long-range interacting models beyond the ones pointed by Refs. 9, 10. Third, by playing with the polarization degree of freedom, we also find positions where the nature of the collective interactions change from being coherent to incoherent ones, which might lead to strong super/sub-radiant effects [37].
The manuscript is divided as follows: First, we introduce the formalism that connects the classical Green The main focus of the work will be to obtain the two-point Green function, G(r1, r2), which characterizes the interaction between emitters. We will restrict to situations where the emitters are placed at positions r1 and r2 = r1 + R, being R is a lattice displacement vector that can be written as a linear combination of the primitive vectors a1,2.
functions of the photonic structure with the effective quantum emitter interactions appearing in a master equation description of the emitters' dynamics. Second, we describe the Dirac-photonic structure that we will consider along this manuscript, and characterize its band-structure and the properties of their eigenmodes around the Dirac point. Then, we explain the basic ingredients of the theory we develop to analyze the effective emitter interactions, that are, the guided-mode expansion technique [38] and the k · p method [39][40][41]. Next, we use our theory to study the positiondependence of the coupling strength, the range, and the coherent/incoherent nature of the photon-mediated interactions appearing in these systems. Finally, we will summarize our findings and point to other possible applications of our method.

PHOTON-MEDIATED INTERACTIONS IN PHOTONIC-CRYSTALS
The goal of this manuscript is to describe the emergent photon-mediated interactions when several quantum emitters couple to Dirac-like light-matter interfaces, as schematically depicted in Fig. 1(a). A suitable formalism to describe these interactions in these dielectric media is macroscopic QED [42][43][44][45][46][47], where the light-matter interactions are expressed in terms of the classical Green function, G(r, r , ω), obtained from the following electromagnetic wave equation: with ε(r, ω) being the permittivity of the medium. Within this formalism, and assuming that Born-Markov conditions are satisfied such that the photonic degrees of freedom can be adiabatically eliminated [42][43][44][45][46][47], the resulting emitter dynamics can be described by the following master equation: where, i) H 0 corresponds to the independent emitters' Hamiltonian. For this manuscript, we will assume that the emitters have a single optical transition from an optically excited state, e, to the ground state g, with frequency ω A and dipole matrix element d, such that H 0 = j ω A σ j ee , using the notation σ j αβ = |α j β| for the emitter operators; ii) H eff represents the unitary (coherent) part of the photon-mediated interactions that for two-level emitters reads: Thus, this term yields coherent excitation exchanges between emitters at a rate J ij , that can be harnessed, e.g., to make SWAP-like gates. Finally, iii) L eff (ρ) describes the non-unitary (incoherent) dynamics induced by the bath that reads: and accounts for both the individual (Γ ii ) and collective (Γ i =j ) dissipation. Despite being non-unitary, these collective dissipative terms yield strong super/subradiance effects [37], which can be harnessed to engineer decoherence-free quantum gates [48,49] or to improve multi-photon generation [50][51][52] and absorption [47] fidelities, among other applications. Remarkably, both the coherent and incoherent terms of the master equation are related to the classical Green function as follows [43][44][45][46][47]: where d i , r i are the optical dipole moment and position of the i-th emitter, respectively. Thus, to know the photonmediated interactions, J ij , Γ ij induced by a photonic media, it suffices to calculate the two-point Green function, G(r i , r j , ω A ) ≡ G(r i , r j ), of the particular structure. Before explaining the method we develop to do it, let us first introduce the particular photonic structure we will use to benchmark our results. In panel (a), we show the numerical calculation of the GME method for dispersion around for a particular path of the first Brillouin zone, depicted in the inset. In panel (b), we compare the k · p approximation with the GME method around the Dirac cone at K, where the light-blue surface correspond to the k · papproximation and the yellow-orange tones correspond to the GME calculations. (c-d) Electric fields in the whole unit cell for an eigenmode in the upper band with momentum k = K + q, with |q| = 0.03|K| and φq = 0. In the two panels we compare the results obtained numerically with GME method (c) and the fields obtained by the k·p approximation (d). The geometric parameters for these calculations are: rs = 0.0833a, r l = 0.15a, l = 0.4a, d = 0.25a; we use 1165 different g vectors and consider 3 guided modes for each g, so we have a basis with total 3480 guided modes.

DIRAC PHOTONIC STRUCTURE
The photonic structure that we will consider is the photonic-crystal (PhC) slab depicted in Fig. 1, which was first introduced in Ref. 10. The slab consists of a hexagonal lattice of GaP (ε GaP = 10.5625) with primitive vectors a 1/2 = ( √ 3/2, ±1/2)a, being a the lattice constant, and with six inner and six outer air holes, of radii r l and r s , respectively. The unit cell of the PhC slab is shown in detail in Fig. 1(b), where l denotes the distance between the center of the unit cell and the inner circles, and d the slab thickness. As shown in Ref. 10, an interesting property of this structure is that it features energetically isolated Dirac-cone dispersions at the high symmetry points K = 2π 3a ( √ 3, 1) and K = 2π 3a ( √ 3, −1) for certain parameters regime that are the ones we will focus along this manuscript.
In Ref. 10 the numerical analysis of the structure was done using a plane-wave expansion method [53]. With the results from these numerical calculations, they proposed an empirical ansatz of the band-structure (ω (n) k ) and their associated eigenmodes at the central position of the unit cell (E k,n (r = 0)) around the Dirac point K. Then, using that ansatz they were able to construct the two-point Green function, G(r i , r j ) for emitters placed at the central point of the unit cell. Here, we will use an alternative approach that consists on a combination of the Guided-Mode Expansion (GME) technique [38] and the k · p method [39][40][41]. The former allows one to find efficient solutions for the eigenfrequencies and eigenmodes of slab geometries, whereas the later will enable us to systematically expand the energy dispersions and eigenmodes around the Dirac point, and find the associated electromagnetic field at all positions of the unit cell. This will be the key advantage of our approach, since we will be able to compute the two-point Green function, and thus the emergent photon-mediated interactions, for emitters placed at any point of the unit cell.
Let us now briefly explain the main steps of our theory, and apply it to characterize the photonic structure of Fig. 1. The first step consists in using the GME technique [38]. This method initially finds the solutions of the homogeneous slab waveguide with an effective permittivity, that we denote by h g,µ (r, z). Using Bloch-theorem we can rewrite these solutions as h g,µ (r, z) = e ik·r U g,µ (r, z), separating the part e ik·r given by Bloch theorem, from U g,µ (r, z), that provides the real space distribution within each unit cell associated to each guided mode with momentum g, and which has the same periodicity than the original PhC lattice. Here, g = k + G is a combination of k inside the first Brillouin zone, G a vector of the reciprocal lattice with information about the periodicity, and µ denotes the different modes of the effective homogeneous waveguide for a given g, which is related with the z-quantization of the modes. Note also that here we use r to express a vector in the xy plane and we will write the dependence with the z coordinate when needed. Using that form for the eigenmodes of the homogeneous slab, the magnetic field of the complete photonic structure is expanded in terms of these guided modes of the homogeneous waveguide (see Supporting Information), as follows where c n (k + G, µ) are the coefficients of expansion, and n denotes the different modes of the PhC slab that appear with frequency ω (n) k . Inputting this expansion into the eigenvalue problem of the Maxwell equations, one arrives to a linear matrix eigenvalue problem that can be solved numerically by truncating the size of the guided mode basis. Since we are interested in studying emitters placed at the z = 0 plane of the slab, we can restrict to the study of transverse-electric TE-like (even) modes, because transverse-magnetic TM-like ones have an electric field zero at z = 0 and thus will not couple to the emitter.
In Fig. 2(a) we plot the band-structure of the TE-like modes of the structure of Fig. 1 using the following geometric parameters: r s = 0.0833a, r l = 0.15a, l = 0.4a, d = 0.25a, and calculated with the GME method using a truncation for the homogeneous waveguide basis of 3480 guided modes. There, we observe how indeed energetically isolated Dirac-cone dispersions at the K-point appear, that will be the main focus of this article. The next step consists in making a perturbative expansion around the K-point, i.e., k = K + q, for small q, using the k · p method. The key idea of the perturbative expansion is to use the semi-analytical solutions obtained from the GME method at the K point as the basis to express the problem at k = K + q. For that, we use agin the fact that the magnetic field within the PhC also satisfies Bloch theorem such that we can rewrite it also as H k,n (r, z) = e ik·r u k,n (r, z), where u k,n (r, z) is the Bloch periodic function which contains the dependence of the electric field within the unit cell. Using the GME results, we can obtain the basis of periodic function {u K,n (r, z)} around K as follows: (8) with which we can write the magnetic field for momenta around the Dirac point, k = K + q, as Putting this expansion into the eigenproblem of the magnetic field, one arrives in the following equation where the matrix elements H l,j can have terms of order q 0 , q 1 to q 2 , as explicitly shown in the Supporting Information. Since the Dirac points K ( ) have degenerate modes, we only consider the contribution of the two first bands (l, j = 1, 2). Additionally, as we find that the energy dispersion is approximately linear, we can also neglect the order q 2 terms. With these considerations the general eigenvalue problem of Eq. (10) reduces to a 2 × 2 matrix eigenproblem: where ∆λ ± = . The subindex ± is related to the upper and lower bands of the cone, and ω D corresponds to the frequency at the Dirac point, which we define numerically as the average between ω (1) K ( ) and ω (2) K ( ) . These particular relations between the matrix elements result from the crystal symmetry called "deterministic degeneration" [54]. Details about P l,j and its calculation are given in the Supporting Information.
Solving the simplified eigenvalue problem of Eq. (11), we indeed obtain a linear dependence with |q| of the eigenvalues, as ∆λ ± = ±2 v c ω D c |q|, with v being the group velocity at the Dirac points, with which we can approximate the energy dispersion of the two bands as follows: In Fig. 2(b) we plot these analytical approximations (in shaded blue surface) together with the numerically obtained energy bands using the GME method (orange surface), showing a good agreement between the two. Interestingly, from the simplified eigenvalue problem obtained through the k · p approximation, we can also obtain the following magnetic and electric field expansions around the Dirac points where k 0 can be K or K ; and the subindices 1, 2 indicate the two degenerate modes at K ( ) . The parameters ξ ± and η ± have the following values for K ( ) : where the symbols ± and ∓ are for K/K , respectively, tan(φ q ) = q y /q x , and the phases δ K ( ) = ±π/6. In Figs. 2(c-d) we plot a comparison of the electric field amplitude corresponding to a particular k value of the upper band along the whole unit cell obtained numerically through the GME method (c) and semi-analytically using the k · p-approximation (d). There, we can see how the approximation captures indeed very well the emergent physics along the whole unit cell. Finally, let us also note that despite the different appearance of the expressions of the electric field in Eq. (14) with respect to the empirical ansatz proposed in Ref. 10 for r = 0, they reproduce the same physics if we restrict our approximated fields to this position, as we explicitly show in Supporting Information.

CONSTRUCTING GREEN FUNCTIONS
Now, we will use the analytical approximations of Eqs. (12)- (14) to construct the Green function of the problem, as also did in Ref. 10. However, the advantage of our method is that the electric field expansion of Eq. (14) is valid for any position of the unit cell, not only for the center r = 0, and thus, we will be able to calculate G(r 1 , r 2 ) in a more general fashion. For this work, we restrict to the situation in which the position of the two emitters differs only by primitive lattice displacement, that is, Thus, we only target to calculate: G αβ (r 1 , r 1 + R) := G αβ (r 1 ; R). This Green function can be calculated integrating the momentum-space Green function as follows: with: where E (n) p,α (r i ) is the electric-field α-component at the r i position associated to the eigen-energy ω (n) p of the n-th band.
Note that due to the periodicity of the PhC (see Eq. (14)), E p ), and the integral in Eq. (19) is mostly given by the contributions around the K ( ) -points. Then, we can use the analytical approximations of Eqs. (12)- (14) to obtain the Greenfunction components: where j (r) are the first kind Hankel function of order j, and are coefficients defined in terms of the direction of R and the electrical fields at position r 1 for the two first modes of K ( ) . More details about the explicit form of these coefficients can be found in the Supporting Information. We do not consider G αz because it is strictly zero for emitters placed at the symmetry plane z = 0. Finally, let us note that the Green functions associated to the circular polarized components can be reconstructed from Eq. (21) as follows:  In the next section, we analyze these functions in detail, putting special emphasis on the dependence of G αβ (r 1 , R) at different places of the unit cell, since it is the main strength of our method.

DIRAC-PHOTON-MEDIATED INTERACTIONS ALONG THE WHOLE UNIT CELL
Let us start by considering the photon-mediated interactions of emitters with linearly polarized optical transitions oriented along one of the Cartesian components (x,ŷ). As aforementioned, these interactions are given by the Green functions of Eq. (21), that are expressed as sums of several decaying terms scaling with different power-laws. In previous works [9,10] the dominant contribution was shown to be given by a 1/r-decay law. However, we find that when the position r 1 of the emitters is varied, the interference between the different terms can give rise to even longer-ranged interactions at certain points. This is observed in Fig. 3(a), where we plot |G xx (r 1 ; R)| as a function of R for a particular primitive lattice direction (other components and directions lead to qualitatively similar conclusions) for three different emitter positions in different colors, highlighted with arrows in Fig. 3(b). At all emitter positions chosen, |G xx (r 1 ; R)| have an oscillatory and power-law decay behaviour. However, the decay range is faster in some positions with respect to the others. To make this more evident, we take only the maximum value of the oscillations within each period and plot them in the inset of Fig. 3(a) in logarithmic scale, where one clearly observes that they follow a different power-law behaviour.
To make a more detailed analysis of the change of this exponent we make a numerical fit of the envelopes of |G xx (r 1 ; R)| to a power-law ∝ 1/|R| γ at all positions of the unit cell r 1 and plot the results of the fitting in Fig. 3(b) in a color scale. In this figure, the white color denotes regions where the 1/r behaviour dominates (γ = 1), whereas red denotes scalings with longer-range decays (γ < 1). There, we observe how indeed most of the regions display the expected 1/r behaviour predicted in previous works. However, there also appear other regions where the interference between the different terms lead to longer-ranged interactions.
Apart from the range of the interaction, another very relevant magnitude is their strength. This strength also depends on the emitter position r 1 , since the mode function profile also changes along the unit cell. A way of characterizing this strength is by plotting the absolute value of the photon-mediated interactions between nearest-neighbouring atoms placed at different r 1 positions, |G xx (r, r + a 1 )| ≡ |G (1) xx (r)|. This is what we plot in Fig. 3(c), multiplying it by the dielectric index ε so that the strength at the air/hole regions appear on a similar color scale. From this figure, we can make two important observations: first, the optimal position to couple the emitter is not at the center of the unit cell. For emitters within the dielectric, the largest coupling strength is obtained at the regions between the air holes. For emitters lying in the holes, the coupling is maximum at regions very close to the interface. However, trapping atoms close to surfaces is generally challenging due to Casimir-Polder forces [26,55]. Thus, in that case it would still be easier to place the atoms at the center of the holes, even if the coupling strength is smaller. The second important observation is that there is an apparent trade-off between the range and strength of interactions, since the red regions of panel (b) appears to be whiter (no coupling) in panel (c). We will explore this trade-off in more detail in Fig. 5.
Finally, let us note that in Figs. 3(a-c) we plot the absolute value without differentiating the contributions from the real and imaginary part of the Green Function. However, as we introduced in the first part of the manuscript,   both terms give rise to very different quantum dynamics: the real part leads to purely coherent exchanges (J ij ), whereas the imaginary part (Γ ij ) yields super/subradiant effects. Thus, in Fig. 3(d) we characterize which term dominates by plotting the ratio between the real and imaginary part of the nearest-neighbour position: As shown in the legend accompanying the panel we use a logarithmic colorscale for W 1 (r) where the red (blue) color denotes W 1 (r) > (<)1, while white denotes the transition where W 1 (r) = 1. From the figure we can clearly observe how linearly polarized emitters lead to decoherence-free photon-mediated interactions (W 1 (r) > 1) along the unit cell. This is in accordance to the results of the literature [9,10] and it is expected because of the vanishing density of states of the photonic bath around the Dirac point, i.e, D(E) ∝ |E|.
Let us now repeat this analysis with emitters with cir-cularly polarized transitions. For the sake of illustration we focus again on a particular component of the Green Function, i.e., G σ+σ+ , although the conclusions can be extrapolated to the other components. We start by plotting again |G σ+σ+ (r 1 , r 1 + R)| in Fig. 4(a) at several representative positions r 1 , with different power-law scaling. This justifies again making the numerical fit to a powerlaw ∝ 1/|R| γ at all unit cell positions, whose results we plot in Fig. 4(b). There, we observe how although the dominant exponent is the expected γ = 1, there are several regions with longer-ranged interactions, as it occurs for the linearly polarized transition. Then, in Fig. 4(c) we plot the first-neighbour coupling strength multiplied by the permittiviy of the material, i.e., ε|G 1 σ+σ+ (r 1 )|. Differently from the linearly polarized case, the center of the unit cell is one of the least efficient places to couple the emitter. Again, the more favorable regions are the dielectric parts within the holes, although there also appear other places around the central part of the unit cell. At the holes, the pattern is very similar to the one of linear polarizations of Fig. 3(c). The main difference with respect to the linearly polarized emitters is the behaviour of the real over imaginary part ratio, W 1 (r), defined now for G (1) σ+σ+ (r 1 ), which we plot in Fig. 4(d). There, we observe how the photon-mediated interactions in most of the unit cell are dominated by the collective dissipative terms (blue regions), differently from the linearly polarized transitions (see Fig. 3(d)). This is a consequence that for certain atomic positions of this structure Re [G xy ] = Re [G yx ], which is known to lead to decay in circularly polarized transitions [56][57][58][59]. This change of the coherent nature of the interactions with polarization is something that was not captured in previous analysis because they either neglected polarization effects [9] or placed the emitters in positions where G σ+σ+ ≈ 0 [10].
One important result from the previous analyses is the possibility of increasing interaction range (decreasing γ) by placing the emitter at different positions (see panel (b) of Figs. 3, 4). However, as we also see in panels (c) of the same figures, this increase comes at the price of reducing the absolute coupling strength. Since the color scales does not let us appreciate in detail this trade-off, in Figs. 5(a-b) we plot both γ (black) and ε|G (1) αβ (r)| (red) for the linear and circular polarization, respectively, and following a particular path of the unit cell depicted in the insets of this panels. There, we notice that indeed the minima of γ coincide with positions where |G  observe in spite of this trade-off there are regions within this path where the regions with γ < 1 leads to an overall better coupling at long distances than the one of γ = 1.

CONCLUSIONS AND OUTLOOK
Summing up, we have developed a semi-analytical theory based on the Guided-Mode Expansion and k · p method to calculate the photon-mediated interactions in Dirac light-matter interfaces. An important advantage of our theory with respect to existing ones is that it enables us to calculate such interactions when emitters are placed at all positions of the unit cell. We benchmark our theory on a particular photonic structure, and find that it is possible to tune the interaction range of the emergent interactions by placing the emitter at different positions of the unit cell. Besides, we find the position that optimize the interactions between emitter at short and long-distances. Thus, we believe our theory and results can become a useful guide for future experimental designs of such Dirac light-matter interfaces.
AGT acknowledges support from CSIC Research Platform on Quantum Technologies PTI-001 and from Spanish project PGC2018-094792-B-100(MCIU/AEI/FEDER, EU). EPNB thanks financial support from the "Programa de becas de excelencia doctoral del bicentenario -MINCIENCIAS 2019". HVP gratefully acknowledges funding by COLCIENCIAS under the project "Impact of phonon-assisted cavity feeding process on the effective light-matter coupling in quantum electrodynamics", HERMES 47149. AGT acknowledge discussions with P. A. Huidobro about the decay of circularly polarized transitions.

Supplemental Material: Photon-mediated interactions near a Dirac photonic crystal slab
In this Supporting Information, we provide more details on: i) the Guided-Mode Expansion method to calculate the eigenfrequencies and eigenmodes of the Dirac photonic structure; ii) the k · p method to approximate the eigenfrequencies and eigenmodes around the Dirac cones; iii) the derivation of the Green tensor by analytical integration using the expressions obtained by the k · p method; iv) additional results for other components of the Green function.

Guided-Mode Expansion method
Time-independent electromagnetic waves satisfy the following eigenvalue problem where the eigenvalues are ω 2 c 2 , and the eigenfunctions are the magnetic field H(r, z). Additionally, in materials with discrete symmetry translation, like photonic crystals (PhC), the eigenfunctions fulfill the following Bloch condition where k is the pseudo-momentum, and u k (r, z) is a periodic vector function with the same periodicity of the material, containing the real space dependence of the eigenfunctions within the unit cell. Photonic crystal slabs can always be schematically represented as in Fig. SM1(a): they are structures with three different periodic materials in each region: z < −d/2, −d/2 ≤ z ≤ d/2 and z > d/2. To solve Eq. (SM1) in these structures, one can use the Guided-Mode Expansion (GME) method [38], which initially considers a homogeneous slab waveguide with an effective permittivity given by the mean value of the permitivityε j in the xy-plane for each region in z, as shown in Fig. SM1(b). To ensure the possibility of total internal reflection inside the region 2, the structures have to satisfy the conditionε 2 >ε 1(3) . This homogeneous problem has two types of solutions: the guided modes, that are quantized and evanescent in the z direction at the regions |z| > d/2; and the radiative ones, that are a set of continuous modes and oscillating in z direction. The GME method use the former, plus Bloch theorem, to compose a new basis for solving the complete PhC slab problem. The basis can be separated into two orthogonal families, transverse electric (TE) and transverse magnetic (TM) modes. The TE modes can be written as: and TM modes: g,µ and D (j) g,µ are coefficients determined by the boundary conditions between the three spatial regions in the z-direction and by a normalization condition [38]; g is the projection of the wavevector in the xy-plane, z is the unitary vector in z-direction, S = √ 3a 2 /2 is the area of the unit cell, g = |g|,ĝ = g/g, andê g =ẑ ×ĝ. The constants χ (j) g,µ and q g,µ are related with g and the frequencies modes by The eigenfrequencies of the TE polarized modes are obtained by solving the following transcendental equation: whereas the TM polarized frequencies are given by: For convenience in the derivation of k · p-approximation, we can introduce the guided periodic function, U g,µ (r, z), of each homogeneous guided modes (h g,µ (r, z)) as follows: h g,µ (r, z) = e ik·r U g,µ (r, z) → U g,µ (r, z) = e −ik·r h g,µ (r, z) , Using these guided modes for the homogeneous slab, we can expand the magnetic field of the PhC as where we write g as the sum of k a vector in the first Brillouin zone and G a vector of the reciprocal lattice, such that we have the information about the periodicity of the structure already incorporated in the expansion. Here, the index µ runs over the different homogeneous guided modes, and we also introduce the index n related to denote the different modes (bands) of the PhC slab. Then, using the expansion of H k,n of Eq. (SM10) and the orthogonality of the guided modes, the eigenvalue Eq. (SM1) is rewritten as 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 H g,µ;g ,ν read where V c is the integration volume corresponding to the unit cell in the xy-plane and (−∞, ∞) in the z-direction. An explicit expression for the matrix elements in terms of the guided modes can be found in Ref. 38. Eq. (SM11) is, in principle, an infinite linear matrix eigenproblem that must be truncated to be solved numerically. Solving this truncated eigenproblem, one can calculate the eigenfrequencies ω (n) k and the eigenvectors containing the coefficients c n (k + G, µ) that expand the magnetic field of the PhC modes. Due to the reflection symmetry at the horizontal plane at z = 0 of our structure (z → −z), it is possible to classify the solutions into two families, odd (TM-like) and even (TE-like) eigenmodes. We are interested in the TE-like modes because they have non zero electric fields at z = 0 where we put the emitters, so we only take into account the modes with even symmetry under reflection in z. k · p method for Dirac cones In the main text, the k · p method is used to obtain analytical expressions of ω (n) k , H k,n (r, z) and E k,n (r, z) for momenta k around the Dirac points K ( ) . Here, we provide more details about this procedure.
According to Bloch theorem, the solution to the magnetic field at the K point has the form of H K (r) = e iK·r u K (r). Using the GME method, we have an expression of the periodic function {u K,n } which reads: We use the set {u K,n } as a basis to expand u k,n for k = K + q, with which the magnetic field can be written as: By putting this magnetic field in Eq. (SM1), we obtain the following eigenproblem [39] j where, The H l,j matrix elements have in general terms of order q 0 , q 1 and q 2 , such that solving Eq. (SM15) yields a perturbative solution of the eigenvalues and eigenmodes. To solve this equation numerically. we use a finite set of u K,j (r), and obtain the {C j } coefficients. With them, we can reconstruct the magnetic field of the PhC slab using Eq. (SM14), and finally, obtain the corresponding electric field using Maxwell equations.
For the particular Dirac-cone scenario, this general perturbative eigenvalue problem can be further simplified. First, because for Dirac cones around the K ( ) points, we have only two degenerate eigenvalues at the two first bands. Thus, we only need to consider two orthonormal modes from this eigenvalue to solve a perturbative degenerate problem. Second, because the dispersion at these symmetric points is approximately linear, such that we only include the terms of H l,j with zero and first order in q. With these two considerations, Eq. (SM15) is simplified to 2 j=1 q · P l,j C j = ω 2 k,n − ω 2 K,n where P l,j = i(−p l,j + p * j,l ) + 2q lj K − (Q l,j + Q T l,j ) · K. Besides, as a consequence of the deterministic degeneracy [54], there are particular relationships for the P l,j 's, that are, P 2,2 = −P 1,1 and P 2,1 = P * 1,2 , which are approximately satisfied in numerical calculations. In this way, Eq. (SM20) has the following matrix form q · P 11 q · P 12 q · P * 12 −q · P 11 where the eigenvalues are ∆λ ± = ± ω 2 k± −ω 2 0 c 2 , and the subindices ± are related to the upper and lower bands of the cone. Then, solving the 2 × 2 matrix eigenproblem, we obtain ∆λ ± and (ξ ± , η ± ). The former gives us the correction to the squared frequencies in a linear way where m is the slope of the cone, which depends of the geometric parameters of the slab. Around the Dirac cones m|q| 1, such that the square root of Eq. (SM22) can be approximated by: where v is the group velocity of the Dirac cone. Using this velocity, we define m = 2 v c ω D c and the eigenvalues as The coefficients (ξ ± , η ± ) give us the contribution of each orthonormal u K,j (r) to the magnetic fields (H k,n (r)). With these coefficients, the magnetic fields are written as follows Remembering that k = K + q and H K,j (r, z) = e iK·r u K,j (r, z), we can finally express The electric field can be derived from the relation For Dirac cones, these coefficients can be expressed in trigonometric functions, as shown in the main text: where the symbols ± and ∓ are for K/K , tan(φ q ) = q y /q x , and the phases δ K ( ) = ±π/6. Additionally, we can use guided modes to express and calculate the vectors P l,j and the matrix elements H l,j . Using the expansion of u K,j (r, z) in terms of the guided periodic functions U g,j (r, z), we obtain the following expressions.
where, The elements P g ,ν;g,µ and Q g ,ν;g,µ are obtained by integration of guided modes with TE or TM polarization, and each polarization case have different analytical expressions, so we add the superscripts (TE-TE), (TE-TM), (TM-TE) and (TM-TM), for example: The explicit expressions for each element can also be analytically calculated but we do not write them here because they are so large that do not add any further physical insight.
To conclude this section, let us benchmark the approximated electric fields obtained through the k · p expansion doing two sanity checks. First, let us compare our results (Eq. (SM26)) with the ones derived by Perczel and Lukin for the central position of the unit cell, which read: for k around K ( ) . Since the direct comparison between formulas is not obvious, in Fig. SM2 we plot the dependence of the electric field components with the angular variable φ q , calculated by k·p-method and by the ansatz of Eq. (SM36).
To do it, we fix the position to r = 0 and |q| = 0.03|K|. We plot in blue (red) the x(y)-component of the electric field according to Ref. 10 (in empty circles) or to k · p-approximated expressions (in solid-lines) showing a very good agreement for both the K (panels (a-b)) and K (panels (c-d)) points. Apart from this comparison, it is important to ensure the capability of our approximation to capture the electric field dependence at other positions of the unit cell. A way to benchmark it consists in fixing a particular q and calculate the field mode distribution along the whole unit cell using both a full numerical approach (GME) and the approximated expression. This is what we show in Fig. SM3, where we compare in panels (a)-(b) [(c)-(c)] both approaches for k close to the K [K ] points. In the two panels we compare the results obtained numerically with GME method, panels (a) and (c), and the approximated fields obtained by the k · p approximation, panels (b) and (d). The geometric parameters of the PhC slab are the same as in Fig. SM2.

Green functions: analytical integration using k · p method
The calculation of the two-point Green function is made using the definitions: and We are interested in emitters frequencies very close to the Dirac point, so we make the following approximations: • We only consider the contributions of the first two modes (n = 1, 2), that are, the two bands that compose the Dirac cone.
• We only consider k's close to the high symmetry points K ( ) ; these k's are delimited by a circular cutoff around K ( ) , i. e., |k − K ( ) | ≤ q c . In this region, the dispersion relation is linear with q = k − K, as the k · p method shows, and the electrical field for the two bands are given by Eq. (SM26).
• As ω A ≈ ω D , we consider that • We are interested in points separated by lattice displacements, i.e., r 2 = r 1 + R, so the fields at the two points are related by the Bloch theorem as follows: E k,± (r 2 ) = e ik·R E k,± (r 1 ) Using these approximations, the integral of Eq. (SM37) are written as where, the subindex K ( ) in the integral denotes that the integrals are made around the Dirac point at K ( ) . Here, we have considered a change of the integration variable p = K ( ) +q and written q = (q cos(φ q ), q sin(φ q )) in polar coordinates. The limits for the integration variables are 0 ≤ q ≤ q c and 0 ≤ φ q ≤ 2π. Now, we use the explicit form of the electrical field as found in Eqs. (SM26)-(SM30) to write the following product of the components: Considering the trigonometric expression for the coefficients ξ ± and η ± (Eqs. (SM27)-(SM30)), we have three types of terms: In the case of E (K,+) * q,α (r 1 )E (K,+) q,β (r 2 ), we group the terms in the following way: Putting this product in the integral I (K,+) α,β , we obtain three types of integrals: where, and the integrals are: To evaluate the integral, we use the polar representation of R = R(cos(φ), sin(φ)). Additionally, we do the change of variableq = vq/δ A , with δ A = ω D − ω A . When δ A is small, the upper limit for the integral in this new variable is q c = vq c /δ A 1; then, we can extend this limit to infinity.
is the second kind Bessel function, and H j (x) is the Struve function of order j. For the lower band, we have the following contributions where the integrals are: J j (x) is the first kind Bessel function of order j. These integrals have poles, so we add a small imaginary part in the denominator of the integrand and make the integration using the Sokhotski-Plemelj theorem for the real line: where P represents the Cauchy principal value, and a < 0 < b. The sign +(−) in the denominator finally results in the first (second) kind Hankel function, representing the outgoing (incoming) waves. In contrast with Ref. 10, we have chosen the first kind Hankel function, so that the imaginary part at the R = 0 position is positive. Similar to K, we can write the contribution for K in terms of three integrals: where the explicit forms of a K αβ , b K αβ , C K αβ and I ± K ,m are analogous to the case of K, the only difference is to change K to K .
Adding all the contributions we obtain where H (1) j (x) is the first kind Hankel function of order j, and the information of the electric field modes are given by the A and B coefficients which read: remembering δ K ( ) = ±π/6. Additionally, K and K are conected by a rotation and inversion, so the fields have the following symmetry relations: using these symmetry relations, it is possible to link the a K αβ (r 1 ), b K αβ (r 1 ) and c K αβ (r 1 ), with their analogous in the symmetry point K, as follow: In the numerical calculation, these relations are not exactly satisfied, especially close the discontinuities in the dielectric index due where the convergence of GME method is more challenging. This induces a small numerical error when calculating the Green function, that leads to unphysical results, like having a non-zero value of the imaginary part in the limit R → 0 (single point Green function) when it should be exactly 0 due to the vanishing density of states. To circumvent this problem, we do not calculate the a K αβ (r 1 ), b K αβ (r 1 ) and c K αβ (r 1 ) directly from the GME method, but instead we use the results of a K αβ (r 1 ), b K αβ (r 1 ) and c K αβ (r 1 ) to calculate their analogous at K using the Eqs. SM68-SM70.

Dirac-Photon-mediated interactions: additional results
As shown in Eqs. (SM63) the dependence of the Green function on the distance between emitters stems from the combination of Hankel functions H (1) j (x). In the case where the emitters are tuned close to the Dirac point, i.e., δ A → 0, these functions display an asymptotic scaling with: for |x| 1, where γ E is the Euler constant. Depending on the emitters' position, the coefficients A and B accompanying will have different values, and thus lead to a different combination. As shown in the main text for G xx (r 1 ; R) and G σ+σ+ (r 1 ; R), the resulting spatial decay can be always shown to be oscillating and with a power-law decay. Thus, the Green function can be fitted to a power-law envelope function ∝ 1/|R| γ , where the γ gives us a quantifier of the effective range of the interactions emerging from the combination of these Hankel functions. Additionally to the range of interactions, another relevant magnitude is their strength; we quantified this by the absolute value of the two-point Green function between two nearest-neighbouring emitters, |G αβ (r 1 ) (1) |. Finally, it is essential to characterize their coherent (incoherent) nature; to do this, we use the quantifier W 1 (r 1 ) defined in the main text.
Here, we display some additional results for other components that are not shown in the main text. First, we present the results for linear polarization, and second for circular polarization.

Linear polarization
For linear polarization, we show the interaction parameters for the components G yy in Fig. SM4 and G xy in   observe that γ is close to 1 in almost all the unit cell for G yy and G xy , except for some regions in red that have γ < 1 with longer-ranged interactions. In the case of the component G xy , some small regions where the γ is slightly greater than one are in light-blue color in Fig. SM5(a). Remembering the trade-off between the range and strength of the interactions discussed in the main text, we also characterize the strength in Figs SM4(b) and SM5(b) for both components. We find similar results as in the main text: the regions with longer-range interactions have less interaction strength. However, it is possible to choose the position that tunes an interaction with γ < 1 and relevant strength. We also observe that positions with maximum strength are between the air holes for emitters inside the dielectric and close to the edge of holes for emitters lying in the air. For cross-polarization (G xy component), the region around the center has a good strength similar to the maximum at the optimal position. Finally, in Figs SM4(c) and SM5(c), we characterize the nature of the interactions by plotting in logarithmic scale W 1 , finding that the interactions are coherent along most of the unit cell, as also occurred for the G xx component shown in the main text.

Circular Polarization
For the case of circularly polarized transitions, we also find a trade-off between the range and strength for G σ−σ− and G σ+σ− components. This is illustrated in Figs. SM6(a-b) and SM7(a-b), for each polarization, respectively. Although we find the optimal position very close to the holes for emitters inside the dielectric, some regions close to the center have large coupling strength for both components. As it occurs for the G σ+σ+ component shown in the main text, the center is not an optimal position to couple emitters with σ − polarization, since the strength at this position is approximately zero. This does not occur for the cross polarization situation, i.e., G σ+σ− component. In  the case of emitters in air holes, the maximum strength is found very close to the interface air-dielectric for all the components in both polarization bases. One important difference with respect to the linearly polarized cases is that the nature of the interactions can change from being coherent to incoherent, as shown already for the main text for G σ+σ+ . This also occurs for the other circularly polarized components as shown in detail in Figs. SM6(c) and SM7(c).