Quantum electrodynamics in anisotropic and tilted Dirac photonic lattices

The extraordinary electronic properties of graphene are directly related to the emergence of isotropic Dirac-cone 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–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–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 types II and III Dirac points [25–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 solid-state 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 engineered to obtain unconventional photon transport properties [18, 33, 36, 41–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 non-exponential 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–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–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 when many emitters couple to the photonic bath, as the ones found in photonic Van-Hove singularities [61–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–69].

The manuscript is structured as follows: in sections 2 and 3 we describe the system under study and the theoretical framework that we will use to characterize it, respectively. In section 4, we review the results known for isotropic Dirac dispersions, and then in sections 5 and 6 we analyze the change of the photon-mediated interactions with the degree of anisotropy and tilting of the photonic bath energy dispersion. Then, in section 7, 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 8.

The system that we consider along this manuscript (see scheme of figure 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 ℏω₀. The corresponding Hamiltonian then reads (assuming ℏ ≡ 1 for the rest of the manuscript):

$$\begin{equation}{H}_{\text{S}}={\omega }_{0}\sum\limits _{j}{\sigma }_{ee}^{j},\end{equation} \tag{ 1 }$$

where we use the notation ${\sigma }_{\alpha \beta }^{j}=\vert {\alpha \rangle }_{j}\langle \beta {\vert }_{j}$ for the operator acting on the jth emitter. Through their optical transition, these emitters couple to a photonic bath site, described by a bosonic annihilation (creation) operator ${O}_{\mathbf{n}}^{({\dagger})}$, through the standard light–matter Hamiltonian:

$$\begin{equation}{H}_{\text{I}}=\sum\limits _{j}\left({g}_{j}{O}_{{\mathbf{n}}_{j}}{\sigma }_{eg}^{j}+\mathrm{h}.\mathrm{c}.\right),\end{equation} \tag{ 2 }$$

where g_j denotes the strength of the coupling for the jth emitter, and n_j denotes the position of the photonic bath where the emitter couples to. For simplicity we will assume this coupling to be emitter-independent, g_j ≡ g. Note also that in equation (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}_{{\mathbf{n}}_{j}}{\sigma }_{ge}^{j}$,...). 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–73], that is one of the most promising implementations to observe these effects.

Zoom In Zoom Out Reset image size

The non-trivial part of the problem stems from the photonic bath. Instead of taking a continuum description for the photonic bath, we will use a discrete photonic lattice model described as an array of coupled cavities. Like this, it will be easier to extract analytical and numerical insight of the physics. Since we want to describe Dirac-like dispersions the array is disposed in a honeycomb geometry as depicted in figure 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₁, n₂) with ${n}_{i}\in \mathbb{Z}$ in terms of the primitive lattice vector displacements, i.e. $\mathbf{n}\to {\sum }_{i=1}^{2}{n}_{i}{\mathbf{c}}_{i}$, with c_i being the primitive vectors of the lattice ${\mathbf{c}}_{1/2}=3/2{\hat{e}}_{x}\pm \sqrt{3}/2{\hat{e}}_{y}$ (taking the lattice constant as the unit of length). Imposing periodic boundary conditions, we can define ${O}_{\mathbf{k}}=1/\sqrt{N}{\sum }_{\mathbf{n}}{\text{e}}^{-\text{i}\mathbf{k}\cdot \mathbf{n}}{O}_{\mathbf{n}}$, for O = a, b, being N the total number of sites of each sublattice, and $\mathbf{k}\to {\sum }_{i=1}^{2}{k}_{i}{\mathbf{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 photonic bath Hamiltonian can always be written as:

$$\begin{equation}{H}_{\text{B}}=\sum\limits _{\mathbf{k}}({a}_{\mathbf{k}}^{{\dagger}}\enspace {b}_{\mathbf{k}}^{{\dagger}}){\tilde{H}}_{\text{B}}(\mathbf{k})\left(\begin{matrix}\hfill {a}_{\mathbf{k}}\hfill \\ \hfill {b}_{\mathbf{k}}\hfill \end{matrix}\right),\end{equation} \tag{ 3 }$$

where ${\tilde{H}}_{\text{B}}(\mathbf{k})$ is a 2 × 2 matrix that can be expanded in terms of the Pauli matrices: ${\tilde{H}}_{\text{B}}(\mathbf{k})=({\omega }_{a}+{h}_{\text{I}}(\mathbf{k}))\mathbf{1}+{\sum }_{\alpha =x,y,z}{h}_{\alpha }(\mathbf{k}){\sigma }_{\alpha }$. Here, the functions h_α (k) contains the information about the hopping structure of the photonic 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 ω₀ → Δ = ω₀ − ω_a in equation (1), which will now denote the detuning of the optical transition of the emitters with respect to the photonic bath energy modes. This model can be easily diagonalized as follows:

$$\begin{equation}{H}_{\text{B}}=\sum\limits _{\mathbf{k}}\left({\omega }_{u}(\mathbf{k}){u}_{\mathbf{k}}^{{\dagger}}{u}_{\mathbf{k}}+{\omega }_{l}(\mathbf{k}){l}_{\mathbf{k}}^{{\dagger}}{l}_{\mathbf{k}}\right),\end{equation} \tag{ 4 }$$

yielding two bands with energies ${\omega }_{u/l}(\mathbf{k})={h}_{\text{I}}(\mathbf{k})\pm \sqrt{\vert {h}_{x}(\mathbf{k}){\vert }^{2}+\vert {h}_{y}(\mathbf{k}){\vert }^{2}+\vert {h}_{z}(\mathbf{k}){\vert }^{2}}={h}_{\text{I}}(\mathbf{k})\pm \omega (\mathbf{k})$. The relation between the original operators a_k , b_k and the eigen-operators u_k , l_k is obtained through the unitary matrix that diagonalizes ${\tilde{H}}_{\text{B}}(\mathbf{k})$, such that

$$\begin{equation}\left(\begin{matrix}{a}_{\mathbf{k}}\\ {b}_{\mathbf{k}}\end{matrix}\right)={\tilde{U}}_{\mathbf{k}}\left(\begin{matrix}{u}_{\mathbf{k}}\\ {l}_{\mathbf{k}}\end{matrix}\right)=\left(\begin{matrix}\mathrm{cos}\left(\frac{{\theta }_{\mathbf{k}}}{2}\right){\text{e}}^{-\text{i}{\phi }_{\mathbf{k}}}& \mathrm{sin}\left(\frac{{\theta }_{\mathbf{k}}}{2}\right){\text{e}}^{-\text{i}{\phi }_{\mathbf{k}}}\\ \mathrm{sin}\left(\frac{{\theta }_{\mathbf{k}}}{2}\right)& -\mathrm{cos}\left(\frac{{\theta }_{\mathbf{k}}}{2}\right)\end{matrix}\right)\left(\begin{matrix}{u}_{\mathbf{k}}\\ {l}_{\mathbf{k}}\end{matrix}\right),\end{equation} \tag{ 5 }$$

where the angles θ_k , ϕ_k are given by:

$$\begin{equation}\mathrm{cos}\left({\theta }_{\mathbf{k}}\right)=\frac{{h}_{z}(\mathbf{k})}{\omega (\mathbf{k})},\end{equation} \tag{ 6 }$$

$$\begin{equation}\mathrm{tan}\left({\phi }_{\mathbf{k}}\right)=\frac{{h}_{y}(\mathbf{k})}{{h}_{x}(\mathbf{k})}.\end{equation} \tag{ 7 }$$

Once this diagonalization is done, typically it is more convenient to write H_I of equation (4) in terms of the eigen-operators u_k , l_k :

$$\begin{equation}{H}_{\text{I}}=\sum\limits _{j}\left({B}_{{\mathbf{n}}_{j}}{\sigma }_{eg}^{j}+\mathrm{h}.\mathrm{c}.\right),\end{equation} \tag{ 8 }$$

with

$$\begin{equation}{B}_{{\mathbf{n}}_{j}}=\frac{g}{\sqrt{N}}\sum\limits _{\mathbf{k}}{\text{e}}^{\text{i}(\mathbf{k}\cdot {\mathbf{n}}_{j}-{\phi }_{\mathbf{k}})}\left[\mathrm{cos}({\theta }_{\mathbf{k}}/2){u}_{\mathbf{k}}+\mathrm{sin}({\theta }_{\mathbf{k}}/2){l}_{\mathbf{k}}\right],\end{equation} \tag{ 9 }$$

if the jth emitter is coupled to the A-sublattice, i.e. n_j ∈ S_A , whereas:

$$\begin{equation}{B}_{{\mathbf{n}}_{j}}=\frac{g}{\sqrt{N}}\sum\limits _{\mathbf{k}}{\text{e}}^{\text{i}\mathbf{k}\cdot {\mathbf{n}}_{j}}\left[\mathrm{sin}({\theta }_{\mathbf{k}}/2){u}_{\mathbf{k}}-\mathrm{cos}({\theta }_{\mathbf{k}}/2){l}_{\mathbf{k}}\right]\end{equation} \tag{ 10 }$$

if n_j ∈ S_B .

From all these previous expressions, we see that the information of both the emergent band-structure, and how the emitter couples to the photonic bath is hidden in the h_α (k), which depend on the particular hopping structure of the photonic bath. Since we want to explore both anisotropic and tilted Dirac cones we will consider the following hopping structure (see figure 1(b)):

• (a)  
  We consider nearest neighbour hoppings (in green in the figure) of overall strength −J₁, with the possibility of a weakest one with strength −β₁ J₁ for the one corresponding within the same unit cell.
• (b)  
  We also consider the possibility of having next-to-nearest neighbour hoppings with overall strength −J₂, with the possibility of having a weaker link, −β₂ J₂ for the hoppings.

This hopping choice leads to the following h_α (k):

$$\begin{equation}{h}_{\text{I}}(\mathbf{k})=-2{J}_{2}\left[\mathrm{cos}({k}_{1})+\mathrm{cos}({k}_{2})+{\beta }_{2}\enspace \mathrm{cos}({k}_{1}-{k}_{2})\right],\end{equation} \tag{ 11 }$$

$$\begin{equation}{h}_{z}(\mathbf{k})=0,\end{equation} \tag{ 12 }$$

$$\begin{equation}{h}_{x}(\mathbf{k})=-{J}_{1}\left[{\beta }_{1}+\mathrm{cos}({k}_{1})+\mathrm{cos}({k}_{2})\right],\end{equation} \tag{ 13 }$$

$$\begin{equation}{h}_{y}(\mathbf{k})=-{J}_{1}\left[\mathrm{sin}({k}_{1})+\mathrm{sin}({k}_{2})\right].\end{equation} \tag{ 14 }$$

Notice that since h_z (k) = 0 → θ_k = π/2, which simplifies the calculations.

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:

$$\begin{equation}\vert {\Psi}\rangle =\left[\sum\limits _{j}{C}_{e,j}{\sigma }_{eg}^{j}+\sum\limits _{O=a,b}\sum\limits _{\mathbf{p}}{C}_{O,\mathbf{p}}{O}_{\mathbf{p}}^{{\dagger}}\right]\vert \mathrm{v}\mathrm{a}\mathrm{c}\rangle ,\end{equation} \tag{ 15 }$$

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 photonic bath for two purposes:

• (a)  
  To find the emergent qubit-photon bound-states of the system [74–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:

  $$\begin{equation}{E}_{\text{BS}}={\Delta}+{{\Sigma}}_{e}^{\alpha }({E}_{\text{BS}}),\end{equation} \tag{ 16 }$$

  where ${{\Sigma}}_{e}^{A/B}(z)$ is self-energy of the emitter, which for bipartite photonic bath reads:

  $$\begin{equation}{{\Sigma}}_{e}^{A/B}(z)=\frac{{g}^{2}}{N}\sum\limits _{\mathbf{k}}\left[\frac{z-{h}_{\text{I}}(\mathbf{k})}{{(z-{h}_{\text{I}}(\mathbf{k}))}^{2}-\omega {(\mathbf{k})}^{2}}\right].\end{equation} \tag{ 17 }$$

  Here, the superindex denotes the sublattice to which the emitter is coupled to. However, since the photonic 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. ${{\Sigma}}_{e}^{A/B}(z)\equiv {{\Sigma}}_{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):

  $$\begin{equation}{C}_{a,\mathbf{p}}^{A}=\frac{1}{N}\sum\limits _{\mathbf{k}}\left[\frac{{E}_{\text{BS}}-{h}_{\text{I}}(\mathbf{k})}{{({E}_{\text{BS}}-{h}_{\text{I}}(\mathbf{k}))}^{2}-\omega {(\mathbf{k})}^{2}}\right]{\text{e}}^{\text{i}\mathbf{k}\cdot \mathbf{p}},\end{equation} \tag{ 18 }$$

  $$\begin{equation}{C}_{b,\mathbf{p}}^{A}=\frac{1}{N}\sum\limits _{\mathbf{k}}\frac{\omega (\mathbf{k})}{{({E}_{\text{BS}}-{h}_{\text{I}}(\mathbf{k}))}^{2}-\omega {(\mathbf{k})}^{2}}{\text{e}}^{\text{i}(\mathbf{k}\cdot \mathbf{p}-{\phi }_{\mathbf{k}})},\end{equation} \tag{ 19 }$$

  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 equation (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 photonic 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–79]:

  $$\begin{equation}{H}_{\text{eff}}\approx \sum\limits _{ij}{G}_{ji}^{\alpha \beta }{\sigma }_{eg}^{j}{\sigma }_{ge}^{i},\end{equation} \tag{ 20 }$$

  where ${G}_{ij}^{\alpha \beta }$, which represents the interaction between the ith emitter coupled to the α-sublattice and the jth emitter coupled to the β-sublattice, is given (up to a constant) by ${G}_{ij}^{AA/BB}\propto {C}_{a/b,{\mathbf{r}}_{i}-{\mathbf{r}}_{j}}^{A/B}$ and ${G}_{ij}^{AB}\propto {C}_{b,{\mathbf{r}}_{i}-{\mathbf{r}}_{j}}^{A}$. 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.
• (b)  
  We will also use the ansatz of equation (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:

  $$\begin{equation}\vert {\Psi}(t)\rangle ={\text{e}}^{-\text{i}Ht}\vert {\Psi}(0)\rangle ,\end{equation} \tag{ 21 }$$

  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]:

  $$\begin{equation}{C}_{e}^{\alpha }(t)=-\frac{1}{2\pi \text{i}}{\int }_{-\infty }^{\infty }\mathrm{d}E\frac{{\text{e}}^{-\text{i}Et}}{E+\text{i}{0}^{+}-{\Delta}-{{\Sigma}}_{e}^{\alpha }(E+\text{i}{0}^{+})},\end{equation} \tag{ 22 }$$

  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: ${{\Sigma}}_{e}^{\alpha }(E+\text{i}{0}^{+})={{\Sigma}}_{e}^{\alpha }({\Delta}+\text{i}{0}^{+})=\delta {\omega }_{e}^{\alpha }({\Delta})-\text{i}{{\Gamma}}_{e}^{\alpha }({\Delta})/2$ in the integrand of equation (22), defining as $\delta {\omega }_{e}^{\alpha },{{\Gamma}}_{e}^{\alpha }$ the renormalization of the emitter's energy and lifetimes due to the interaction with the photonic bath. Within that approximation, the excited state population will decay exponentially to the ground state:

  $$\begin{equation}\vert {C}_{e}^{\alpha }(t){\vert }^{2}\approx {\text{e}}^{-{{\Gamma}}_{e}({\Delta})t},\end{equation} \tag{ 23 }$$

  with Γ_e (Δ) being exactly the one predicted by Fermi's golden rule:

  $$\begin{equation}{{\Gamma}}_{e}^{\alpha }({\Delta})=\pi {g}^{2}\sum\limits _{\mathbf{k},p=u/l}\delta ({\Delta}-{\omega }_{p}(\mathbf{k})),\end{equation} \tag{ 24 }$$

  that is, proportional to g² and the density of states of the photonic bath, D(E) = ∑_k,p=u/l δ(E − ω_p (k)). As the emitter decays to the ground state, it launches propagating excitations into the photonic 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 photon-induced interactions will eventually generate (individual and) collective dissipative terms, Γ_ij , that can lead to strong super/subradiant effects when many emitters couple to the photonic bath [64]. The spatial shape of these terms, Γ_ij , is ultimately the same as the one of the propagating fields into the photonic 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 photonic bath.

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₂ = 0, β₁ = 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:

$$\begin{equation}{\omega }_{\text{iso}}(\mathbf{k})={J}_{1}\sqrt{3+2\left[\mathrm{cos}({k}_{1})+\mathrm{cos}({k}_{2})+\mathrm{cos}({k}_{1}-{k}_{2})\right]}.\end{equation} \tag{ 25 }$$

Note that for these coordinates, (k₁, k₂), the two inequivalent Dirac points appear at K, K' = (2π/3) (±1, ∓1) (see figure 2(a)), and the energy dispersion around it is not isotropic ${\omega }_{\text{iso}}({\mathbf{K}}^{(\prime )}+\mathbf{q})\approx {J}_{1}\sqrt{{q}_{1}^{2}+{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})/\sqrt{3})$ obtained by using that the reciprocal vectors in terms of real space axis read ${\mathbf{d}}_{1/2}={\hat{e}}_{x}/3\pm {\hat{e}}_{y}/\sqrt{3}$. In these coordinates, there are six Dirac points at K⁽'⁾ (only two inequivalent ones) around which the energy dispersion is isotropic: ${\omega }_{\text{iso}}({\mathbf{K}}^{(\prime )}+\mathbf{q})\approx \frac{3{J}_{1}}{2}\vert \mathbf{q}\vert$.

Zoom In Zoom Out Reset image size

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 3, it was found that [49]:

• (a)  
  Since ${{\Sigma}}_{e}^{\alpha }(0)=0$, then equation (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}_{\alpha ,\mathbf{p}}^{\alpha }=0$, while:

  $$\begin{equation}{C}_{b/a,\mathbf{p}}^{A/B}=\left[{\text{e}}^{\text{i}\mathbf{K}\cdot \mathbf{p}}{I}_{\mathbf{K}}(\mathbf{p})+{\text{e}}^{\text{i}{\mathbf{K}}^{\prime }\cdot \mathbf{p}}{I}_{{\mathbf{K}}^{\prime }}(\mathbf{p})\right],\end{equation} \tag{ 26 }$$

  where ${I}_{{\mathbf{K}}^{(\prime )}}(\mathbf{p})$ can be shown to scale in the limit of large distances as:

  $$\begin{equation}\vert {I}_{{\mathbf{K}}^{(\prime )}}(\mathbf{p})\vert \propto \frac{1}{\vert \tilde{\mathbf{p}}\vert },\end{equation} \tag{ 27 }$$

  with $\tilde{\mathbf{p}}=\left(3({p}_{1}+{p}_{2})/2,\sqrt{3}({p}_{1}-{p}_{2})/2\right)$ (see reference [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 equation (26), which does not modify the asymptotic scaling of the $\vert {C}_{a,b,\mathbf{p}}^{A/B}\vert$. This very slow decay makes this bound-state very different from the standard bound-states found in other photonic band-gaps [74–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–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.
• (b)  
  The excited state probability amplitude of an initially single excited emitter coupled to the photonic bath is approximately given by:

  $$\begin{equation}{C}_{e}(t)\approx {R}_{0}+{C}_{\text{NM}}(t),\end{equation} \tag{ 28 }$$

  where R₀ is the overlap of the initial state with the (quasi)-bound state, and C_NM(t) is a non-exponential 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–78], C_NM(t) is typically hidden by the constant contribution R₀ 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₀ depends on system size as follows [49]:

  $$\begin{equation}{R}_{0}(N)=\frac{1}{1+\frac{{g}^{2}}{{J}_{1}^{2}}g(N)},\end{equation} \tag{ 29 }$$

  with $g(N)=\frac{{J}_{1}^{2}}{N}{\sum }_{\mathbf{k}}\frac{1}{{\omega }_{\text{iso}}{(\mathbf{k})}^{2}}\approx 0.2+0.37\enspace \mathrm{ln}(N)$. This means that in the thermodynamic limit R₀(N → ∞) → 0. Thus, in an infinite-size photonic 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.
• (c)  
  Besides, when many emitters couple at these conditions to the photonic bath, their associated quasi-bound state can mediate coherent long-range interactions, ${G}_{ij}^{\alpha \beta }\propto 1/\vert {\mathbf{r}}_{ij}\vert$, as the ones predicted by equation (20). However, it was also shown that the finite overlap with the quasi-bound state, R₀(N), ultimately renormalizes these interactions as follows:

  $$\begin{equation}{G}_{ij}^{\alpha \beta }\propto \frac{{g}^{2}{R}_{0}(N)}{{J}_{1}\vert {\tilde{\mathbf{r}}}_{ij}\vert }.\end{equation} \tag{ 30 }$$

  Thus, for an infinite photonic bath ${G}_{ij}^{\alpha \beta }\equiv 0$. Fortunately, since R₀(N) goes to zero logarithmically with system size, i.e. R₀(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 quasi-bound states, and their quantum optical consequences, get affected by the inclusion of anisotropies and tilting in the Dirac-energy dispersion.

Let us first consider what happens when we include an anisotropy in the nearest-neighbour hoppings, e.g. by taking β₁ ≠ 1 in the hopping model of figure 1(b), and J₂ = 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:

$$\begin{equation}\frac{{\omega }_{1}(\mathbf{k})}{{J}_{1}}=\sqrt{2+{\beta }_{1}^{2}+2\enspace \mathrm{cos}({k}_{1}-{k}_{2})+2{\beta }_{1}\left[\mathrm{cos}({k}_{1})+\mathrm{cos}({k}_{2})\right]}.\end{equation} \tag{ 31 }$$

The main consequence of the anisotropy β₁ 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(−β₁/2), ∓ arccos(−β₁/2)), without altering its linear dispersion. This can be done until a critical anisotropy value, β₁ = 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 figure 2(b)). This can be better understood making a rotation of the axis ${q}_{1/2}=({k}_{2}\pm {k}_{1})/\sqrt{2}$ and expanding the H_B(k) around the K, K' to obtain:

$$\begin{equation}{\tilde{H}}_{\text{B}}(\mathbf{K}/{\mathbf{K}}^{\prime }+\mathbf{q})\approx {v}_{1}{\sigma }_{y}{q}_{1}\mp {v}_{2}{\sigma }_{x}{q}_{2},\end{equation} \tag{ 32 }$$

for |q_i | ≪ 1. Here, σ_x,y are Pauli matrices, and v_i the effective velocities appearing in the different q_i -directions:

$$\begin{equation}{v}_{2}={J}_{1}\sqrt{2-\frac{{\beta }_{1}^{2}}{2}},\qquad {v}_{1}={J}_{1}\frac{{\beta }_{1}}{\sqrt{2}}.\end{equation} \tag{ 33 }$$

Like this, it is trivial to see that the energy dispersion maintains its linear dispersion around the Dirac point ${\omega }_{1}({\mathbf{K}}^{(\prime )}+\mathbf{q})\approx \sqrt{{v}_{1}^{2}{q}_{1}^{2}+{v}_{2}^{2}{q}_{2}^{2}}$, until β₁ = 2, where v₂ vanishes, such that the linear expansion of equation (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 photonic bath density of states. This is what we show in figure 2(c) where we plot D(E) for increasing values of anisotropy ranging from β₁ = 1 (isotropic case) to β₁ = 2.1, including the semi-Dirac case, β₁ = 2 highlighted in red. There, we observe several important features:

• (a)  
  For β₁ < 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 equations (32) and (33) one can show that:

  $$\begin{equation}D(E)\simeq \frac{\vert E\vert }{2\pi {v}_{1}{v}_{2}},\end{equation} \tag{ 34 }$$

  for |E| ≪ J₁. This expression is the one we plot in dashed blue in figure 2(c), finding a very good agreement with the numerically obtained D(E) for this energy range.
• (b)  
  Another important difference in the regime 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 β₁ → 2, until they collapse at the critical value where the two Dirac energies merge, i.e. β₁ = 2. At this semi-Dirac point it can be shown that D(E) loses its linear energy dependence and transforms to $D(E)\propto \sqrt{\vert E\vert }$, as highlighted in red in figure 3(c), in accordance to other works [4, 8, 16].
• (c)  
  For β₁ > 2, the bands do not touch anymore, opening thus a band-gap around the Dirac energy, as clearly shown by D(E).

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 ${{\Sigma}}_{e}^{\alpha }(0)=0$ for all β₁, it is guaranteed that a bound state exists at E_BS = 0 (see equation (16)) for the whole anisotropy range. We can then characterize numerically its spatial shape by solving equations (18) and (19) for several β₁ and a fixed system size. The result of this analysis is shown in figures 3(a)–(f) where we plot the bound state shape ${C}_{b,\mathbf{n}}^{A}$ for several β₁ ranging from β₁ = 1 (isotropic case, in panel (a)) to β₁ = 2.1 (gapped situation, in panel (b)), passing through the semi-Dirac situation β₁ = 2 (in panel (e)). From this analysis we can extract several conclusions:

• (a)  
  As it occurs in the case of 3D photonic Weyl points [58, 59], since the density of states maintains its singular nature for β₁ ∈ [1, 2] one can modify the shape of the bound state without altering its power-law dependence. The main change of the spatial shape is due to the modification of the interference factor of ${I}_{\mathbf{K},{\mathbf{K}}^{\prime }}$ due to the movement of the Dirac points K⁽'⁾ with the anisotropy parameter β₁. In fact, the asymptotic decay for large distances still shows a similar 1/r dependence:

  $$\begin{equation}\vert {I}_{{\mathbf{K}}^{(\prime )}}(\mathbf{p})\vert \propto \frac{1}{\sqrt{{v}_{1}^{2}{\tilde{p}}_{1}^{2}+{v}_{2}^{2}{\tilde{p}}_{2}^{2}}},\end{equation} \tag{ 35 }$$

  with a modification compared to the isotropic case of equation (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 figure 4(a) where we plot a line cut of the quasi-bound-state $\vert {C}_{b,(-n,-n)}^{A}\vert$ for several β₁'s. There, we observe how for the β₁'s in that regime, i.e. β₁ = 1, 1.3, 1.6, 1.9, the asymptotic scaling of the quasi-bound 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 equation (26)).
• (b)  
  This situation changes dramatically at the semi-Dirac point β₂ = 2, see figure 3(e), where the interference between ${I}_{\mathbf{K},{\mathbf{K}}^{\prime }}$ 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₂ velocity vanishes. Evidently, a different choice of the anisotropy direction in the hopping model of figure 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 figure 4(a), where we observe how the bound state decay with β₁ = 2 (yellow triangles) is significantly longer ranged than the other ones. In fact, one can make a numerical fit to a power-law decay $\propto 1/{n}^{{\gamma }_{f}}$ obtaining a value γ_f ≈ 0.45 for the system size chosen in figure 4(a), i.e. N = 1000². Through a numerical analysis we demonstrate that this exponent γ_f depends on the system size. This is shown in figure 4(b) where we plot the quasi-bound state cut for β₁ = 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², 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/\sqrt{x}$ dependence resembles the results related to electronic and photonic transport found in the literature [15, 18].
• (c)  
  Finally, when β₁ > 2 the quasi bound-state transforms into a standard one recovering its exponential decay dependence, as clearly shown in purple rhombus in figure 4(a). Remarkably, though, this exponential localization preserves the one-dimensional directional behaviour from the semi-Dirac dispersion, as shown clearly in figure 3(f).

Zoom In Zoom Out Reset image size

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 figure 5 by plotting the numerically calculated excited state population, |C_e (t)|² of a single initially excited emitter coupled to a photonic bath with strength g = 0.1J₁, 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 β₁ parameters chosen in figure 3. Besides, we perform the simulation using a photonic bath with periodic boundary conditions to compare with the results of reference [49]. On top of the dynamics we plot in dashed lines the expected steady-state value, |R₀|², obtained from calculating the overlap with the quasi-bound state given by equation (29) but replacing ω_iso(k) → ω₁(k). From this analysis, we can extract several conclusions:

• (a).  
  For β₁ ∈ [1, 2], and for times in which the photonic excitations do not see the border, i.e. tJ₁/N < 1, all β₁ show a similar non-exponential decay dynamics which resembles the one of the isotropic case. Thus, for an infinite system, N → ∞, we expect then the dynamics to be dominated purely by the non-Markovian dynamics, C_NM(t) in equation (28), induced by the non-analytical behaviour of the density of states around the Dirac energy.
• (b).  
  For β₁ ∈ (1, 2) and long-enough times, i.e. tJ₁/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₀|². As expected, R₀ depends now not only on system size, N, but also on the anisotropy parameter β₁, i.e. R₀ ≡ R₀(N, β₁). However, as we show in figures 5(b)–(e), an important difference arises with respect to the isotropic case: for fixed β₁, R₀ does not monotonically decrease to 0 with N, but rather with some oscillations that depend on the β₁ parameter. The same occurs with the dependence with β₁ for fixed system size. As a general rule, we find that R₀ tends to zero as β₁ gets closer to the critical value, β₁ → 2, as indicated by the larger oscillations found in figures 5(c) and (d).
• (c).  
  At the critical point, β₁ = 2, we find that R₀(β₁ = 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 figure 5(e), and it ultimately comes back coherently to the emitter at sufficiently long times (not shown).
• (d).  
  Finally, for β₁ > 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₀|² ≈ 1, and independent of system size. This is very clear in figure 5(e) where we see the collapse of all the lines for the dynamics of different system sizes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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}^{\alpha \beta }$ appearing in equation (20). For the regime β₁ ∈ (1, 2) one can tune the spatial shape without losing its power-law behaviour, something only found so far for three-dimensional topological photonic baths [58, 59]. For β₁ > 2 one can obtain a directional and longer-ranged interaction ($1/\sqrt{r}$), corrected by an exponential decay length which decreases with the value of β₁. At the critical point β₁ = 2, however, the overlap with the quasi-bound state goes to zero, such that it cannot eventually mediate the interactions ${G}_{ij}^{\alpha \beta }$ for very long times.

For completeness in the study of the possible Dirac-like 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₂ ≠ 0 and β₂ ≠ 1 (see figure 1(b)). The reason is that this simplified model will allow us to tilt the energy dispersion [19–22], and explore the quantum optical phenomena induced by the so-called types II and III Dirac points [25–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 photonic bath Hamiltonian, ${\tilde{H}}_{\text{B}}(\mathbf{k})$, the addition of next-nearest neighbour hoppings enters in the diagonal part h_I(k) ≠ 0 written in equation (11). Thus, the energy dispersion will be given by two bands with energies:

$$\begin{equation}{\omega }_{2,u/l}(\mathbf{k})={h}_{\text{I}}(\mathbf{k})\pm {\omega }_{1}(\mathbf{k}).\end{equation} \tag{ 36 }$$

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 ω₁(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:

$$\begin{equation}{\tilde{H}}_{\text{B}}(\mathbf{K}/{\mathbf{K}}^{\prime }+\mathbf{q})\approx {\tilde{E}}_{\text{D}}+({v}_{0,1}{q}_{1}\mp {v}_{0,2}{q}_{2})\mathbf{1}+\end{equation} \tag{ 37 }$$

$$\begin{equation}+{v}_{1}{\sigma }_{y}{q}_{1}\mp {v}_{2}{\sigma }_{x}{q}_{2},\end{equation} \tag{ 38 }$$

for |q_i | ≪ 1, where ${\tilde{E}}_{\text{D}}$ is the shifted energy of the Dirac points given by:

$$\begin{equation}{\tilde{E}}_{\text{D}}={J}_{2}\left[2{\beta }_{2}+{\beta }_{1}(2-{\beta }_{1}{\beta }_{2})\right],\end{equation} \tag{ 39 }$$

v_1/2 the slopes around the Dirac energy defined in equation (33), and v_0,1/2 the new effective diagonal velocities, which can be shown to be given by:

$$\begin{equation}{v}_{0,1}=0\enspace ,\qquad {v}_{0,2}=2\frac{{J}_{2}}{{J}_{1}}{v}_{2}\left({\beta }_{1}{\beta }_{2}-1\right).\end{equation} \tag{ 40 }$$

As explained in references [19, 33], with these velocities, v_0,i , v_i , one can define a parameter, $\tilde{v}$:

$$\begin{equation}\tilde{v}=\sqrt{{\left(\frac{{v}_{01}}{{v}_{1}}\right)}^{2}+{\left(\frac{{v}_{02}}{{v}_{2}}\right)}^{2}}=\frac{2{J}_{2}}{{J}_{1}}\enspace \vert {\beta }_{1}{\beta }_{2}-1\vert ,\end{equation} \tag{ 41 }$$

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. reference [33] for an explanation of the classification). For $\tilde{v}\in (0,1)$ the energy dispersion is tilted but still has a single-point band-touchings (type I Dirac points). At the critical value $\tilde{v}=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 $\tilde{v} > 1$ the band is over-tilted showing the so-called type II Dirac points characterized by nodal line crossings.

In figures 6(a)–(d) we show an example of this transition for the energy dispersions ω_2,u/l (k) of equation (36). For this plot, we fix J₂ = 0.1J₁, β₁ = 1, and move the anisotropy parameter for positive β₂ ranging from β₂ = 3 to 7.5 in 1.5 steps. From equation (41) and with this choice of parameters, it can be shown that the critical tilting $\tilde{v}=1$ should occur for β₂ = 6 (corresponding to figure 6(c)). However, what happens in this anisotropic model is more intricate. For small enough anisotropy β₂ the energy dispersions are indeed tilted and only touch at a single point corresponding to the Dirac energy ${\tilde{E}}_{\text{D}}$, see figure 6(a). For bigger β₂, though still smaller than the critical value for which $\tilde{v}=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 figure 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, β₂ = 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 figure 6(c)). Finally, for larger β₂, 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.

Zoom In Zoom Out Reset image size

Before considering what happens when an emitter is coupled to this photonic bath, it is instructive to understand first how this transition is manifested in the density of states of the photonic bath, D(E), that we know plays a very relevant role in the emitter's dynamics. In figure 7(a) we study the evolution of D(E) as a function of β₂ for the same parameters than figure 6. We take a range of β₂ ∈ (3, 7.5) in steps of 0.5, which includes the transition through the critical value β₂ = 6, highlighted in red in figure 7(a). There, we observe several features: (i) as expected from equation (39), the Dirac point shifts to larger values as β₂ increases. We highlight this movement in dashed blue lines. There, we observe that for β₂ ⪅ 4, the density of states at the Dirac energy is approximately zero, as expected from type I Dirac points. However, for β₂ ⪆ 4 the density of states starts to acquire a non-zero value. This is more clear at the inset panel of figure 7(a), where we plot the value of the density of states precisely at the Dirac energy, $D({\tilde{E}}_{\text{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 β₂ goes to a critical point, one of the Van-Hove singularities of the density of states approaches and crosses the Dirac point, which is one of the signatures of type III-Dirac points. For β₂ > 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.

Zoom In Zoom Out Reset image size

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 types II and III Dirac points. We start by plotting in figure 7(b) the excited state population of an emitter coupled, with strength g = 0.1J₁ and energy ${\Delta}={\tilde{E}}_{\text{D}}$, to the central A-site of a photonic lattice with the parameters of figures 6(a)–(d). There, we observe how the slow non-Markovian relaxation expected for type I Dirac points is only observed for the smallest β₂ 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({\tilde{E}}_{\text{D}})$. Thus, the presence of other resonant k-modes 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 photonic bath, given by the wavefunctions ${C}_{a/b,\mathrm{n}}^{A}(t)$ of equation (15). This is what we plot in figure 8 by choosing the same parameters of figure 7(b) and looking at the wavefunction at the time $t=0.3\sqrt{N}/{J}_{1}$. We only plot the population in the b-modes since the one of the a sites displays a qualitatively similar behaviour (except for β₂ = 3 where the a modes shows no population). In these panels, we can observe how for β₂ ⪅ 4 (figure 8(a)) the light becomes strongly localized around the emitter, being thus dominated by the quasi-bound state physics discussed in section 5. When the other k-modes become resonant, like in β₂ = 4.5 (figure 8(b)), the population into the photonic bath is a superposition between the quasi-bound state localization and the propagating emission coming from the other modes available. At the critical point (figure 8(c)), the emission becomes strongly directional due to the presence of a Van-Hove singularity around the Dirac point [61–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 β₂ = 7.5 corresponding to the type-II situation (see figure 8(d)), 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 figure 6(d)).

Zoom In Zoom Out Reset image size

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–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 reference [34], especially designed to have a perfect type-III energy dispersion. The model can be written as a bipartite photonic bath as follows:

$$\begin{equation}{\tilde{H}}_{\text{B}}(\mathbf{k})=\left(\begin{matrix}\hfill d(\mathbf{k})+a(\mathbf{k})\hfill & \hfill a(\mathbf{k})\hfill \\ \hfill a(\mathbf{k})\hfill & \hfill a(\mathbf{k})-d(\mathbf{k})\hfill \end{matrix}\right),\end{equation} \tag{ 42 }$$

where d(k) = 2J₁(cos k₁ + cos k₂) and a(k) = 2J₂(cos k₁ − cos k₂). This model has two bands ${\omega }_{u/l}(\mathbf{k})=a(\mathbf{k})\pm \sqrt{a{(\mathbf{k})}^{2}+d{(\mathbf{k})}^{2}}$. Since d(k) vanishes at k₁ ± (∓)k₂ = ±π, and a(k) vanishes when k₁ = ±k₂, 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 figure 9(a) where we plotted ω_u/l (k) for J₂ = 0.3J₁.

Zoom In Zoom Out Reset image size

Apart from having a perfect singular band-touching, the energy dispersion along the straight lines k₁ ± (∓)k₂ = ±π 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 figure 9(b), where we plot D(E) for this model with J₂ = 0.3J₁. This divergent density of states has strong impact in the quantum emitter dynamics, that we plot in figure 9(c), for an emitter coupled to the central A site of a photonic lattice with J₂ = 0.3J₁ and N = 1024² sites. Differently from the critical point dynamics of figure 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 photonic bath, as plotted in figure 9(d). Thus, both the reversible dynamics and directional emission can be considered as smoking guns of the critical Dirac transition.

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–69]. The idea of these proposals is sketched in figure 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–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–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–103], have also been obtained using different laser configurations to generate the optical trapping potentials for the atomic arrays.

Zoom In Zoom Out Reset image size

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 figure 10(a): it is described by a two-site Bravais lattice with primitive vectors ${\mathbf{c}}_{1/2}=d\sqrt{3}/2(\sqrt{3}{\hat{\mathbf{e}}}_{x}\pm {\hat{\mathbf{e}}}_{y})$, and basis vector ${\mathbf{d}}_{\text{bas}}=-{d}_{\text{intra}}{\hat{\mathbf{e}}}_{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 ($\sqrt{3}d$) 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 figure 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 Hamiltonian [87, 96]:

$$\begin{equation}\frac{{H}_{m}}{\hslash }=\sum\limits _{j=1}^{{N}_{\text{A}}}\left({\omega }_{a}-\text{i}\frac{{{\Gamma}}_{a}}{2}\right){\sigma }_{ee}^{j}+\sum\limits _{i\ne j=1}^{{N}_{\text{A}}}\left({J}_{ij}-\text{i}\frac{{{\Gamma}}_{ij}}{2}\right){\sigma }_{eg}^{i}{\sigma }_{ge}^{j},\end{equation} \tag{ 43 }$$

where j is an index running over all atoms in the metasurface (N_A), placed at positions r_j , and ${{\Gamma}}_{a}=\vert {\boldsymbol{\wp }}_{\boldsymbol{a}}{\vert }^{2}{\omega }_{a}^{3}/(3\pi \hslash {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₀(r_i − r_j ) [104, 105]:

$$\begin{equation}{J}_{ij}-i\frac{{{\Gamma}}_{ij}}{2}=-\frac{3\pi {{\Gamma}}_{a}c}{{\omega }_{a}}{\hat{\boldsymbol{\wp }}}_{i}^{{\ast}}\cdot {\mathbf{G}}_{0}({\mathbf{r}}_{i}-{\mathbf{r}}_{j})\cdot {\hat{\boldsymbol{\wp }}}_{j},\end{equation} \tag{ 44 }$$

where ${\hat{\boldsymbol{\wp }}}_{i}={\boldsymbol{\wp }}_{i}/\vert {\boldsymbol{\wp }}_{i}\vert$, and:

$$\begin{equation}{\mathbf{G}}_{0}(\mathbf{r})=\frac{1}{4\pi }\left[\mathbb{1}+\frac{\nabla \otimes \nabla }{{k}_{0}^{2}}\right]\frac{{\text{e}}^{\text{i}{k}_{0}\vert \mathbf{r}\vert }}{\vert \mathbf{r}\vert }.\end{equation} \tag{ 45 }$$

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 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,

$$\begin{equation}{S}_{\mathbf{k}}^{{\dagger}}=\frac{1}{\sqrt{N}}\sum\limits _{n=1}^{N}\sum\limits _{m=1}^{2}{\sigma }_{eg}^{n,m}{\text{e}}^{\text{i}\mathbf{k}\cdot {\mathbf{r}}_{n}},\end{equation} \tag{ 46 }$$

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 eigen-energies of the Bloch modes from the above Hamiltonian then reduces to diagonalising the following matrix [96]:

$$\begin{align}{\mathbf{M}}_{\mathbf{k}}^{\alpha \beta ,\mu \nu }& =\left({\omega }_{a}-\text{i}\frac{{{\Gamma}}_{a}}{2}\right){\delta }_{\alpha \beta }{\delta }_{\mu \nu }\\ & \quad -\frac{3\pi {{\Gamma}}_{a}c}{{\omega }_{a}}\left[\sum\limits _{m=1}^{2}\sum\limits _{{\mathbf{R}}_{n}\ne 0}{\text{e}}^{-\text{i}\mathbf{k}\cdot {\mathbf{R}}_{n}}{G}_{0}^{\alpha \beta }({\mathbf{R}}_{n}){\delta }_{m\mu }{\delta }_{m\nu }\right.\\ & \quad +\sum\limits _{{\mathbf{R}}_{n}}{\text{e}}^{-\text{i}\mathbf{k}\cdot {\mathbf{R}}_{n}}{G}_{0}^{\alpha \beta }({\mathbf{R}}_{n}+{\mathbf{d}}_{\text{bas}}){\delta }_{1\mu }{\delta }_{2\nu }\\ & \quad +\left.\sum\limits _{{\mathbf{R}}_{n}}{\text{e}}^{-\text{i}\mathbf{k}\cdot {\mathbf{R}}_{n}}{G}_{0}^{\alpha \beta }({\mathbf{R}}_{n}-{\mathbf{d}}_{\text{bas}}){\delta }_{2\mu }{\delta }_{1\nu }\right],\end{align} \tag{ 47 }$$

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 ${\omega }_{\mathbf{k}}-\text{i}\frac{{\gamma }_{\mathbf{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].

Figures 10(b) and (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) and (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'MK 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) and (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 references [67–69], in order to probe the physics explored along this manuscript one can place additional impurity atoms (as depicted in figure 10(a)) with their optical transitions ω₀ tuned to the spectral region where these non-trivial band-crossings occur, e.g. using Raman-assisted processes.

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-types 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–86] to be further engineer the quantum emitter dynamics.

AGT acknowledges support from CSIC Research Platform on Quantum Technologies PTI-001 and from Spanish project PGC2018-094792-B-100(MCIU/AEI/FEDER, EU). AGT also acknowledges initial discussions of this project with U Zabaleta. PAH acknowledges funding from Fundação para a Ciência e a Tecnologia and Instituto de Telecomunicações under project UIDB/50008/2020 and the CEEC Individual program from Fundação para a Ciência e a Tecnologia with reference CEECIND/02947/2020. M.B.P. acknowledges support from the Spanish Ministerio de Ciencia e Innovación (PID2019-109905GA-C2) and from Eusko Jaurlaritza (IT1164-19 and KK-2019/00101). AGT acknowledges finantial support from the Proyecto Sinérgico CAM 2020 Y2020/TCS-6545 (NanoQuCo-CM).

All data that support the findings of this study are included within the article (and any supplementary files).