Quantum electrodynamics in anisotropic and tilted Dirac photonic lattices

One of the most striking predictions of quantum electrodynamics is that vacuum fluctuations of the electromagnetic field can lead to spontaneous emission of atoms as well as photon-mediated interactions among them. Since these processes strongly depend on the nature of the photonic bath, a current burgeoning field is the study of their modification in the presence of photons with non-trivial energy dispersions, e.g., the ones confined in photonic crystals. A remarkable example is the case of isotropic Dirac-photons, which has been recently shown to lead to non-exponential spontaneous emission as well as dissipation-less long-range emitter interactions. In this work, we show how to further tune these processes by considering anisotropic Dirac cone dispersions, which include tilted, semi-Dirac, and the recently discovered type II and III Dirac points. In particular, we show how by changing the anisotropy of the lattice one can change both the spatial shape of the interactions as well as its coherent/incoherent nature. Finally, we discuss a possible implementation where these energy dispersions can be engineered and interfaced with quantum emitters based on subwavelength atomic arrays.


I. INTRODUCTION
The extraordinary electronic properties of graphene are directly related to the emergence of isotropic Diraccone energy dispersions in the electronic band-structure of these systems [1]. These properties have sparkled the interest in finding other materials where these energy dispersions appear [2,3], as well as in finding ways of modifying them to tune the emergent behaviour. For example, by introducing tunneling anisotropies, the Dirac cones can be displaced in reciprocal space, and even be annihilated by making two of them overlap [4][5][6][7][8][9][10][11][12][13][14][15][16][17]. This overlap generates a new type of energy band touching linear in one direction and quadratic in the other, labeled as semi-Dirac point, which leads to unusual magnetic field dependence of electron transport [7] or anisotropic polariton transport [18]. Another example is the case of tilted Dirac dispersions that has been predicted to appear in quinoid and hydrogenated graphene [19][20][21][22], and which can be harnessed for valley filtering [23] or generating photocurrents [24]. Besides, when the tilting reaches a critical value, its associated Fermi surface changes from being a point to a line becoming what have been labeled as Type II and III Dirac points [25][26][27][28][29][30][31][32][33][34][35][36], which are believed to be instrumental to enhance the superconducting gap [37], to probe flat-band physics [38] and quantum chaos [39], or even study analogue black holes in solidstate environments [30,40].
Triggered by these exciting properties, the interest in engineering Dirac energy dispersions has expanded to fields beyond the electronic realm. One paradigmatic example is the case of photonics, where such Dirac energy dispersions have been proposed and engi- * a.gonzalez.tudela@csic.es neered to obtain unconventional photon transport properties [18,33,36,[41][42][43][44][45][46][47][48], among other phenomena. One of the latest frontiers of the field is to consider the interplay of these structured photons with quantum emitters, since they are expected to modify substantially quantum optical phenomena such as spontaneous emission or photon-mediated interactions. These changes have already been explored in two recent works for isotropic Dirac photons [49,50]. There, it was shown how a single emitter tuned to the Dirac point can display a nonexponential relaxation in spite of the vanishing density of states at this point [49]. Besides, this unconventional relaxation was linked to the emergence of a power-law localized photonic mode around the emitter, labeled as quasi-bound state [15,51,52], which can ultimately mediate power-law interactions between emitters with no associated dissipation [49,50]. These long-range interactions are especially attractive because they can be harnessed for quantum information and simulation purposes [53][54][55][56][57]. Since more complex Dirac-energy dispersions have started being designed and built within the photonic context [33,36,47,48], a timely question is to understand their impact in these quantum optical phenomena when quantum emitters interact with such photonic lattices.
In this work, we target precisely this problem by considering a system where quantum emitters interact with a honeycomb photonic lattice with a controlled degree of anisotropy and tilting of the Dirac points. For the case with anisotropic Dirac points we observe how indeed the spatial shape of the photon-mediated interactions can be modified without losing neither its long-range nor its decoherence-free character, as it occurs with the case of 3D Weyl points [58][59][60]. Something similar occurs for the tilted case, until a critical value of the tilting when the density of states becomes finite and the resonant nodal lines lead to directional radiation patterns. The latter will translate into anisotropic collective dissipative terms arXiv:2106.10743v1 [quant-ph] 20 Jun 2021  Figure 1. a) Scheme of the system: one or many two-level emitters are coupled locally with strength g to a photonic lattice described by a bipartite lattice with primitive vectors c 1/2 (b) Hopping structure of the bipartite bath: we consider nearest (green) and next-to-nearest (orange) hoppings with overall strength J1,2, respectively. The strength of the nearest (and next-nearest neighbour) hoppings along one direction can be modified through an anisotropy parameter β1 (β2), respectively. when many emitters couple to the bath, as the ones found in photonic Van-Hove singularities [61][62][63], that will lead to strong super/subradiant effects [64]. We illustrate all these effects using a minimal tight-binding model of the photonic bath where these energy dispersions appear, that allows to extract numerical and analytical insight of these phenomena. However, we also finally provide a possible way of probing these effects using a combination of subwavelength atomic arrays [65,66] and additional impurity atoms [67][68][69].
The manuscript is structured as follows: in Sections II and III we describe the system under study and the theoretical framework that we will use to characterize it, respectively. In Section IV, we review the results known for isotropic Dirac dispersions, and then in Sections V and VI we analyze the change of the photon-mediated interactions with the degree of anisotropy and tilting of the bath energy dispersion. Then, in Section VII, we analyze the band-structure of several subwavelength atomic arrays configurations where some of these energy dispersions appear. Finally, we conclude and summarize our findings in Section VIII.

II. SYSTEM
The system that we consider along this manuscript (see scheme of Fig. 1(a)) consists of a collection of N e quantum emitters that, for simplicity, we consider as two-level systems with a ground, g, and an excited state e with an energy difference ω 0 . The corresponding Hamiltonian then reads (assuming ≡ 1 for the rest of the manuscript): where we use the notation σ j αβ = |α j β| j for the operator acting on the j-th emitter. Through their optical transition, these emitters couple to a photonic bath site, described by a bosonic annihilation (creation) operator O ( †) n , through the standard light-matter Hamiltonian: where g j denotes the strength of the coupling for the j-th emitter, and n j denotes the position of the bath where the emitter couples to. For simplicity we will assume this coupling to be emitter-independent, g j ≡ g. Note also that in Eq. (2) we are implicitly assuming that the light-matter coupling is i) local, since the emitter can only couple a single bosonic mode at a given position; and ii) excitation-conserving, since we are neglecting the counter-rotating terms (O nj σ j ge ,. . . ). Both approximations are generally a good description in the optical regime [70], or even an exact one in the case of simulated light-matter couplings with matter-waves [71][72][73], that is one of the most promising implementations to observe these effects.
The non-trivial part of the problem stems from the photonic bath. To be able to extract analytical and numerical insight of the physics, it is convenient to adopt a discrete photonic lattice model where the bath is described as an array of coupled cavities. Since we want to describe Dirac-like dispersions the array is disposed in a honeycomb geometry as depicted in Fig. 1(a). This array can be described by a bipartite lattice with bosonic operators a n , b n for the A/B sublattices, respectively. We label the positions n = (n 1 , n 2 ) with n i ∈ Z in terms of the primitive lattice vector displacements, i.e., n → 2 i=1 n i c i , with c i being the primitive vectors of the lattice c 1/2 = 3/2ê x ± √ 3/2ê y (taking the lattice constant as the unit of length). Imposing periodic boundary conditions, we can define O k = 1/ √ N n e −ik·n O n , for O = a, b, being N the total number of sites of each sublattice, and k → 2 i=1 k i d i with d i being the primitive vectors of the reciprocal lattice, i.e, c i · d j = δ ij , with k i ∈ [−π, π) in 2π/N steps. With these operators the bath Hamiltonian can always be written as: whereH B (k) is a 2 × 2 matrix that can be expanded in terms of the Pauli matrices:H B (k) = (ω a + h I (k))1 + α=x,y,z h α (k)σ α . Here, the functions h α (k) contains the information about the hopping structure of the bath, and ω a is the energy of cavities describing the lattice. In what follows, we will set this energy as the energy reference of the problem, i.e., ω a ≡ 0, and replace ω 0 → ∆ = ω 0 − ω a in Eq. (1), which will now denote the detuning of the optical transition of the emitters with respect to the bath energy modes. This model can be easily diagonalized as follows: yielding two bands with energies . The relation between the original operators a k , b k and the eigenoperators u k , l k is obtained through the unitary matrix that diagonalizesH B (k), such that where the angles θ k , φ k are given by: Once this diagonalization is done, typically it is more convenient to write H I of Eq. 4 in terms of the eigenoperators u k , l k : with if the j-th emitter is coupled to the A-sublattice, i.e., n j ∈ S A , whereas: From all these previous expressions, we see that the information of both the emergent band-structure, and how the emitter couples to the bath is hidden in the h α (k), which depend on the particular hopping structure of the bath. Since we want to explore both anisotropic and tilted Dirac cones we will consider the following hopping structure (see Fig. 1(b)): • We consider nearest neighbour hoppings (in green in the figure) of overall strength −J 1 , with the possibility of a weakest one with strength −β 1 J 1 for the one corresponding within the same unit cell.
• We also consider the possibility of having next-tonearest neighbour hoppings with overall strength −J 2 , with the possibility of having a weaker link, −β 2 J 2 for the hoppings.
This hopping choice leads to the following h α (k): Notice that since h z (k) = 0 → θ k = π/2, which simplifies the calculations.

III. THEORETICAL FRAMEWORK
Along this work we are interested in studying how individual and collective spontaneous emission are modified by the interaction with the photons supported by such anisotropic Dirac photonic lattices. Since the total Hamiltonian describing the global system, i.e., H = H S + H B + H I , conserves the number of excitations, the dynamics and eigenstates of the system can be characterized independently in each excitation manifold. For example, in the single excitation subspace, the following ansatz: describes both the shape of all possible eigenstates, as well as the dynamics of the system when initialized with a single excitation. In this work, we will use this ansatz for the case when a single emitter is coupled to the bath for two purposes: • To find the emergent qubit-photon bound-states of the system [74][75][76]. These can be obtained by finding the solutions of H |Ψ BS = E BS |Ψ BS , with E BS = ω u/l (k). Using this ansatz one can find that the energy of these bound states is given by: where Σ A/B e (z) is self-energy of the emitter, which for bipartite bath reads: Here, the superindex denotes the sublattice to which the emitter is coupled to. However, since the bath has chiral symmetry, h z (k) = 0, the single emitter self-energy is the same no matter to which sublattice the emitter is coupled to, i.e., Σ A/B e (z) ≡ Σ e (z). An analytical expression for the wavefunction coefficients C O,p can also be obtained. For example, if the emitter couples to the A sublattice, the photonic component of the bound state reads (up to a normalization coefficient): where the super-index denotes the sublattice that the emitter is coupled to, and the sub-index the mode to which the spatial coefficient belongs in Eq. (15). The spatial distribution of these wavefunctions is very relevant because it eventually dictates the shape of the emitter-emitter interaction when many emitters couple to the bath. In particular, it can be shown that when the emitters' energy lie within a vanishing density of states region and the photons can be adiabatically eliminated (assuming Born-Markov conditions are satisfied [70]), the effective emitters' dynamics induced by the photons is given by [77][78][79]: where G αβ ij , which represents the interaction between the i-th emitter coupled to the α-sublattice and the j-th emitter coupled to the β-sublattice, is given (up to a constant) by G Thus, by studying already the bound state wavefunctions of a single emitter, one can already characterize the emergent coherent dipole interactions appearing when many emitters couple to the photonic bath.
• We will also use the ansatz of Eq. (15) to calculate the dynamics of the system after initializing the system in a given state, |Ψ(0) . For that one must solve the time-dependent Schrödinger equation: of the global system either by numerical or analytical means. For example, in the case of a single emitter initially excited, |Ψ(0) = σ eg |vac , the excited state population C e (t) can be expressed analytically as the following Laplace transform [80]: depending on whether the emitter couples to the α = A, B sublattice. There, we observe how the single emitter self-energy not only determines the emergence of qubit-photon bound states, but also ultimately determines its dynamics. In fact, it is possible to recover the expected Markovian result by approximating: Σ α e (E + i0 + ) = Σ α e (∆ + i0 + ) = δω α e (∆) − iΓ α e (∆)/2 in the integrand of Eq. (22), defining as δω α e , Γ α e the renormalization of the emitter's energy and lifetimes due to the interaction with the bath. Within that approximation, the excited state population will decay exponentially to the ground state: with Γ e (∆) being exactly the one predicted by Fermi's Golden Rule: that is, proportional to g 2 and the density of states of the bath, D(E) = k,p=u/l δ(E − ω p (k)). As the emitter decays to the ground state, it launches propagating excitations into the bath which can be characterized by studying C p,n (t). These dynamical spatial distributions can be also linked to the effective emitter's interactions induced by the photons. When the photons eventually propagate away from the emitter, as it is the case if their energies lie within the photonic band structure, the photoninduced interactions will eventually generate (individual and) collective dissipative terms, Γ ij , that can lead to strong super/subradiant effects when many emitters couple to the bath [64]. The spatial shape of these terms, Γ ij , is ultimately the same as the one of the propagating fields into the bath. This is why, again, by looking into the spontaneous emission of a single emitter, we will be able to extract information about the collective dissipative terms when many emitters couple to the bath.

IV. EMITTER DYNAMICS AND INTERACTIONS IN ISOTROPIC GRAPHENE
Before digging into more complex energy dispersions, let us summarize the main findings obtained for the situation of isotropic Dirac-cone dispersions [49,50]. These are obtained by setting J 2 = 0, β 1 = 1 in our model. This implies that h I (k) = 0, such that the spectrum is symmetric around the reference energy, i.e., ω u/k (k) = ±ω iso (k), with: Note that for these coordinates, (k 1 , k 2 ), the two inequivalent Dirac points appear at K, K = (2π/3)(±1, ∓1) (see Fig. 2(a)), and the energy dispersion around it is not isotropic ω iso (K ( ) + q) ≈ J 1 q 2 1 + q 2 2 − q 1 q 2 . However, the isotropy can be immediately recovered by changing to the real space coordinates (k x = (k 1 + k 2 )/3, k y = (k 1 − k 2 )/ √ 3) obtained by using that the reciprocal vectors in terms of real space axis read d 1/2 =ê x /3 ±ê y / √ 3. In these coordinates, there are six Dirac points at K ( ) (only two inequivalent ones) around which the energy dispersion is isotropic: Focusing on the regime where the energy of the emitters is tuned to the Dirac point, i.e., ∆ = 0, and using the theoretical framework explained in section III, it was found that [49]: • Since Σ α e (0) = 0, then Eq. (16) has always a trivial solution at E BS = 0. When looking at the spatial distribution of this bound-state solution, one can find that C α α,p = 0, while: where I K ( ) (p) can be shown to scale in the limit of large distances as: [49]). This means that when the emitter couples to the A [B] sublattice, a localized mode emerges around it in the opposite sublattice with an isotropic power-law 1/r-decay, plus some oscillatory factor that accounts for the interference of two integrals I K,K in Eq. (26), which does not modify the asymptotic scaling of the |C A/B a,b,p |. This very slow decay makes this bound-state very different from the standard bound-states found in other photonic band-gaps [74][75][76][77][78], since its wavefunction is not integrable in the thermodynamic limit N → ∞ (it diverges logarithmically). Since this bound state is similar to the one found in the context of graphene electronic transport [15,51,52] we will use the same nomenclature and label it as quasi-bound state, to distinguish it from conventional qubit-photon bound states [74][75][76][77][78]. In what follows, we will see how this quasi-bound state will also lead to distinctive qualitative features of the individual and collective quantum emitter dynamics.
• The excited state probability amplitude of an initially single excited emitter coupled to the bath is approximately given by: where R 0 is the overlap of the initial state with the (quasi)-bound state, and C NM (t) is a nonexponential contribution (∝ 1/ log(t)) appearing due to the non-analytical nature of the density of states around the Dirac point, i.e., D(E) ∝ |E|.
In standard photonic band-gaps [74][75][76][77][78], C NM (t) is typically hidden by the constant contribution R 0 which is generally very close to 1 and independent of system size. However, the non-integrable nature of the bound states appearing in Dirac photonic environments leads to a dramatically different behaviour. In particular, here R 0 depends on system size as follows [49]: with g(N ) = . This means that in the thermodynamic limit R 0 (N → ∞) → 0. Thus, in an infinite-size bath, one should observe only a purely non-Markovian relaxation C e (t) ≈ C NM (t), instead of the no decay expected from the vanishing density of states at the emitter's frequency, i.e., D(∆ = 0) = 0.
• Besides, when many emitters couple at these conditions to the bath, their associated quasi-bound state can mediate coherent long-range interactions, G αβ ij ∝ 1/|r ij |, as the ones predicted by Eq. (20). However, it was also shown that the finite overlap with the quasi-bound state, R 0 (N ), ultimately renormalizes these interactions as follows: Thus, for an infinite bath G αβ ij ≡ 0. Fortunately, since R 0 (N ) goes to zero logarithmically with system size, i.e., R 0 (N ) ∝ 1/ ln(N ), the quasi-bound state can effectively mediate interactions still in very large systems.
In what follows, we will consider how these quasibound states, and their quantum optical consequences, get affected by the inclusion of anisotropies and tilting in the Dirac-energy dispersion

V. PHOTON-MEDIATED INTERACTIONS IN ANISOTROPIC AND SEMI-DIRAC POINTS
Let us first consider what happens when we include an anisotropy in the nearest-neighbour hoppings, e.g., by taking β 1 = 1 in the hopping model of Fig. 1(b), and J 2 = 0. Note that here we are already choosing a particular direction of the anisotropy. However, the results can be extrapolated to the other choices just by considering a rotation of the real space coordinates. In that case, it is still satisfied that h I (k) = h z (k) = 0, such that one can use all the general expressions of the self-energies, bound-state wavefunctions, etc, derived for the isotropic case, but replacing ω iso (k) by the following anisotropic energy dispersion: The main consequence of the anisotropy β 1 is that it provides a knob to move the position of the Dirac points in the Brillouin zone [4-12, 14, 16, 17], K/K = (± arccos(−β 1 /2), ∓ arccos(−β 1 /2)), without altering its linear dispersion. This can be done until a critical anisotropy value, β 1 = 2, where the Dirac points merge forming a new type of band-crossing labeled as semi-Dirac point, linear in one direction, and parabolic in the other (see Fig. 2(b)). This can be better understood making a rotation of the axis q 1/2 = (k 2 ± k 1 )/ √ 2 and expanding the H B (k) around the K, K to obtain: for |q i | 1. Here, σ x,y are Pauli matrices, and v i the effective velocities appearing in the different q i -directions: Like this, it is trivial to see that the energy dispersion maintains its linear dispersion around the Dirac point where v 2 vanishes, such that the linear expansion of Eq. (32) is not valid anymore.
This transition from a linear to a semi-Dirac band touching translates into important qualitative differences in magnitudes that govern the quantum optical phenomena of the system, such as the bath density of states. This is what we show in Fig. 2(c) where we plot D(E) for increasing values of anisotropy ranging from β 1 = 1 (isotropic case) to β 1 = 2.1, including the semi-Dirac case, β 1 = 2 highlighted in red. There, we observe several important features: • For β 1 < 2, the density of states around the Dirac energy maintains both its singular nature (D(0) = 0) and linear energy dependence, D(E) ∝ |E|. In fact, using Eqs. (32)-(33) one can show that: for |E| J 1 . This expression is the one we plot in dashed blue in Fig. 2(c), finding a very good agreement with the numerically obtained D(E) for this energy range.
• Another important difference in the regime 1 < β 1 ≤ 2 is the appearance of two additional Van-Hove singularities between the main ones of the isotropic case and the Dirac energy. These Van-Hove singularities move closer to the Dirac energy as β 1 → 2, until they collapse at the critical value where the two Dirac energies merge, i.e., β 1 = 2. At this semi-Dirac point it can be shown that D(E) loses its linear energy dependence and transforms to D(E) ∝ |E|, as highlighted in red in Fig. 3(c), in accordance to other works [4,8,16].
• For β 1 > 2, the bands do not touch anymore, opening thus a band-gap around the Dirac energy, as clearly shown by D(E).  cut for fixed β1 = 2, and increasing system size N as depicted in the legend. Solid lines are numerical fits to a power-law decay A/n γ f . In the inset, we plot the exponent γ f obtained from the numerical fit of the asymptotic shape of the quasi-bound state. Dashed black line is set at 0.5, which is the value where γ f converges for large system sizes.
We expect that this transition from a gapless Dirac point to a gapped one passing through the semi-Dirac situation, also translates in important changes of the quasi-bound state appearing for the isotropic case. Since Σ α e (0) = 0 for all β 1 , it is guaranteed that a bound state exists at E BS = 0 (see Eq. (16)) for the whole anisotropy range. We can then characterize numerically its spatial shape by solving Eqs. (18)- (19) for several β 1 and a fixed system size. The result of this analysis is shown in Fig. 3(a-f) where we plot the bound state shape C A b,n for several β 1 ranging from β 1 = 1 (isotropic case, in panel (a)) to β 1 = 2.1 (gapped situation, in panel (b)), passing through the semi-Dirac situation β 1 = 2 (in panel (e)). From this analysis we can extract several conclusions: • As it occurs in the case of 3D photonic Weyl points [58,59], since the density of states maintains its singular nature for β 1 ∈ [1, 2] one can modify the shape of the bound state without altering its powerlaw dependence. The main change of the spatial shape is due to the modification of the interference factor of I K,K due to the movement of the Dirac points K ( ) with the anisotropy parameter β 1 . In fact, the asymptotic decay for large distances still shows a similar 1/r dependence: with a modification compared to the isotropic case of Eq. (27) provided by the effective velocities v i in the different directions. This perseverance of the 1/r decay in this regime can be better observed in Fig. 4(a) where we plot a line cut of the quasibound-state |C A b,(−n,−n) | for several β 1 's. There, we observe how for the β 1 's in that regime, i.e., β 1 = 1, 1.3, 1.6, 1.9, the asymptotic scaling of the quasibound state follows perfectly a 1/n-law, as plotted in dashed blue line. Something similar can be found along other directions (up to the oscillatory factors introduced by the interference of the two integrals of Eq. (26)).
• This situation changes dramatically at the semi-Dirac point β 2 = 2, see Fig. 3(e), where the interference between I K,K has more striking consequences. On the one hand, it leads to a directional bound state, that is, the quasi-bound state localizes mostly along one particular direction, that is, the (−n, −n) direction where the v 2 velocity vanishes. Evidently, a different choice of the anisotropy direction in the hopping model of Fig. 1(b) will result in a different localization direction, but will not alter the rest of the conclusions. Besides this directionality, the localization of the bound state can be shown to be longer ranged. This is more clear when we look at a line cut along this direction, as shown in Fig. 4(a), where we observe how the bound state decay with β 1 = 2 (yellow triangles) is significantly longer ranged than the other ones. In fact, one can make a numerical fit to a power-law decay ∝ 1/n γ f obtaining a value γ f ≈ 0.45 for the system size chosen in Fig. 4(a), i.e., N = 1000 2 . Through a numerical analysis we demonstrate that this exponent γ f depends on the system size. This is shown in Fig. 4(b) where we plot the quasi-bound state cut for β 1 = 2 and several system sizes (markers), together with its corresponding fittings (solid lines). The different slope in the logarithmic scale indicates a different exponent of the power-law decay. To make it more evident, we complement this figure with the inset, where we plot the exponent of the power-law, γ f , obtained from the numerical fit. There, we observe how γ f ranges from 0.3 for the smallest system size considered, i.e., N = 200 2 , to 0.5, which is the value it converges to in the limit of very large systems. We want to note that both the directional behaviour and 1/ √ x dependence resembles the results related to electronic and photonic transport found in the literature [15,18].
• Finally, when β 1 > 2 the quasi bound-state transforms into a standard one recovering its exponential decay dependence, as clearly shown in purple rombus in Fig. 4(a). Remarkably, though, this exponential localization preserves the one-dimensional directional behaviour from the semi-Dirac dispersion, as shown clearly in Fig. 3(f).
After this exploration of the modification of the (quasi)-bound state spatial shape, a natural question arises: will the dynamical features be similar to the ones of the isotropic model explained in the previous section, or are there important qualitative differences? This is explored in Fig. 5 by plotting the numerically calculated excited state population, |C e (t)| 2 of a single initially excited emitter coupled to a bath with strength g = 0.1J 1 , energy tuned to the Dirac point, ∆ = 0, and for increasing system sizes, N , depicted in the different colors. For comparison, we plot in the different panels (a-f) the dynamics for the same β 1 parameters chosen in Fig. 3. Besides, we perform the simulation using a bath with periodic boundary conditions to compare with the results of Ref. [49]. On top of the dynamics we plot in dashed lines the expected steady-state value, |R 0 | 2 , obtained from calculating the overlap with the quasi-bound state given by Eq. (29) but replacing ω iso (k) → ω 1 (k). From this analysis, we can extract several conclusions: • For β 1 ∈ [1, 2], and for times in which the photonic excitations do not see the border, i.e., tJ 1 /N < 1, all β 1 show a similar non-exponential decay dynamics which resembles the one of the isotropic case. Thus, for an infinite system, N → ∞, we expect (e) (f) Figure 5. (a-f) Excited state population, |Ce(t)| 2 , for an initially excited emitter for the same β1 parameters of Fig. 3. In each panel, we plot the dynamics for increasing system sizes in different colors, as depicted in the legend. The bath is simulated with periodic boundary conditions and the emitter-bath coupling strength is g = 0.1J1. In dashed lines we plot the expected steady-state values, |R0| 2 , obtained from the overlap with the quasi-bound state given by Eq. (29).
then the dynamics to be dominated purely by the non-Markovian dynamics, C NM (t) in Eq. (28), induced by the non-analytical behaviour of the density of states around the Dirac energy.
• For β 1 ∈ (1, 2) and long-enough times, i.e., tJ 1 /N > 1, the decay quenches and the population starts to oscillate around a steady-state value, like in the isotropic case. This value is determined by the overlap of the quasi-bound state with the emitter, |R 0 | 2 . As expected, R 0 depends now not only on system size, N , but also on the anisotropy parameter β 1 , i.e., R 0 ≡ R 0 (N, β 1 ). However, as we show in Figs. 5(b-e), an important difference arises with respect to the isotropic case: for fixed β 1 , R 0 does not monotonically decrease to 0 with N , but rather with some oscillations that depend on the β 1 parameter. The same occurs with the dependence with β 1 for fixed system size. As a general rule, we find that R 0 tends to zero as β 1 gets closer to the critical value, β 1 → 2, as indicated by the larger oscillations found in Figs. 5(c-d).
• At the critical point, β 1 = 2, we find that R 0 (β 1 = 2, N ) ≡ 0 for all evaluated system sizes. The underlying reason is that for this case a zero-energy extended mode appears, for the boundary conditions chosen, that resonantly couples the emitter. Thus, the emitter excitation gets resonantly transferred to such extended mode, as shown in Fig. 5(e), and it ultimately comes back coherently to the emitter at sufficiently long times (not shown).
• Finally, for β 1 > 2 one recovers the dynamical features of the standard qubit-photon bound states appearing in conventional band-gaps [81]: that is, the excitation decays initially up to a constant value |R 0 | 2 ≈ 1, and independent of system size. This is very clear in Fig. 5(e) where we see the collapse of all the lines for the dynamics of different system sizes.
Summing up, we have shown that the anisotropy of Dirac lattices provides a very useful knob to tune the spatial shape of the (quasi)-bound state, and thus, on the effective coherent interactions, G αβ ij appearing in Eq. (20). For the regime β 1 ∈ (1, 2) one can tune the spatial shape without losing its power-law behaviour, something only found so far for three-dimensional topological baths [58,59]. For β 1 > 2 one can obtain a directional and longer-ranged interaction (1/ √ r), corrected by an exponential decay length which decreases with the value of β 1 . At the critical point β 1 = 2, however, the overlap with the quasi-bound state goes to zero, such that it can not eventually mediate the interactions G αβ ij for very long times.

VI. PHOTON-MEDIATED INTERACTIONS IN TILTED DIRAC POINTS
For completeness in the study of the possible Diraclike dispersions that can appear, in this section we will consider a generalization of the hopping model studied in the previous section by including anisotropic next-nearest neighbour hoppings, i.e., making J 2 = 0 and β 2 = 1 (see Fig. 1(b)). The reason is that this simplified model will allow us to tilt the energy dispersion [19][20][21][22], and explore the quantum optical phenomena induced by the so-called Type II and III Dirac points [25][26][27][28][29][30][31][32][33][34][35][36], which occur when a critical tilting leads to a linear band-touching along the Brillouin zone, instead of a single point as in the standard (Type I) Dirac points.
In the 2 × 2 version of the bath Hamiltonian,H B (k), the addition of next-nearest neighbour hoppings enters in the diagonal part h I (k) = 0 written in Eq. (11). Thus, the energy dispersion will be given by two bands with energies: Note that h I (k) provides a k-dependent shift common to the upper and lower bands. Thus, the appearance (or not) of band-touching points is still dominated by the ω 1 (k) studied in the previous section. The main effect of h I (k) is to introduce both a tilting and a shift of the energy dispersion around the Dirac points. This is more clear when expanding the Hamiltonian around the Dirac points: for |q i | 1, whereẼ D is the shifted energy of the Dirac points given by: v 1/2 the slopes around the Dirac energy defined in Eq. (33), and v 0,1/2 the new effective diagonal velocities, which can be shown to be given by: As explained in Refs. [19,33], with these velocities, v 0,i , v i , one can define a parameter,ṽ: that characterizes the degree of tilting of the lattice, and allows to classify the different type of Dirac band touching that occur (see, e.g., Ref. [33] for an explanation of the classification). Forṽ ∈ (0, 1) the energy dispersion is tilted but still has a single-point band-touchings (Type I Dirac points). At the critical valueṽ = 1, the energy dispersion undergoes a transition from a single-point band touching to a combination of linear band touching and flat band (Type III). Finally, whenṽ > 1 the band is over-tilted showing the so-called Type II Dirac points characterized by nodal line crossings. In Figs. 6(a-d) we show an example of this transition for the energy dispersions ω 2,u/l (k) of Eq. (36). For this plot, we fix J 2 = 0.1J 1 , β 1 = 1, and move the anisotropy parameter for positive β 2 ranging from β 2 = 3 to 7.5 in 1.5 steps. From Eq. (41) and with this choice of parameters, it can be shown that the critical tiltingṽ = 1 should occur for β 2 = 6 (corresponding to Fig. 6(c)). However, what happens in this anisotropic model is more intricate. For small enough anisotropy β 2 the energy dispersions are indeed tilted and only touch at a single point corresponding to the Dirac energyẼ D , see Fig. 6(a). For bigger β 2 , though still smaller than the critical value (a) (b) (c) (d) for whichṽ = 1, the energy dispersion still have linear band touchings at the K ( ) points, however, now there are additional k-modes from the upper band that become resonant with the Dirac energy, as shown in the corners of Fig. 6(b). The latter will be important when we couple emitters since a point dipole generally couples to all k-modes which are resonant to the emitter's optical transition. At the critical value, β 2 = 6 for this choice of parameters, the energy dispersion indeed undergoes a transition to a linear band-touching around K ( ) points, although again accompanied by other resonant modes at the same energy (see Fig. 6(c)). Finally, for larger β 2 , the energy bands are overtilted and there are resonant modes at the Dirac energy for both the upper/lower bands, as it is expected for type II Dirac points.
Before considering what happens when an emitter is coupled to this bath, it is instructive to understand first how this transition is manifested in the density of states of the bath, D(E), that we know plays a very relevant role in the emitter's dynamics. In Fig. 7(a) we study the evolution of D(E) as a function of β 2 for the same parameters than Fig. 6. We take a range of β 2 ∈ (3, 7.5) in steps of 0.5, which includes the transition through the critical value β 2 = 6, highlighted in red in Fig. 7(a). There, we observe several features: i) As expected from Eq. (39), the Dirac point shifts to larger values as β 2 increases. We highlight this movement in dashed blue lines. There, we observe that for β 2 4, the density of states at the Dirac energy is approximately zero, as expected from type I Dirac points. However, for β 2 4 the density of states starts to acquire a non-zero value. This is more clear at the inset panel of Fig. 7(a), where we plot the value of the density of states precisely at the Dirac energy, D(Ẽ D ). The reason for this finite value at the density of states are precisely the k-resonant modes that start to appear at other areas of the Brillouin zone far from the linear band-touching at the K ( ) -points. As β 2 goes to a critical point, one of the Van-Hove singular-ities of the density of states approaches and crosses the Dirac point, which is one of the signatures of Type III-Dirac points. For β 2 > 6, the density of states around the Dirac point is also finite, as expected for Type II Dirac points, and increasing since more k-modes are available for the Dirac energy E D .
Let us now see what happens when emitters couple to these type of photonic baths, and how the transition through the critical Type III Dirac point affects quantum optical phenomena. As for semi-Dirac points, we will use spontaneous emission of a single initially excited emitter as it will allow to unravel phenomena which has to do both with the localization of photon emission, as it is expected in Type I Dirac points, and with the emission of propagating photons, as it is expected for Type II-III Dirac points. We start by plotting in Fig. 7(b) the excited state population of an emitter coupled, with strength g = 0.1J 1 and energy ∆ =Ẽ D , to the central Asite of a photonic lattice with the parameters of Fig. 6(ad). There, we observe how the slow non-Markovian relaxation expected for Type I Dirac points is only observed for the smallest β 2 chosen, where there are no additional k-resonant modes. For the other cases, the emitter dynamics displays the standard exponential relaxation to the ground state, with a decay rate proportional to the value of the density of states at the emitter's energy, D(Ẽ D ). Thus, the presence of other resonant kmodes does not allow to observe any strong signature of the crossing of the critical point in the quantum emitter dynamics.
On the contrary, one can indeed observe some signatures by looking the radiation pattern in the bath, given by the wavefunctions C A a/b,n (t) of Eq. (15). This is what we plot in Fig. 8 by choosing the same parameters of Fig. 7(b) and looking at the wavefunction at the time t = 0.3 √ N /J 1 . We only plot the population in the bmodes since the one of the a sites displays a qualitatively similar behaviour (except for β 2 = 3 where the a  Fig. 6. The emitter is always tuned to the Dirac energy positionẼD so that it always probes the Dirac energy spectrum modes shows no population). In these panels, we can observe how for β 2 4 (Fig. 8a) the light becomes strongly localized around the emitter, being thus dominated by the quasi-bound state physics discussed in section V. When the other k-modes become resonant, like in β 2 = 4.5 (Fig. 8b), the population into the bath is a superposition between the quasi-bound state localization and the propagating emission coming from the other modes available. At the critical point (Fig. 8c), the emis- sion becomes strongly directional due to the presence of a Van-Hove singularity around the Dirac point [61][62][63]. As we will see next, this directional emission is a strong indication of the crossing of the critical point of Type III Dirac points. Here, we want to note that one also gets emission from the other k-resonant modes, however, since the density of states associated to the Van-Hove singularity is typically much larger, it eventually dominates the emission. Finally, for β 2 = 7.5 corresponding to the Type-II situation (see Fig. 8d), the emission starts to lose its directional character due to the increasing number of k-modes that become resonant with the emitter's energy (see Fig. 6d).
As we have just seen, for the simple model we have chosen the type I-III-II transition is hindered by the presence of additional k-modes at the emitter's energy. A natural question then arises: how will spontaneous emission look like if the emitter could couple perfectly only to k modes where the linear band-touching appear? One option to obtain this would be to keep the same photonic lattice model, but get of rid of the coupling to these spurious modes. This could be done, for example, using giant atoms [82][83][84][85][86] and designing their couplings such that the emitter only couple to the linear band-touching modes.
Here, however, we will follow a cleaner approach just for the sake of illustration. We will take an alternative photonic lattice model proposed in Ref. [34], especially designed to have a perfect Type-III energy dispersion. The model can be written as a bipartite bath as follows: where d(k) = 2J 1 (cos k 1 +cos k 2 ) and a(k) = 2J 2 (cos k 1 − cos k 2 ). This model has two bands ω u/l (k) = a(k) ± a(k) 2 + d(k) 2 . Since d(k) vanishes at k 1 ±(∓)k 2 = ±π, and a(k) vanishes when k 1 = ±k 2 , the type III Dirac cones appear at the intersection of those lines, i.e, at the points k = (±π/2, ±π/2), as can be seen in Fig. 9(a) where we plotted ω u/l (k) for J 2 = 0.3J 1 .
Apart from having a perfect singular band-touching, the energy dispersion along the straight lines k 1 ±(∓)k 2 = ±π are perfectly flat, which ultimately results in a highly divergent density of states, i.e., Van-Hove singularity, at the Dirac energy. This is clearly seen in Fig. 9(b), where we plot D(E) for this model with J 2 = 0.3J 1 . This divergent density of states has strong impact in the quantum emitter dynamics, that we plot in Fig. 9(c), for an emitter coupled to the central A site of a photonic lattice with J 2 = 0.3J 1 and N = 1024 2 sites. Differently from the critical point dynamics of Fig. 7(b), the emitter now displays a strong non-Markovian relaxation showing even some reversible oscillations. Besides, these oscillations are also accompanied by a highly directional emission into the bath, as plotted in Fig. 9(d). Thus, both the reversible dynamics and directional emission can be considered as smoking guns of the critical Dirac transition.

VII. ANISOTROPIC AND TILTED ENERGY DISPERSIONS IN SUBWAVELENGTH ATOMIC ARRAYS
In this section, we discuss a possible implementation to observe the phenomena predicted along this manuscript based on the recent proposals for exploring quantum optical phenomena using subwavelength atomic arrays [67][68][69]. The idea of these proposals is sketched in Fig. 10(a): a few impurity atoms (in blue) are placed near a subwavelength atomic array (in green), so that they can couple to the subradiant atomic excitations that propagate within the array with an energy dispersion, ω(k), which depends on the geometry/atomic parameters [87][88][89][90][91][92][93][94][95][96][97][98]. Such subwavelength atomic arrays have been recently implemented using Rubidium atoms in square geometries [65] for inter-atomic distances d/λ a ≈ 0.68, with λ a being the wavelength associated to the optical transition. However, one can reach smaller distance regimes using Alkaline-Earth atoms [67][68][69], since they feature a combination of short-and long-wavelength transitions [99]. Besides, let us also note that other geometries, such as hexagonal or triangular lattices [100][101][102][103], have also been obtained using different laser configurations to generate the optical trapping potentials for the atomic arrays.
The goal of this section is to find a subwavelength atomic array configuration where the non-trivial energy dispersions discussed along this manuscript (semi-and tilted-Dirac points) appear for the subradiant modes of the array. For that purpose, we consider an hexagonal geometry as the one depicted in Fig. 10(a): it is described by a two-site Bravais lattice with primitive vectors c 1/2 = d √ 3/2( √ 3ê x ±ê y ), and basis vector d bas = −d intraêx . For standard honeycomb lattices like the one depicted in the figure, the distance between the nearest-neighbours in all directions are equal, i.e., d intra ≡ d. However, in order to realise honeycomb lattices with uniaxial anisotropy we let d intra = d, while keeping the distance between the atoms within the same sublattice ( √ 3d) fixed. We characterize the degree of anisotropy by a parameter defined as β = d intra /d inter , where d inter is the intercell nearest neighbour distance as indicated in Fig. 10(a). Hence, β > (<)1 represents an anisotropic lattice where the two atoms in the unit cell are pushed together (apart), respectively. In what follows, we will show how these anisotropic honeycomb atomic arrays display dispersion relations with semi-Dirac and tilted Dirac cones.
We initially consider a 2D atomic array arranged in a honeycomb lattice β = 1. Each atom is assumed to be a two level system with resonance frequency ω a = 2πc/λ a and polarization dipole ℘ a . The dynamics of this atomic array can be described within an stochastic wavefunction approach through the following effective non-Hermitian (f) (g) Figure 10. (a) Subwavelength atomic array with tunable honeycomb geometry hosts subradiant atomic excitations with nontrivial Dirac energy dispersions ω(k). When additional atoms are placed nearby the structure, they can exchange interactions with these subradiant modes, mimicking light-matter couplings. We tune the honeycomb geometry keeping the primitive lattice vectors fixed: c 1/2 = d √ 3/2( √ 3êx ±êy), and modifying the ratio between the intra/inter-cell distance given by β = dintra/dinter. Like this, one generates lattices with uniaxial anisotropy and tune the energy dispersions of the subradiant excitations. (b) Dispersion relation of an isotropic honeycomb lattice, β = 1, with dinter = dintra = 0.15λa, for out of plane modes. (c,d) We keep d = 0.15λa fixed but change β < 1, and plot the dispersion of the out of plane modes for this anisotropic lattice. A semi-Dirac point emerges at M for β = 0.85, as shown by plotting the two paths in the Brillouin zone sketched in the inset: while the dispersion is linear along kx, it is quadratic along ky. (e) Dispersion relation of the in plane modes of the same isotropic honeycomb lattice as in (b). (f,g) We now tune β = 1.15 > 1, and show how a tilted Dirac cone emerges for the in-plane modes of this anisotropic lattice. We plot the dispersion along kx (f) and ky (g) following the two different paths in the Brillouin zone sketched in the inset. In all panels the color scale represents the decay rate, γ k , of the modes as indicated in the legend.
Hamiltonian [87,96]: where j is a index running over all atoms in the metasurface (N A ), placed at positions r j , and Γ a = |℘ a | 2 ω 3 a /(3π c 3 ) is the individual free-space decay rate. The coherent (J ij ) and incoherent (Γ ij ) photon-mediated interactions among emitters are given by the free space Green dyadic G 0 (r i − r j ) [104,105]: where℘ i = ℘ i /|℘ i |, and: Note that the couplings within the atomic array have two important differences with respect to the simplified coupled resonator models we consider along the manuscript: i) the photon-exchange interactions within the atomic array are long-ranged; ii) they depend on the polarization of the optical transitions. Thus, finding the same energy energy dispersions than in the nearest neighbour model is not guaranteed. In order to find such energy dispersions we restrict to the single-excitation subspace and the infinite size limits, where the eigenstates of the Hamiltonian H m are Bloch functions, where k = (k x , k y ) is the Bloch wavevector in the first Brillouin zone, and we explicitly take into account that we have a non-Bravais lattice, with the sums now running over n, up to the total number of unit cells (N ) and over the two sites per unit cell (m). Obtaining the eigenenergies of the Bloch modes from the above Hamiltonian then reduces to diagonalising the following matrix [96]: Rn =0 e −ik·Rn G αβ 0 (R n )δ mµ δ mν + Rn e −ik·Rn G αβ 0 (R n + d bas )δ 1µ δ 2ν where {R n } represent the in-plane position vector of all the unit cells in the lattice. Additionally, {α, β} run over the 3 spatial degrees of freedom, and {µ, ν} over the 2 lattice sites, such that M k is a 6×6 matrix with eigenvalues ω k − i γ k 2 , from which we obtain the photonic band structure and the radiative decay of the Bloch modes. The lattice sums of the Green's tensor are performed using the Ewald method [106,107]. Figure 10(b,e) shows the band structure of a subwavelength isotropic honeycomb lattice (d inter = d intra = 0.15λ a ) with the decay rate in color code, for modes where the dipole moments lie out of the plane (b) and in the plane (e), and for the path in the Brillouin zone sketched in the inset. In both cases we see Dirac points at K and K' [95], as well as the characteristic features (polariton-like bends and divergences) when the modes interact with the light line. While modes within the light cone (gray area) couple to the free-space photonic continuum and hence present a non-zero decay rate, the Dirac cones appear outside the light cone where modes are confined to the array and thus subradiant.
Next we introduce uniaxial anisotropy in the metasurface, by allowing β = 1. By calculating the band structures, we find two critical values of the anisotropy parameters where generalised Dirac dispersions emerge. First, for β = 0.85 a semi-Dirac cone emerges for the out-of-plane modes. This is shown in panels (c,d), where we show the bands along the ΓM direction displaying a Dirac crossing with linear dispersion at M (panel c), while along the K M K direction the dispersion away from the degeneracy is quadratic (panel d). On the other hand, for β = 1.15, two of the in-plane modes intersect forming a tilted Dirac cone (f,g). This occurs along the ΓM line but away from any high symmetry points, as shown in (panel f). The dispersion in the perpendicular direction is plotted in (panel g), showing also the linear dispersion away from the degeneracy.
Overall, these calculations show that both semi-Dirac and tilted Dirac energy dispersion can emerge in such subwavelength arrays. Then, as explained in Refs. [67][68][69], in order to probe the physics explored along this manuscript one can place additional impurity atoms (as depicted in Fig. 10(a)) with their optical transitions ω 0 tuned to the spectral region where these non-trivial bandcrossings occur, e.g., using Raman-assisted processes.

VIII. CONCLUSION & OUTLOOK
Summing up, we have characterized the spontaneous emission properties of emitters spectrally tuned to several types of Dirac points using discrete photonic lattice models. For the case of anisotropic Dirac energy dispersions, we have shown how the shape of the emergent qubit-photon bound states and photon-mediated interactions can be modified without altering the asymptotic decay of the isotropic Dirac situation. This tuning can be done until the critical point (semi-Dirac case) where a directional localized quasi-bound state appears with a longer-ranged spatial decay. Besides, we have also explored the dynamical and radiation features of the case of tilted Dirac points, including a model that features a perfect critical tilting (type III Dirac point). There, we observe how the Dirac-Type I-III-II transition manifests in reversible decay dynamics and directional emission patterns. Finally, we demonstrate that these energy dispersions can be obtained using subwavelength atomic arrays despite the intrinsically long-range nature of their interaction when compared to the simplified models. Thus, this platform appears as a very suitable one to probe these phenomena using additional impurity atoms.
An interesting outlook of this work is to consider other types of two-dimensional band-touchings, such as purely parabolic ones that appear in bilayer systems [108], to study the interplay of these phenomena with different boundary conditions, or to harness non-local light-matter couplings [82][83][84][85][86] to be further engineer the quantum emitter dynamics.