Non-Markovian Quantum Optics with Three-Dimensional State-Dependent Optical Lattices

Quantum emitters coupled to structured photonic reservoirs experience unconventional individual and collective dynamics emerging from the interplay between dimensionality and non-trivial photon energy dispersions. In this work, we systematically study several paradigmatic three dimensional structured baths with qualitative differences in their bath spectral density. We discover non-Markovian individual and collective effects absent in simplified descriptions, such as perfect subradiant states or long-range anisotropic interactions. Furthermore, we show how to implement these models using only cold atoms in state-dependent optical lattices and show how this unconventional dynamics can be observed with these systems.

Quantum emitters coupled to structured photonic reservoirs experience unconventional individual and collective dynamics emerging from the interplay between dimensionality and non-trivial photon energy dispersions. In this work, we systematically study several paradigmatic three dimensional structured baths with qualitative differences in their bath spectral density. We discover non-Markovian individual and collective effects absent in simplified descriptions, such as perfect subradiant states or long-range anisotropic interactions. Furthermore, we show how to implement these models using only cold atoms in state-dependent optical lattices and show how this unconventional dynamics can be observed with these systems.

I. INTRODUCTION
Initially motivated by overcoming the figures of merit of standard quantum optical setups, there exists a growing interest in integrating quantum emitters with nanophotonic structures [1][2][3][4][5][6][7][8][9]. The confined photons in these structures display highly structured energy dispersions, whose interplay with the dimensionality induces qualitatively new phenomena in the individual and collective quantum emitter (QE) dynamics.
To name a few examples, these structured photonic reservoirs give rise to the emergence of atom-photon bound states [10][11][12][13][14][15], novel super/subradiant phenomena [16][17][18][19][20][21][22][23] or coherent non-local interactions [24][25][26], among others. On top of that, these structured baths also display non-analytical spectral regions in their density of states, e.g., band edges, where nontrivial dynamics emerge [27][28][29][30][31][32][33] beyond what could be observed in unstructured baths [34]. This phenomena is typically referred to as non-Markovian dynamics, since standard perturbative treatments like Born-Markov master equations [35] fail. The failure of conventional methods poses a great theoretical challenge for the characterization of these systems, which together with the possibility of observing exotic phenomena, makes the study of such non-Markovian dynamics still nowadays a very active research arena (see Ref [36] and references therein for a recent review on the subject).
The exciting prospects of finding qualitative different phenomena compared to standard quantum optical systems have also attracted the attention of other communities beyond the optical regime, like circuit QED [37][38][39] or atoms in statedependent optical lattices [40][41][42], allowing one to enter in a completely new parameter regime as compared to the optical implementations. As an illustration of the potential of these new platforms, they have already succeeded observing singlephoton bound-state physics associated to 1D band-edges in frequency [37] and time domain [42], a long-standing prediction in the optical regime [11,12,27], which took more than 25 years to be observed [43].
So far, most studies of these unconventional setups are devoted to one and two-dimensional baths, while the physics emerging from three-dimensional (3D) structured baths remains largely unexplored due to the theoretical and experimental challenges to study them. From the theoretical side, * alejandro.gonzalez-tudela@mpq.mpg.de 3D structured baths are challenging due to the large Hilbert space to be considered. This is why most theoretical studies of 3D baths use simplified models, e.g., isotropic dispersions [11,12,40,41]. From the experimental point of view, they also represent an outstanding challenge, since the integration and addressing of QEs embedded in three dimensional systems will be extremely difficult in nanophotonics or circuit QED implementations. For all those reasons, the characterization and implementation of these unconventional 3D structured baths represent one of the current frontiers of quantum optical studies.
In this manuscript, we both i) analyze the individual and collective dynamics of QEs coupled to several paradigmatic 3D baths, discovering non-Markovian phenomena qualitatively different from other reservoirs; and ii) we propose a setup based on state-dependent optical lattices where to observe them [40,41]. In particular, we consider several nearestneighbour tight-binding baths with different geometries: cubic simple (CS), body-centered cubic (BCC), face-centered cubic (FCC) and diamond lattice [44]. All these baths show distinctive features in their density of states, which translates in very different individual and collective QE dynamics. We discover a plethora of effects absent in simplified descriptions [11,12,40,41], such as: • Long-time reversible dynamics for a single QE spectrally tuned at the middle of the band for CS geometries; • Directional emission and perfect subradiance with QEs in BCC geometries; • Robust 3D bound states leading to anisotropic coherent interactions between QEs in FCC lattices; • Purely long-range coherent interactions between QEs in diamond geometries; among many other effects. We use both numerical and analytical techniques, which allows us to extract the scalings with the system parameters of these features analytically. Regarding the implementation with state-dependent optical lattices, i) we provide the laser configurations required to implement these structured baths with optical dipole traps, ii) we analyze the estimated timescales of the dynamics, and finally iii) we study the deviations from the idealized nearest neighbour description of the bath dynamics, focusing on their impact in the bath density of states. Overall, we show how the observation of the discovered non-Markovian 3D phenomena can potentially be observed with three-dimensional state-dependent optical lattices. The outline of this manuscript reads as follows: in Section II we discuss the general features of the setup that we study along this manuscript. Here, we also explain how to implement these models using cold atoms in state-dependent optical lattices, since it is currently the more realistic platform where to implement them. Then, to make the manuscript selfcontained, in Section III we briefly review the theoretical techniques we use to study the problem. In Section IV, we give an overview of the properties of the different baths considered along the manuscript, explaining both the lattice geometry and the expected density of states in the thermodynamic limit. Then, in Sections V-VIII we calculate the quantum emitter dynamics for the aforementioned baths of interest, focusing on non-Markovian phenomena and emphasizing the differences with respect to other types of structured reservoirs. After presenting the emergent phenomena, in Section IX we go back to the implementation discussion and show the laser configuration required to implement the bath Hamiltonians studied along this manuscript in optical lattices. Furthermore, we calculate the band structure of the proposed lattices to estimate how well they describe the toy models considered along the manuscript, which we characterize by looking at the density of states of the obtained model. Finally, in Section X we summarize the main results of the manuscript, and point to future work directions.

II. GENERAL QUANTUM OPTICAL SETUP AND ITS IMPLEMENTATION WITH STATE-DEPENDENT OPTICAL LATTICES
In this Section we describe the general model that we consider along this manuscript and its natural implementation with cold atoms in state dependent optical lattices.

A. General model
We are interested in the dynamics of N e QEs, which we describe as two-level systems, {|g j , |e j } N e j=1 , whose intrinsic dynamics is given by the following Hamiltonian (we useh ≡ 1 along the manuscript): where we use the notation σ j αβ = |α j β | j for the spin operator of the j-th QE, and ω e is the transition frequency of the QE optical transition that couples to the bath modes.
For the bath description, we take a discretized version with a finite number of bosonic modes. We consider baths of linear size N, such that the total number of sites is N × N × N. For the baths that we consider along this manuscript, one has to distinguish two situations, namely, the single and two-band configurations. The single-band model corresponds to a configuration where the bath geometry can be described as a simple Bravais lattice [44]. This means that the position of each bosonic mode, with creation/annihilation operators a † n , a n , is given by three integer numbers: n = (n 1 , n 2 , n 3 ) = ∑ 3 i=1 n i c i , with n i ∈ (0, . . . , N − 1), which describes the displacement in terms of the primitive vectors, c i , expanding the lattice in real space. Unless stated otherwise, we assume that each bosonic mode only couples to its nearest-neighbour, with strength J, that we will take as the unit of energy. The number and positions of the nearest neighbours depend on the particular bath geometry, as we will show in Section IV. Using all these assumptions, the bath Hamiltonian can be generally written as follows: where ω a is the bosonic mode frequency which we assume to be equal for all the bath sites. In the last step we have assumed periodic boundary conditions for the bath modes and introduced the operatorsâ k = 1 N 3/2 ∑ n e ik·n a n , which diagonalize the bath Hamiltonian in momentum space. Notice, we are using the hat notation,·, to distinguish the bath operators in position/momentum space. The momenta k are also described by three numbers k = (k 1 , k 2 , k 3 ) = ∑ i k i d i , being d i are the primitive vectors of the reciprocal lattice, which satisfy: c i · d j = δ i j . As we show in Section IV, this single-band situation describes the cases of the CS, BCC and FCC lattices.
There are bath geometries, however, which cannot be described as simple Bravais lattices [44]. For example, the diamond lattice is described by two interspersed FCC lattices (we explain this in detail in Section IV). In those situations, one needs to upgrade the single-band model to a more general case which allows one to capture the exact dynamics of the bath. For this manuscript, we consider the simplest upgrade, namely, a two-band model consisting of two discretized N × N × N lattices (denoted as A/B sublattices), with associated bosonic modes, a n , b n , with the same energy ω a , interacting through nearest neighbour coupling as follows: Hereû k ,l k are the operators that diagonalize the bath Hamiltonian in momentum space, with eigenenergies: ω u,l (k) respectively. This is obviously not the most general two-band configuration, since it may also occur that there exist hoppings between the AA/BB sites, or there is some energy off-set between sublattices. However, since the idealized diamond lattice is described by a Hamiltonian like in Eq. 3, we restrict ourselves to this situation in our manuscript.
In both the single and two band models, it is convenient to move to a rotating frame with ω a , where H B → H B − ω a ∑ n a † n a n (−ω a ∑ n b † n b n ) and H S = ∆ ∑ j σ j ee , with ∆ = ω e − ω a . Figure 1. Two dimensional view of the implementation of QED Hamiltonians with state-dependent optical lattices: two internal atomic states, labeled as f and a, are trapped in very different trapping potentials. The f /a atoms are trapped in a deep/shallow lattice playing the role of QEs/photon-like modes. They are connected through a two-photon Raman transition, or direct optical coupling for Alkaline-Earth atoms, at a rate g.
The last ingredient is the system-bath interaction. For the sake of simplicity, we just assume a local coupling between the QEs and the bath sites, since the spread of the QE wavefunction ultimately depends on the particular implementation considered. Mathematically, this local coupling assumption means that the j-th QE couples only to a bosonic mode at site n j , such that the general light-matter interaction Hamiltonian generally reads: where c = a, b denotes a general bath operator, depending on whether the atom couples to an A/B site. Notice, that we have only written the excitation conserving terms in the lightmatter Hamiltonian of Eq. 4. This is a safe assumption in the optical regime [45] and it will be an exact description in the cold-atoms implementation that we discuss in the next Section.

B. General implementation with state-dependent optical lattices
In this subsection, we give a brief explanation on how to implement these structured 3D quantum optical models using cold atoms in state-dependent optical lattices. It will be mostly based on the proposal of Refs. [40,41], which has been recently implemented with Rb atoms for 1D reservoirs [42].
The idea to mimic the quantum optical scenario with cold atoms is sketched in Fig. 1: we require two atomic internal states, that we label as f and a, trapped in two very different potentials. One of them, e.g., the f state, is trapped in a very deep potential such that its hopping to nearest neighbours is suppressed. We assume the trapping potential to be harmonic, with trapping frequency ω e , and restrict to the dynamics of the lowest motional state and in the collisional blockade regime where no more than one atom can be the trapped in each site.
In this regime, the creation atomic operators f † n ≈ σ n eg can be approximated by spin operators and therefore play the role of the QEs in our problem.
The other internal state, a, is trapped in a shallower potential, with trapping frequency ω a , such that they can hop to their nearest neighbours sites at a rate J. In the original proposal [40,41], the atoms in this internal state were assumed to be free propagating particles with isotropic dispersion ω(k) ∝ |k| 2 . Here instead, we consider that the a-modes feel the underlying potential, leading to a non-trivial energy dispersion ω(k) which depends on the geometry of the lattice considered.
Finally, the coupling between the two internal states can be obtained through two-photon Raman transitions [40][41][42], or even through a direct optical transition in the case of Alkaline-Earth atoms [46][47][48], where there exists an optically excited state with lifetimes even longer than seconds. To mimic the dynamics of one or few excitations, we assume the laser addressing to be local as depicted in Fig. 1. As shown in the original proposal [40,41], these laser fields induce hopping between a/f states, where the Hamiltonian reads: where g k contains a k-dependence emerging from the finite size of the Wannier functions. This finite size of the wavefunction gives a natural cut-off in momentum space, whose main effect was shown to be a renormalization of the frequencies [40,41]. In the limit of very localized Wannier functions, and local laser addressing,ḡ k ≈ g/N 3/2 , such that H las ≈ H int .
Remarkably, H las conserves the total number of excitations of the system irrespective of g, unlike what happens in the optical regime. This allows one to explore parameter regimes very difficult to access with other platforms including nonperturbative ones, such as g > J.
Finally, let us mention other advantages like: i) the possibility of single-site detection and addressing [49,50] to, e.g., observe in-situ the bath dynamics; ii) low decoherence rates (<Hz), compared to the Hamiltonian timescales g, J ∼ 10 KHz; iii) importantly for this work, the possibility to simulate 3D models, something very difficult to implement with other platforms; and iv) the possibility to change ω(k) just by changing the laser configuration. For example, simulating CS baths can be done by simply sending 6 counterpropagating lasers in the X/Y/Z directions with orthogonal polarizations [51], which results in a separable optical potential V (R) ∝ cos 2 (x) + cos 2 (y) + cos 2 (z). The configurations required for the rest of the lattices are less trivial. Thus, we leave the discussion on how to implement these optical lattices, and their resulting band structure, to Section IX.

III. THEORETICAL TOOLS
In this manuscript we mainly consider the spontaneous emission of a single excitation into the bath being emitted from one or few QEs, as a first step to unravel the dynamics

BCD
Branch cut detour contribution: appearing for every detour due the non-analyticity in the Green Function induced by the branch cut. Figure 2. (a) Example of a possible contour of integration to calculate C α (t) as defined in Eq. 9; One closes the contour with a semiarc in the lower complex plane. In the continuum limit, G(z) develop nonanalyticities and branch cuts that one must avoid with detours (in yellow) to be able to apply Residue Theorem. In this example there are four detours, although in general can be more or less (at least two for finite bands). Real (BS) and complex poles (UP), depicted in blue/red, will also contribute to the dynamics. (b) Summary of acronyms of the different contributions to the dynamics of C α (t).
emerging from these unconventional 3D baths. This allows us to characterize both the individual and collective response in the linear regime. Importantly, the system-bath Hamiltonian, H = H S + H bath + H int , conserves the total number of excitations, such that we restrict to the single-excitation subspace where the calculations are more accessible. To calculate the dynamics we use two different and complementary approaches, which we describe in what follows.

A. Analytical techniques
The first one is an analytical approach based on the resolvent operator technique [45,52]. This technique consists in calculating the Laplace transform of the time-evolution operator, U(t) = e −iHt , which is the so-called Green-Function: If {|α } is complete basis expanding H and we are interested only in the dynamics of a particular subspace, e.g., P = β |β , one can show that the Green-Function in this subspace can be calculated as: where P is the projection in the subspace we are interested in, and Σ(z) is the self-energy of the problem, that contains the effect of the interaction with the bath: where Q = 1 − P. This method has two challenges: first, calculating the resolvent operator, which for 3D baths is given, in the thermodynamic limit, by a 3D integral. For example, when only a single QE is coupled to a single-band bath it reads: where in the last equality we used: 1 , which is the usual prescription to go to the thermodynamic limit. An important asset of this work is that we have analytical formulas for Σ e (z) for our baths of interests, which allows us to extract valuable information from the dynamics, such as scaling of decays, bound-state energies,. . .
The second challenge consists in moving from the Laplace to the time domain to get the dynamics. This can be done by calculating the inverse Laplace transform, which is done by solving an integral of the type: where C α (t) is the probability amplitude of a given state |α . As sketched in Fig. 2(a), one possibility to solve these integrals is to find a closed a contour that contains the integral of Eq. 9 and apply complex integral techniques, such as Residue Theorem. However, in the continuum limit, depending on ω(k), Σ(z), and consequently G(z), may develop non-analytical regions and posses branch cuts that have to be avoided with detours in the contour of integration (depicted in yellow in Fig. 2). Depending on ω(k) one must introduce as many detours as required to guarantee the analyticity of G(z) in the whole closed region to be able to apply Residue Theory. Using this technique, one is able to decompose the dynamics of C α (t) in different contributions (summarized in Table 2(b)): • Bound State (BS) contributions. They arise from real poles, E BS ∈ R, of G(z), whose energy lies out of the continuum. They give a contribution to the dynamics of the form: R BS e −iE BS t , where R BS is the residue of the pole. The choice of the BS notation is because the origin of these real poles are the emergence of photon bound states which localize around the QEs [11].
• Unstable pole (UP) contribution. When taking the detours to avoid non-analytical regions in the domain of integration, one might need to analytically continue G(z) to other Riemann sheets [52]. These analytical continuations may show complex poles, z UP , with Rez UP ∈ ω(k), and Imz UP = 0 which also contribute to the dynamics as: R UP e −iz UP t , being R UP the residue associated to them.
• Branch cut detour (BCD) contribution. As we explained, the branch cuts of G(z) will force us to take detours at certain frequencies E BC . Typically, only bandedge detours are considered in the literature. However, we will see along this manuscript that they can also emerge in the middle of the band [18,53]. Since at both sides of the detour one must generally use a different analytical continuation of G e (z), this contribution can be generally written as: where E ± BC = lim ε→0 E BC ± ε, where the definition of G(z) must be adapted depending on the value of z.
Even though each contribution will be calculated numerically, the separation in different terms gives valuable information about the underlying contribution dominating in each parameter regime. Moreover, in some situations these formulas can be used to obtain asymptotic scalings at long times or perturbative couplings, that is, g J.

B. Numerical techniques
The most straightforward method is to solve directly the dynamics of the total system+bath Hamiltonian, that is, e −iHt |Ψ 0 , by trotterizing the time and using split methods. This means that we apply the evolution of the bath, H B , in kspace, whereas the one of H int + H S in position space, as we explained in Ref. [53].
Another option, very useful for 3D models, is to work in the space of frequencies by discretizing ω(k) in steps of dω. Like this, one can group the evolution of the bath Hamiltonian in different isofrequency subspaces [53]: ω n ,αâωn,α . (11) where N ω is the total number of frequency modes considered, i.e., N ω = maxω(k)−minω(k) dω , and ω n = minω(k) + ndω. The index α of the frequency creation/destruction operatorŝ a † ω n ,α /â ω n ,α denotes the degeneracy of the modes at ω n . This index runs from 1 to N (ω n ), where N (ω n ) is the number of modes at the frequency step ω n .
To illustrate how it simplifies the evolution, let us consider a single QE coupled to the bath, and choose the convention where the index α = 1 corresponds to the symmetric combination of the k modes, which is the one in this case that couples the QE dipole transition. With this choice, the interaction Hamiltonian reads: As the interaction Hamiltonian only couples σ eg to the α = 1 indices, for each ω n , one can neglect the dynamics of the N (ω n ) − 1 modes that are uncoupled from the QE. With this method, we have empirically observed one can obtain the dynamics up to a given time, T 1/dω. This means, the longer the time we want to explore, the more frequency modes we need to accurately describe the dynamics.
C. Dynamics in the Markovian regime.
To conclude, it is instructive to remind the reader the predictions for the dynamics obtained in the Markovian/perturbative regime for the problems that we consider along the manuscript. First, let us mention that the Markovian results are easily recovered in our formalism just by replacing Σ(E + i0 + ) ≈ Σ(∆ + i0 + ) in Eqs. 7-9, and assuming that the only contribution to the integral of Eq. 9 is coming from a single pole: z = ∆ + Σ(∆ + i0 + ). This approximation predicts: • For a single QE initially excited, the probability amplitude of the excited state is given by: where δ ω M and Γ M 2 are the real and imaginary parts of the self energy Σ e (∆ + i0 + ), given by the formula: Notice, this approximation agrees with the prediction of the Fermi's Golden Rule, which predicts an exponential relaxation of the population |C e (t)| 2 = e −Γ FGR t , with a decay rate given by: Γ FGR = 2πg 2 D(∆), where D(∆) is the density of states of the bath, which in the thermodynamic limit reads: • When more than one QE is coupled to the bath, there are two different situations. When we start in an initial state which couples only to a collective bath mode A α , the dynamics is analogous to the single QE case, but with a modified self-energy Σ α (∆ + i0 + ) which now will contain extra terms. This will be the case of two QEs initialized in the symmetric/antisymmetric superposition because these two states coupled to orthogonal bath modes, as long as ω(k) = ω(−k) [53]. In that case, the dynamics of these states is governed by a modified self-energy, which now depends on the relative distance position between the QEs, n 12 , as follows Σ ± (z; n 12 ) = Σ e (z) ± Σ 12 (z; n 12 ) for the symmetric/antisymmetric QE superposition, with: • The other situation occurs when we start in a state that couples to many bath modes, e.g., for the situation with two QEs with one of them initially excited. In that case, one needs to consider the coupling to both the symmetric/antisymmetric bath modes [53]. For example, when only the first is excited, the dynamics of the two QEs reads: which predicts a coherent exchange of excitations at rate J M,+ − J M,− , exponentially damped by Γ M,± .

IV. BATH PROPERTIES: GEOMETRIES AND DENSITY OF STATES
After having introduced the general setup, implementation, and the techniques to deal with these problems, let us present the reservoirs that we consider along this manuscript. They will be mainly nearest neighbours tight-binding models with different geometries well-known in the condensed matter context [44]. However, given the importance that they will have for the discussion of the results, we present in this Section their primitive vectors, energy dispersions, and their corresponding density of states, D(E), in the thermodynamic limit. The latter is particularly relevant since it gives us hints on where non-perturbative dynamics is expected due to the appearance of non-analytical behaviour. We have chosen four different models with qualitative differences in their density of states, as we show in Fig. 3, and therefore, where we expect to obtain different dynamical features.

A. Cubic Simple (CS) lattices
We start describing the most simple 3D structured bath, namely, the one given by a CS geometry depicted in Fig. 3(a). In this case, considering the length of the natural cube L as the unit of length, the primitive vectors spanning the lattice are: c 1,2,3 =ê x,y,z . Each lattice site has 6 nearest neighbours at positions: ±c 1,2,3 . This can be shown to lead to an energy dispersion: where k = (k 1 , k 2 , k 3 ) is the momenta given with respect to the primitive lattice coordinates, d i , which in the CS lattice are: d 1,2,3 = 2πê x,y,z . Notice, the reciprocal space for a CS lattice is another CS lattice. It can be easily shown that this band expands from [−6J, 6J], and its corresponding density of states, shown in Fig. 3(e), has several spectral regions which are candidates to display non-perturbative dynamics (in shaded yellow in the figure). For example, at the bandedges the density of states is continuous but with discontinuous derivative. It is easy to show that around the upper/lower edge the density of states is isotropic, e.g., in the lower one ω CS (k)/J ≈ −6 + |k| 2 , which leads to: for E −6J. This is a feature of most 3D band-edges, which is also captured by simplified isotropic models considered in the literature [40,41]. However, the CS lattice already shows some distinctive features not captured by isotropic models, namely, at E = ±2J, D(E) has again a discontinuous derivative but with a finite value of the density of states. As we show in Section V, this will lead to non-Markovian relaxation dynamics different from other types of reservoirs.

B. Body-Centered Cubic (BCC) lattices
The next lattice that we consider is the BCC lattice, depicted in Fig. 3(b), which is characterized by having an extra bosonic mode in the center of the unit cube. This is also a simple Bravais lattice, which can be expanded by choosing the following primitive vectors: Each bosonic mode has 8 nearest neighbours, which position written in the basis of primitive vectors is given by: In this case, the reciprocal vectors are given by: d 1 = 2π(ê x +ê y ), d 2 = 2π(ê x +ê z ), and d 3 = 2π(ê y +ê z ). Diagonalizing the bath Hamiltonian in the momentum space, one arrives to the following energy band dispersion: (20) This band expands from [−8J, 8J], as shown in the corresponding density of states of Fig. 3(f). Apart from the standard scaling D(E) ∝ √ E at the band-edges, due to the isotropic character of the dispersion at these points, this bath has another prominent feature at E ≈ 0, where D(E) diverges. In Section VI, we will prove that this is a square logarithmic divergence, D(E) ∝ ln(E) 2 , which translates into prominent features in the QE individual and collective dynamics.

C. Face-Centered Cubic (FCC) lattices
The other simple Bravais lattice that we consider is the FCC lattice, depicted in Fig 3(c), which is obtained when there is one extra bosonic mode in the center of each side of the unit cube. The three primitive vectors that expand this lattice are: c 1 = (ê x +ê y )/2, c 2 = (ê x +ê z )/2, and c 3 = (ê y +ê z )/2, whereas the three reciprocal ones are: d 1 = 2π(ê x +ê y −ê z ), d 2 = 2π(ê x −ê y +ê z ), and d 3 = 2π(−ê x +ê y +ê z ). Notice that the reciprocal space of the FCC lattice is a BCC one and viceversa. In real space, each bosonic mode of the FCC lattice interacts  (e) The corresponding density of states of this bath, plotted in Fig. 3(g), expands from [−12J, 4J] and show remarkable qualitative differences with respect to the other 3D baths. First, the upper/lower band edges have very different features, unlike the other type of reservoirs. Whereas the lower edge has the typical D(E) ∝ √ E scaling of 3D reservoirs, the upper edge has a divergence, which we will show in Section VII to scale as ∝ ln(4J − E) 2 , for E 4J. Apart from enhancing the decay rates at this edge, this divergence has important consequences in the emergence of robust atom-photon bound states [11,12].

D. Diamond lattice
Finally, we consider a bath geometry which cannot be described as a simple Bravais lattice, namely, the diamond lattice. As we plot in Fig. 3(d), this lattice can be constructed by two interspersed FCC lattices, whose modes are denoted by a n , b n , displaced by a vector f = 1/4(1, 1, 1). To calculate its corresponding band structure, the relevant information is that each bosonic mode in the A lattice interacts with four nearest neighbours from the B lattice at sites: (0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1), and vice-versa. This results into a bath Hamiltonian written in momentum space: is the coupling between the A/B lattices in k-space. The Hamiltonian can be easily diagonalized by introducing the following operators:û which represent the annhiliation operator of the upper/lower band modes. The resulting eigenenergies, ω u/l (k) = ±| f (k)|, have a close connection with the energy dispersion of the FCC lattice: From here, it is very easy to read that the k-points that give rise to the upper edge of the FCC lattice will be the ones where the upper/lower band of the diamond lattice touch, since ω u,l (k) ≡ 0 at these points. The global upper/band edges In the middle of each band, at E = ±2J there are some kinks at the density of states, similar to the ones obtained for the CS and FCC lattice. The main qualitative difference here is the appearance of a singular band gap point at the position where the two bands touch, i.e., E = 0. This is the 3D analogue of the Dirac cone in 2D lattices [54], and it will be the source of qualitatively new features in the QE dynamical response.
After having explored how the structure of the different baths give rise to non-analytical features in their corresponding density of states, in the next sections we study the exact individual and collective QE dynamics emerging from them. We will emphasize the features that are different from reservoirs of lower dimensions or simplified descriptions to the problem.

V. QUANTUM DYNAMICS IN CUBIC-SIMPLE BATHS: LONG-TIME REVERSIBLE DYNAMICS
In this Section we analyze the QE dynamics emerging from CS baths. Since we use the same procedure to study the different baths, we explain it here explicitly to guide the reader in the discussion: • We always start out by analyzing the relaxation of an initially excited QE. We first give the analytical expression of the self-energy Σ e (z). Then, we show how to analytically continue to the whole complex plane to perform the integral of Eq. 9 distinguishing the different contributions.
• Once we have analyzed the exact dynamics of C e (t), we study its emission into the bath and spot the parameters to observe interesting collective phenomena such as super/subradiance or coherent exchange of excitations. Then, we analyze the exact dynamics of several QEs interacting with the bath focusing only on those regions (if any).

A. Single QE
By using the fact ω(k) = ω(−k), the single QE self-energy of the CS bath can be written as: . (25) An analytical expression of this function can be found [55] in terms of elliptic integrals [56], which reads: where the functions ξ (z) and m(z) read: and where K(m) is the complete elliptic integral of the first kind: The complete elliptic integral is real as long as m < 1 and has a branch cut along Rem ∈ [−1, ∞). Evaluating this function above the real axis: 2 , one obtains both the density of states of Fig. 3(a) (up to a factor is the same as Γ e (E)), and the Lamb-shift δ ω e (E), plotted in dashed red in the same figure.
To integrate C e (t) using Eq. 9 and the Residue Theorem, one needs to avoid the non-analytical regions of Σ e,CS (z) by taking detours in the contour of integration. In this particular case, one possibility is to take four detours at energies E = −6J, −2J, 2J, 6J. As we explained in Section III, when taking these detours it is possible that one needs to analytically continue the function to other Riemann sheets and, therefore, change the expression of Σ e,CS (z). In this case, we will need to define the Σ e,CS (z) in the complex lower plane in six different regions: • The definition of Eqs. 26-29 is valid for regions where Rez ∈ (−∞, −6J) and [6J, ∞), which we label as regions I and VI.
• When −6J < Re(z) < −2J and 2J < Re(z) < 6J, which we denote as regions II and V, the 1 − 36J 2 /z 2 of ξ (z) becomes complex. Thus, in order to analytically continue to another Riemann sheet here, one needs to • When −2J < Re(z) < 0/0 < Re(z) < 2J, which we label as regions III/IV respectively, both square roots 1 − 36J 2 /z 2 and 1 − 4J 2 /z 2 become complex such that one needs to change the sign of both of them in ξ (z) to go to the other Riemann sheet. However, this introduces one complication: with the new definition of ξ (z), the argument of the elliptic integral crosses its branch cut when Rez = 0. Thus, one needs to adapt the definition of Σ e (z) for regions III and IV. In particular, if Im(m(z)) < / > 0, then Σ e,CS (z) is defined in the same way as in Eqs. 26, whereas if Im(m(z)) > / < 0, then If one does not change the definition of K(m), function Σ e,CS (z) will be discontinuous in the detours at E = ±2J.
With this piecewise definition of Σ e,CS (z), we can integrate C e (t) using Eq. 9 separating its different contributions: BS, UP and the ones coming from the BCDs (see Fig. 2(b)). To illustrate the different dynamical regimes that appear in this reservoir, we consider a fixed g/J and plot the (absolute value) of the weight of the different contributions as a function of ∆/J, that we show Fig. 4(a). Let us now explain in detail the different regimes that we observe: Outside of the band. For |∆/J| 6J, the dynamics is dominated by the BS contribution, i.e., |R BS | ≈ 1 (in solid black in Fig. 4(a)). This is expected since there are no bath modes which can lead to the relaxation in the bath and it also happens in lower dimensional baths. As ∆/J gets closer to the band-edge, the BS contribution starts decreasing until a critical value of ∆ crit where it suddenly goes to 0. This is different from lower dimensional baths, where the BS are more robust since they survive for all values of ∆. The origin of this disappearance of the 3D BS can be traced back to the finite value of δ ω e (±6J) at the band edge (see Fig. 3(a)), compared to the divergent behaviour of δ ω e (E) for 1D and 2D [53,57]. Intuitively, the bath can only push the energy of the BS out of the band up to a critical value, which actually can be calculated to ∆/J = 0 (red). In solid lines we plot the result from the complex integration of Eq. 9, and in markers a numerical evolution of the complete Hamiltonian using discretized frequency space. be: for the upper/lower edge respectively. The dynamics at these band-edges is then dominated by a combination of the BS, and the BCD contributions which gives a subexponential decay of the population. This is a natural feature of 3D isotropic bandedges that was already predicted in Refs. [11,15,40,41]. For this reason, we will not discuss it further and rather focus on the qualitatively new regimes emerging from our structured reservoir.
Inside of the band. In the middle of the band the dynamics is generally dominated by the UP contribution, that is, the one due to the contribution of complex poles obtained by using the analytical continuation of Σ e,CS (z) in the different regions like Figure 5. (a-c) Bath population in real space, |C R | 2 , emitted from a single QE coupled to a CS bath with g/J = 0.1 obtained numerically with a finite bath with linear size N = 2 7 and for a time T J = N we defined before. As it also occurs for other reservoirs, in the strongly interacting regime (g ∼ J) the imaginary part of these UPs deviate significantly from the perturbative prediction, as shown in Fig. 4(b).
Beyond this non-perturbative renormalization of the lifetimes, two non-trivial spectral regions emerge in the middle of the band at energies ∆ ≈ ±2J, because of non-analytical behaviour of Σ e,CS (z) around them. This behaviour has two consequences: i) there are two extra BCD contributions in the middle of the band (plotted jointly in purple rhomboids in Fig. 4(a)); ii) there are two critical spectral regions where two UPs can contribute simultaneously to the dynamics. Although this looks similar to the behaviour emerging from 2D Van Hove singularities [18,53], it also shows different features: • For example, the δ ω e (E) in 2D reservoirs experiences discontinuous jump from positive to negative values, which creates a symmetric region around the 2D di-vergence where the UP contributions coexist. In the 3D case, the δ ω e (E) is continuous, with a value of δ ω e (±2J) ≈ ∓0.321g 2 /J, and only has a discontinuous derivative. Consequently, this makes that while the UP from region II/V can contribute for |∆| < 2J, the UP contribution of regions II and III only survive for ∆ ∈ (−2J, 2J). Thus, the coexistence regions of the UPs appear in this case at ∼ (−2J, −2J + 0.321g 2 /J) and (2J − 0.321g 2 /J, 2J). In these regions, the dynamics shows overdamped oscillations followed by slow relaxation dynamics, as shown in Fig. 4(c) for ∆ = −2J, similar to the one reported in 2D reservoirs [53]. The origin of the oscillations is the different real part of the UPs giving rise to the dynamics, while the subexponential decay comes from the BCD contribution.
• Between these two regions, this bath induces a dynamical effect different from other reservoirs. As shown in Fig. 4(c), a QE at ∆ = 0 experiences first an exponential decay, which in the perturbative regime happens at a rate Γ M (0) ≈ 0.89g 2 /J. However, after this initial relaxation the QE displays coherent oscillations with a very slow decay. These oscillations can be attributed to the interference between the two BCD contribution and they are only weakly attenuated by a power-law decay. Since we have the analytical formulas for Σ e,CS (z) and ξ (z) in Eqs. 27 and 26, respectively, we can expand the integrand of Eq. 10 around z = ±2J − iy for |y| 1, which gives the asymptotic limit for t → ∞ of the BCD contribution. Since ξ (±2J − iy) ∼ a + b √ y, with a, b ∈ C, this propagates to the integral to arrive to C BCD (t) ∝ dy √ ye −yt ∝ 1/t 3/2 . Thus, the non-Markovian coherent oscillations decay with |C e (t)| 2 ∝ 1/t 3 in the asymptotic limit.
To conclude the study of this bath, it is illustrative to characterize how the excitation from a single QE gets emitted into the bath, especially around the plateau between [−2J, 2J], which is where the dynamics substantially deviates from the isotropic models. This is what we show in Fig. 5, where we plot the bath population in real space, |C R | 2 (in primitive coordinates), in a finite bath of size N = 2 7 after the QE has decayed for a time T J = N. We observe that in spite of having similar decay rates (because of the plateau in the density of states), the spatial decay at E = −2J differs significantly from the one at the middle of the band (E = 0). In particular, we observe that at the kinks of the density of states, E = ±2J, the excitation gets emitted in stripes around the QE, whereas at E = 0 it is emitted in 4 directions, namely, at (1, 1, 1), (1, 1, −1), (1, −1, 1) and (−1, 1, 1).
Both types of emission at E = 2J and 0 give rise to highly anisotropic collective decays when more QEs are coupled to the lattice. However, there seems not to be a perfect subradiant state with many QEs. Since the other lattice geometries display more interesting collective phenomena, we prefer to stop the discussion of CS bath here, and move to the next considered reservoir. The next bath geometry that we consider is the BCC lattice one, whose density of states (see Fig. 3(f))) differs significantly from the CS one since it displays a divergence in the middle of the band.

A. Single QE
In the continuum limit, the single QE self-energy reads: where ω BCC (k) is given in Eq. 20. Notice, it is not yet justified to change k i → −k i to restrict the range of integration to [0, π] because of the term cos(k 1 + k 2 + k 3 ). It is convenient to make a change of variables: In these new variables the dispersion relation factorizes: in the integral, the new integral values run at most from [−π, π], but actually the integration region is a parallelepiped, R, which goes from [−π, −π, −π] to [π, π, π]. The self-energy then reads: where the factor 4 is coming from the Jacobian of the transformation, i.e., |Det(U BCC )|. One can check that one can transport the parallelepiped into two cubes in [0, π] × [0, π] × [0, π] (and [−π, 0] × [−π, 0] × [−π, 0]), arriving then to the following expression for the self-energy: .
This integral can again be expressed analytically in terms of elliptic integrals [55] as follows: Evaluating the self-energy slightly above the real axis, i.e, Σ e,BCC (E +i0 + ), one recovers the expected perturbative decay rate, Γ e (E), and Lamb-shift, δ ω e (E), plotted in Fig. 3(b). The behaviour of Σ e,BCC (z) in the middle of the band resembles the one appearing in 2D Van Hove singularities [17,18,58]. However, in this case the δ ω e (E), apart from having discontinuous jump at E = 0, it diverges.
To be able to separate the different contributions to the dynamics of C e (t) using Residue Theorem one possibility consists in closing the contour of integration taking three detours at E = −8J, 0, 8J to avoid the non-analyticities of Σ e,BCC (z). This defines four different regions in the lower complex plane depending on Re(z) ∈ (−∞, −8J), (−8J, 0), (0, 8J), (8J, ∞), that we label from I-IV respectively. In regions II-III, the 1 − 64J 2 z 2 becomes complex, such that in order to go to the other Riemann sheet one must change: √ → − √ . The definition of K(m) will be always the one of Eq. 29 since m(z) never crosses the branch cut in all the integration region.
With these prescriptions, we separate the different contri-butions to C e (0) as a function of ∆ for a fixed g/J, as shown in Fig. 6(a). There, we observe that the band-edge effects are similar to the ones obtained in the CS lattice, with a BS contribution that vanishes at a given critical ∆. In the middle of the band, the QE dynamics is mostly governed by the UP, whose lifetimes are also renormalized with respect to the Markov predictions (see Fig. 6(b)). Like in the CS lattice, the divergence leads to the coexistence of the UPs due to the divergence of Σ e (E + i0 + ) around E ≈ 0. Moreover, it gives a non-negligible BCD weight at the middle of the band. Since now δ ω e (E) experiences a discontinuous jump from negative to positive around 0, the coexistence region is now symmetric around the divergence since both the UP in regions II/III can be obtained for ∆ > / < 0 respectively. To obtain the critical ∆ that characterizes the coexistence region of the two UPs, one must check when: z−∆−Σ e,BCC (z) = 0 has still two solutions (using the analytical continuation of the self-energy in the corresponding regions). Expanding Σ e,BCC (x − iy), for x, y J, we find: The complicated shape of the expansion of Σ e,BCC (z) prevented us from finding an analytical solution for the coexistence region of the two BS. Instead, we numerically solved the pole equations for ∆ = 0, and study the non-perturbative scaling of the real/imaginary part of the UP exactly at the middle of the band (shown in Fig. 6(c)). Moreover, in the figure we also give approximated analytical expressions that we found to fit the numerical data in the regions of g/J considered.
After having spotted the most interesting spectral region of this reservoir (around ∆ = 0), let us now study the dynamics at this point. In Fig. 7(a), we plot |C e (t)| 2 for a QE coupled to a BCC bath with g = J and ∆ = 0. One observes two features: first, we observe oscillations in |C e (t)| 2 , emerging from the interference of the two UP contributions. These oscillations are overdamped because the imaginary part of the UP is (slightly) larger than their real ones. Secondly, we observe a slow non-exponential relaxation for long times, whose origin is the extra BCD contribution appearing at ∆ = 0. In principle, one can try to obtain this scaling analytically using Eq. 10, which reads: G BCDII (y) = 2ReΣ e,BCC (0 + − iy) (y + ImΣ e,BCC (0 + − iy)) 2 + ReΣ e,BCC (0 + − iy) 2 .
Since the long-time limit is dominated by the behaviour at y J, we can use Eqs. 37 to obtain G BCDII (y J), and then use it to estimate analytically the asymptotic scaling of C BCDII (t). As we also show numerically in Fig. 7(b), G BCDII (y J) ∼ ln(y) −3 , however, the convergence to this limit is very slow, and for intermediate y's one observes G BCDII (y J) ∼ ln(y) −2 . In both cases, one can use the asymptotic expansion of Fourier integrals in Ref. [59] to obtain: For example, for the ranges of times considered in Fig. 7(a) we are probing only the ln(y) 2 limit.
Last, let us plot in Fig. 7(c) how a QE spectrally tuned to ∆ = 0 emits into the bath. There, we observe how the QE excitation spreads into the bath in a very directional fashion, like it also occurred for CS lattices. However, there is an important difference in this case, namely, the emission occurs only in three directions (instead of four), (1, 1, 0), (1, 0, 1) and (0, 1, 1). As we will see next, this will allow us to find perfect subradiant states when many QEs interact with the bath.

B. Many QEs
Let us now study collective phenomena when many QEs interact at the same time with a common BCC-type reservoir. Based on the knowledge of the single QE case, we focus on the case where the QEs are spectrally tuned to the middle of the band, ∆ = 0, since this is the region where they display a highly anisotropic emission into the bath, which translates into unconventional collective decay terms. Since the goal of the manuscript is to present the most relevant features of each bath geometry, we focus on the emergence of 3D perfect subradiance, something that was thought only to be possible when the QEs are at volume of the order of the wavelength of the QE frequency [60] or coupled equally to a cavity mode.
To have perfect subradiance, one needs to be able to cancel perfectly the emission of all QEs at the same time. Since the emission occurs in 3 directions, it can be shown that the minimal configuration to do so is to have 8 QEs. The QEs must be placed at the positions indicated in Fig. 8(a), namely, (±2n, 0, 0), (0, ±2n, 0), (0, 0, ±2n) and (±2n, ±2n, ±2n), where the parameter n controls the distance between QEs. To prove that a perfect subradiant state emerges in this configuration it is convenient to write H int of Eq. 4 in a collective basis of the QEs. For example, one can use the following rotation between the bare QE basis and a new one: Using this transformation, it is possible to show that the first collective state defined by this transformation, that we denote with annihilation (creation) operators σ sb ( † ), couples only to a bath mode defined by: where: U j → (q 1 , q 2 , q 3 ), (q 1 , q 2 , −q 3 ), (q 1 , −q 2 , q 3 ), . . . . One can show that A 8 is orthogonal to the bath modes coupled to the other collective QE modes. Thus, if we take as initial state |Ψ(0) = σ † sb |g ⊗8 , the only relevant part of the interaction Hamiltonian reads: Using this trick, the complexity of the problem reduces to the one of a single QE, but with a collective self-energy which is given by: where the extra factor 2 comes from transporting the parallelepiped into two cubes of size: [0, π] × [0, π] × [0, π], like it also happened for the single QE self-energy. Exploiting the symmetries of the integrand, it is trivial to check that Σ sb (0) = 0, which means that there is indeed a real pole in the middle of the band at E = 0, which therefore does not decay. The only missing thing to prove is that its associated residue is different from zero. Remarkably, since the integrand is separable at z = 0, this residue can still be calculated analytically, yielding: which is therefore finite, and close to 1 if 2g 2 n 3 J 2 1. In Fig. 8(b-d) we certify numerically all these analytical predictions. To start with, in Fig. 8(b) we plot the dynamics of |C sb (t)| 2 for a situation with g = 0.1J and several distances parametrized by the number n (see Fig. 8(a)) obtained by simulating the full QE-bath Hamiltonian in a system with N = 2 7 sites. We observe both the perfect subradiant effect for small distances, and the correction imposed by retardation at large ones. This retardation effect appears due to the finite propagation speed of the excitations in the bath, which makes that the interference can only take place after a certain time which depends on the distance between QEs [18]. To further explore this retardation effect, we plot |C sb (t → ∞)| 2 for several distances and ratios g/J, comparing both the numerical simulations of the full QE-bath Hamiltonian (markers) and the analytical prediction of Eq. 46, showing very good agreement between them. Finally, in Fig. 8(d), we plot the bath population at a time tJ = N for g = 0.1J and n = 5, showing how the bath emission gets trapped between the 8 QEs due to destructive interference between the different decaying paths.

VII. QUANTUM DYNAMICS IN FACE-CENTERED CUBIC BATHS: ANISOTROPIC DIPOLE-DIPOLE INTERACTIONS
In this section, we consider a bath composed by an FCC lattice, depicted in Fig. 3(c), whose density of states also present qualitative differences with respect to the CS and BCC ones.

A. Single QE
The single QE self-energy of the FCC lattice in the continuum limit: where ω FCC (k) was written in Eq. 21 in terms of the primitive reciprocal coordinates, k = (k 1 , k 2 , k 3 ). Like we did for the BCC lattice, it is convenient to make a change of variables: under which the dispersion relation transforms to: The integration region of the self-energy in the new variables, i.e., q 1 = (k 1 + k 2 − k 3 )/2,. . . , is a polyhedron where each q variable runs at most from −3π/2 < q i < 3π/2. Using that the Jacobian of the transformation |DetU FCC | = 2, and transporting the polyhedron into the cube [0, π] × [0, π] × [0, π], one arrives at Using this transformed expression, it is possible to obtain an analytical solution of the self-energy [55] which reads [61]: Using this analytical expression we can obtain both the imaginary, Γ e (E), and real part, δ ω e (E), of the self-energy above the real axis, as depicted in Fig. 3(g). We observe three spectral regions where Σ e,FCC (z) shows a non-analytical behaviour, namely, E = −12J, 0, 4J, and thus they will force us to take a detour in the integration contour to obtain C e (t). This divides the lower half plane in IV regions, where the definition of Σ e,FCC (z) used must be slightly modified depending on the Riemann sheet we want to move in each one: • Outside of the band, that is, when Re(z) < −12J and Re(z) > 4J, that we label as regions I and IV respectively, one can use Σ e,FCC as defined in Eq. 53.
• Inside of the band one needs to distinguish two regions, from −12J < Re(z) < 0 and 0 < Re(z) < 4J, that we denote as regions II and III respectively. In these regions one of the square roots of ξ (z) in Eq. 53 becomes complex, and thus, one needs to adapt its definition to go to the adequate Riemann sheet in each integration region. In particular: With these prescriptions, we can now separate the different contributions to the dynamics of C e (t) depending on the spectral region, ∆, the QE is coupled to. This is what we summarize in Fig. 9, where we plot the (absolute value) of the weight of the BS, BCD and UP contributions for a QE coupled with g = J. To avoid unnecessary lengthening of the manuscript, we just enumerate common features already appearing in other reservoirs (e.g., CS and BCC) and make an extended discussion on the qualitative new effect emerging in FCC structures: • The lower edge at E = −12J shows the typical behaviour of 3D isotropic baths, with a sudden disappearance of the BS contribution.
• At E = 0, the density of states is finite, but with a dis-continuous derivative. This leads to the coexistence of two UP contributions and an extra BCD contribution in the middle of the band, resulting in non-exponential slow relaxations.
• The differential feature of this reservoir occurs at its upper band-edge, E = 4J, where both the real and imaginary part of Σ e,FCC (z) diverge, something very atypical for 3D reservoirs. This divergence leads to two remarkable consequences: i) enhancement of its decay rates as they get closer to the band edge and, ii) the survival of the BS for all the spectral region [15].
To learn more about this behaviour, we can expand the single QE self-energy around E ≈ 4J + x, obtaining where Θ(x) is the Heaviside function. We have checked numerically (not shown), that these formulas reproduce relatively well the behaviour even for values of |x| J. Using these formulas for x > 0, we can obtain an approximated solution for the energy of the upper BS (UBS) for ∆ = 4J, which reads: where W (x) is the so-called Lambert or product-log function, which satisfies W (x) ≈ ln(1/x) for x 1 and W (x) ≈ x for x 1. We have checked (not shown) that this formula reproduces well the behaviour obtained from numerically solving the pole equation at ∆ = 4J for a wide range of g/J.
Apart from its energy, it is also illustrative to study how the wavefunction of this BS looks like in real space since it will be directly connected to the dipole-dipole interactions when more than one QE is interacting with the bath. The wavefunction can be shown to be given by the following integral: where n = (n 1 , n 2 , n 3 ) are the primitive coordinates. In Figs. 9(b-c), we plot the the numerical evaluation of this wavefunction for a finite system with N = 128 sites and several detunings with respect to the band-edge: δ = E UBS − 4J. We highlight two characteristics: • The BS is highly anisotropic spreading mainly in the directions (n, 0, 0), (0, n, 0), (0, 0, n).
• To explore the decay along these axis, we note that when δ J, the most relevant k-points in the integral of Eq. 57 are the ones that satisfy ω(k) = 4J, since they are energetically closer. It can be shown that this equation defines 12 lines in the three-dimensional k-space, such as k 1 = π and k 3 − k 2 = ±π, making a closed loop in the integration region. If we are interested in one particular direction, e.g., (n, 0, 0), the most relevant contribution will come from points where k 1 ≈ 0, e.g., (0, π, π), since the others will lead to a highly oscillatory integrand. Expanding the dispersion relation around those points we find, e.g., for q i 1. Using this expansion, we can make the q 1 integral by extending the limits to (−π, π) → (−∞, ∞) and applying Residue Theorem to arrive at: which gives us the dependence with δ . The dependence with n is embedded in a function F(n), which unfortunately cannot be calculated using the expansion of Eq. 58 since gets a divergent integral. Through numerical inspection, we find F(n) can be well fitted by a power-law decay F(n) = A/n β , with A ≈ 0.02 and β ≈ 0.5, as shown in solid lines in Fig. 9(c). These features are very different from the models considered in the literature [40,41], where the BSs were always isotropic and with a decay law of the Yukawa type, ∝ e −r/ξ /(r/ξ ), with ξ being the bound state length and r the distance to the impurity, also inversely proportional of the detuning with respect to the band-edge.

B. Many QEs.
Using the intuition from other reservoirs [24,25,57], whenever a BS emerges around a QE, it can be used to mediate coherent interactions when many QEs are present in the structure. These interactions will naturally inherit the characteristics of the BS, which in the case of the FCC lattice we have just shown are very distinct from the known behaviours, giving rise to very exotic spin Hamiltonians.
To prove the emergence of these purely coherent dipoledipole interactions, we consider a pair of QEs spectrally tuned close to the upper edge, that is, ∆ 4J. Their initial state is assumed to be |Ψ(0) = |e 1 g 2 , with one excited and one in the ground state. In Figs. 10(a-b) we show the result of numerical evolution of the complete QE-bath Hamiltonian for two QEs coupled with g = 0.1J to an FCC lattice and ∆/J = 4.1 and 4.3, respectively. Their dynamics show indeed that the excitation transfers coherently from one QE to the other with no exponential attenuation. This is what one expects for the Markovian regime from Eqs. 16 and 17 in Section III, which for this particular case lead to: where: is nothing more than the discrete version of the integral of Eq. 57 which gives C n . This shows that the dipole-dipole interactions will inherit the exotic properties of the upper BS of the FCC lattice. Moreover, we also notice that in order to accurately capture the frequency of the oscillation for the parameters of Figs. 10(a-b), we need to solve the pole equation in the symmetric/antisymmetric subspace exactly, that is, finding: In the solid lines of Figs.10(a-b), we plot the solutions of Eq. 60 with a renormalized frequency J exact = (J − −J + )/2 in solid lines showing very good agreement with the results from a simulation of the full QE-bath dynamics. For completeness, in Fig. 10(c), we compare J exact (markers) versus J M (lines) for several distances, showing how the Markov approximation becomes more accurate the smaller the distance (or the smaller g/J, not shown).

VIII. QUANTUM DYNAMICS IN DIAMOND BATHS: ANISOTROPIC DIPOLE-DIPOLE INTERACTIONS
In this Section, we explore one example of a 3D two-band model: the diamond lattice. As we explained in Section IV, this bath geometry is formed by two FCC lattices displaced by a vector 1/4(1, 1, 1), as depicted in Fig. 3(d). This bath geometry leads to the appearance of a singular band-gap in its density of states, shown in Fig. 3(h), which is not present in the other 3D reservoirs considered. Thus, this spectral region will be the main concern of this Section. • In the regions Re(z) ∈ (−4J, −2J) and (2J, 4J), denoted as II and V, one must replace:

Thermodynamic limit
• Finally, when Re(z) ∈ (−2J, 0) and (0, 2J), denoted as regions III and IV, one must change: 4 − 16J 2 z 2 → − 4 − 16J 2 z 2 (as well as 1 − 16J 2 z 2 → − 1 − 16J 2 z 2 ). With these prescriptions, we can separate the contributions of C e (0) of the BSs, BCDs and UPs in the different spectral regions like we did for other reservoirs. This is what we show in Fig. 11(a) for a QE coupled to a diamond bath with g = J/2. We observe similar phenomena to the other reservoirs at the band edges, ∆ ≈ ±4J, with a sudden disappearance of the BS, and at the kinks, ∆ ≈ ±2J, with coexistence of UPs contributions. The main difference with respect to the other reservoirs occurs at ∆ ≈ 0, where another BCD must be taken at E = 0, that gives an extra BCD contribution, denoted as BCDIII and plotted in orange stars in Fig. 11(a) [62]. To further understand the behaviour around this point, it is enlightening to expand Σ e,diam (E + i0 + ) around it: for |E| J. Since Σ e,diam (0) = 0, E = 0 is a solution of the pole equation when ∆ = 0. However, its derivative ∂ z Σ e,diam (z)| z=0 → ∞ has a logarithmic divergence, and therefore its associated residue in the thermodynamic limit must be zero. This is why as we take ∆ closer to the singular point, the BCDIII contribution becomes more important in Fig. 11(a). To certify that, we plot in Fig. 11(b) the associated dynamics for the parameters of panel (a), and several detunings approaching ∆ → 0. We observe that the dynamics is given first by an exponential decay, given by the UP contribution ,followed by a subexponential relaxation. Since the imaginary part of the UP also goes to 0, as predicted by Fermi's Golden Rule and the expansion of Eqs. 67, the dynamics becomes slower, and eventually becomes fully dominated by the BCD contribution. A similar behaviour occurred for 2D singular band-gaps (Dirac points) [57], where it was predicted a 1/ ln(t) decay for a QE in the thermodynamic limit. For finite ∆, but still close to 0, the C BCDIII (t) ∝ 1/t 2 , which is the power law that we observe in Fig. 11(b) for, e.g., ∆ = 0.1J.
After having spotted similarities with 2D Dirac points, it is instructive to consider the role of finite size effects in the diamond bath, since they were proven to play an important role in the physics emerging at these points at the 2D case. In particular, in 2D it was shown [57] how for finite systems (with periodic boundary conditions) the spontaneous emission gets quenched because of the emergence of a quasibound photonic state around the impurity. In Fig. 11(c), we plot a comparison of |C e (t)| 2 for a QE coupled with g = 0.1J and ∆/J = 0.001 calculated in the thermodynamic limit using resolvent operator techniques (solid line) and using a finite-bath simulation with periodic boundary conditions (marker). We observe how both calculations agree perfectly well until times tJ ≈ N, being N the linear size of the system. However, in contrast to what happens in 2D singular band-gaps, here the dynamics accelerates instead of quenching. Moreover, we observe that this acceleration occurs for times proportional to N.
To further investigate this behaviour we make a longer time simulation of |C e (t)| 2 for a QE coupled to the A lattice, that we plot in black squares in Fig. 12(a-b), together with the dynamics of the total A/B bath population in blue spheres/red triangles, respectively. We observe two different behaviours depending on the timescales: -For short times, the QE relaxes very slowly (as we already show in Fig. 11(c)) decaying mainly into the B sites, which is more clear in the zoom we make in Fig. 12(b). Moreover, inspecting the distribution of these modes in real space, shown in Fig. 12(c), we realize the decay into the B bath has a very anisotropic shape resembling the one of the BS of the FCC lattice at the upper edge. This is expected since the k modes resonant at the QE frequency, ∆ = 0, are the same ones that give rise to the divergence in the upper edge of the FCC lattice. One remarkable difference with respect to the FCC BS is that the spatial decay of the wavefunction, |C B,n |, along the main axis (see, e.g., Fig. 12(d)), does not show an exponential attenuation, but rather seems to follow a power law decay ∼ 1/n, being n the distance from the QE. This behaviour at short times resembles the one of 2D Dirac cone, where a qua-siBS emerges in the B/A lattice for a QE coupled to the A/B lattice site with a power-law localization of its wavefunction.
-For long-times, on the contrary, the dynamics differs significantly from the 2D situation, where the QE decay freeze around a constant value [57]. Here, instead, the excitation from the QE gets completely transferred to the A bath after a certain time. This can be understood from the existence of a collective mode of the A 0 -bath at zero energy, which is able to resonantly transfer excitation from the QE to the bath. This mode can be defined as: where N A 0 is the number of k-modes which satisfy | f (k)| 2 = 0, which can be shown be given by the same closed contour composed of 12 lines that we explained to describe the upper BS of the FCC lattice. Thus, this results into a resonant interaction between the QE and the collective mode of the type: g A (σ eg A 0 + h.c.), which gives rise to Rabi oscillations between the QE and A 0 , with frequency: which reproduces the behaviour of Fig. 12(a). The spatial distribution of this mode is plotted in Fig. 12(e), for a time T J = 9N, where most of the population has been transferred to the A bath. Since the g A ∝ 1/N, in the thermodynamic limit this contribution vanishes and eventually only the relaxation from the B bath appears.
To conclude this section, let us remind that all the results with finite baths are obtained using periodic boundary conditions. Although some of the results may vary with different boundary conditions, e.g., the emergence of the zero energy mode, for the sake of concreteness we leave the detailed discussions with other boundary conditions for further works.

B. Many QEs
As a final illustration, let us consider how the non-trivial dynamics emerging at ∆ = 0 translates to the situation where many QEs are interacting with the diamond bath. In particular, we consider a situation where two QEs, 1 and 2, are coupled respectively to the A/B lattice. Then, we set the first QE in the initial state, and study whether the excitations get coherently transferred to the second one. In Fig. 13 we show the dynamics of the two QEs, |C 1,2 (t)| 2 in solid black/green, and of the total bath population |C A,B (t)| 2 in the A/B lattices in blue/red, respectively. In analogy to what happened in the single QE situation, the behaviour occurs in two different timescales: the initially excited QE starts oscillating back and forth with the collective bath mode it is coupled to. However, as time passes there is a small fraction of population being transferred to the second QE (and its corresponding bath), until it gets transferred completely.
From here, there are many research directions to continue exploring, such as what happens in the limit N → ∞, when the coupling to the collective zero-energy mode vanishes, or when QEs couple to the same sublattice, that we leave open for further works.

IX. BATH IMPLEMENTATION WITH OPTICAL LATTICES
In Section II B we introduced the general ideas on how to simulate quantum optical phenomena with cold atoms in statedependent optical lattices based on Ref. [40]. In the original proposal, however, the bath was considered to be freeparticles with energy dispersion ω(k) ∝ |k| 2 . To observe the phenomenology explored in Sections V-VIII, the key ingredient is to be able to engineer more complicated bath energy dispersions, ω(k), like the ones considered along the manuscript. Some of the bath geometries, like CS lattices, are straightforward to generate using three retro-reflected laser beams and they are used nowadays in most of the 3D state-of-art experiments [51]. The other geometries, however, have received much less attention [63][64][65] such that it is still worth revisiting them and characterize their band structure.
In this Section: i) we provide the laser configurations to obtain the lattices explored along this manuscript; ii) we calculate the associated band structure of the laser configurations considered; iii) we analyze to which extent the dynamics predicted in this manuscript could be observed in these setups by studying the associated energy scales and potential problems, such as the emergence of longer-range hopping, which may alter the density of states of the ideal models. We focus on the study of the density of states since its non-analytical behaviour is crutial in most of the phenomenology predicted in the manuscript.
We also acknowledge there will be other experimental imperfections that may affect the predicted dynamics, e.g., decoherence. These effects have been studied for lower dimensional baths, e.g., in Refs. [53,57], where the rule of the thumb to observe the phenomena was that it occurs in timescales faster than the decoherence ones, something that is within the reach with these cold atoms setups. Since the goal of the manuscript is to uncover new phenomena, we leave for future works the systematic characterization of other experimental restrictions.

A. Tools & Analysis
Atomic optical potentials. The interaction of one or many lasers, with a given wavelength λ , with an atomic optical transition generates optical potentials which can be generally writ- ten as [66]: where V 0 is the overall amplitude of the optical potential, and I(R) its intensity profile. The parameter V 0 can be controlled in magnitude and sign through both the laser intensities and detunings with respect to the optical transitions. Along this manuscript, we will assume V 0 to have units of energy, such that I(R) is dimensionless. The intensity profile, I(R), depends on both the atomic polarizability tensor and its interplay with the total electric field. For simplicity, we assume to work in the regime where the polarizability tensor is isotropic and the intensity profile is just given by the total electric field profile, |E(R)| 2 resulting from the interference of the different lasers interacting with the atom. The first task in the following Sections will be to find laser configurations that give rise to optical potentials with the same Bravais structure than the bath geometries depicted in Figs. 3(a-d). For that purposes, we use two experimental resources: i) On the one hand, we exploit the interference between several laser fields with the same wavelength (or equivalently frequency). Each laser field is characterized by its amplitude/polarization/propagation vector {E i ,ε i , p i }, that is, they have associated a vector electric field : E i (r) = e ip i ·r E iεi . Since we assume they all have the same frequency, the interference leads to an electric field intensity which reads: It contains both a constant shift of the potential coming from the self-interference of each electric field, and a position dependent contribution coming from the interference of each pair of electric fields.
ii) Since sometimes these patterns will not be enough, we will need to add up the contribution of other laser fields with different wavelength (and frequencies, ω α ) such that they result in a time-independent averaged potential: where E α (R) is the total electric coming from the interference of the laser fields with the same wavelength ω α as given by Eq. 71. The cross-interference between the lasers with different frequencies averages out in our timescales of interest when ∆ω αβ = ω α − ω β is bigger than any parameter of our simulated Hamiltonian ( 10 kHZ), which is the regime we will assume to be working in. Since: for λ i ∼ 500 nm, we have that ∆ω ≈ 2π × 10 9 ∆λ KHz/nm, such that they can be easily satisfied even for very closely spaced wavelengths, which we consider to the same for the purposes of the defining the lattice geometry.
Characterizing bath structure Once we find an appropriate laser configuration for each bath, we will characterize the emergent atomic dynamics for atoms hopping in a given potential V (R). This is important since having the same periodicity as the original model does not guarantee that the dynamics will be described by the same energy dispersion. In particular, longer range hoppings may emerge that deviate the dynamics from the nearest neighbours descriptions that we considered in the previous Sections.
The Hamiltonian describing the dynamics of atoms hopping in these optical lattices can be generally written [51]: where M is the atomic mass and ∇ 2 = ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 + ∂ 2 ∂ z 2 is the Laplacian in cartesian coordinates. To make the problem adimensional, it is convenient to use λ = λ /(2π) and the recoil energy E R =¯h 2 2Mλ 2 as the natural units of length and energy respectively [67]. With these units, the Hamiltonian of Eq. 74 is reexpressed:Ĥ where for simplicity we have kept the same notation for H,V (R) and V 0 , even though they are all now written in units of the recoil energy. Since we are interested in calculating the lowest energy part of the spectrum, we assume to have translational invariant system to apply Bloch Theorem [68]. Under this assumption, k, is a good quantum number such that the eigenfunctions can always be written as: |Ψ(R) = e ik·R |u k (R) , where |u k (R) are the so-called Bloch-modes. Using this ansatz, the eigenvalue equation transforms to: where the index n denotes the different bands, E k,n , appearing in the model. The k = ∑ 3 j=1 k j d j , where d j are the reciprocal primitive vectors of each lattice used to define the periodic boundary conditions of the problem. Since the reciprocal primitive vectors are different for each lattice considered, the shape ofĤ k varies for every lattice geometry, like we will show explicitly in the next Sections. Finally, since both the potential V (R) and the Bloch mode, |u k,n (R) , have the same spatial periodicity, it is convenient to expand them in terms of a finite number of reciprocal wavevectors: (77) where q max is the numerical cut-off to expand the Bloch modes in the reciprocal space. Thus, to obtain E k,n for each k-point, one needs to solve a eigenvalue equation for matrices of size (2q max + 1) 3 . This problem simplifies a lot for the case of separable potentials, i.e., V (R) = V (x) +V (y) +V (z), where one can solve directly the 1D problem where the size of the the matrix is just (2q max + 1). Except for the diamond lattice case, that will be treated separately in Section IX E, we focus only on the lowest energy band such that from now on we drop the index n of the discussion. Like we mentioned in the beginning of this Section, we will calculate the following figures of merit to estimate to which extent the physics predicted in this manuscript can be observed within optical lattices: • On the one hand, we calculate the strength of the atom (nearest neighbour) hoppings. It can be easily shown that the hopping rate between two atoms separated a distance n can be obtained from by Fourier transforming E k : where N is the (linear) number of sites that the use to discretize the k i variables.
• On the other hand, we also calculate the numerical density of states of the bath, by discretizing the energy space, which can be calculated from E k , which can be calculated as follows: where Θ(x) is the Heaviside function, ω n = minE k + n maxE k −minE k N ω , being N ω the number of steps in which we discretize the frequency space and n = 1, . . . N ω . We will be especially interested to check whether the non-analytical features observed in Figs. 3(e-h) survive in the real lattices.
After we have given the tools and prescriptions to analyze the problems, let us study the results obtained for the different lattice geometries considered.

B. Cubic Simple lattices
This lattice is simple to generate [51]: one just needs three retro-reflected lasers, with orthogonal polarizations, propagating in the x,y,z directions, that is, E i = E 0 , p 1,2/3,4/5,6 = ± 2π λê x/y/z andε 1,2/3,4/5,6 =ê z/x/y . The interference of these six electric fields results in a potential: Since the potential and the kinetic energy terms are both separable in this case, one can treat the problem of each dimension separately. This allows us to make calculations for very large lattices with little numerical effort.
The first thing we are interested in is the scaling of the n-th neighbour hopping rate, that we denote as J n , with the lattice potential V 0 . In this geometry there are six nearest neighbours at positions (±1, 0, 0), (0, ±1, 0), (0, 0, ±1). Due to the separability of the potential, the longer range hoppings also appear only in these directions, such that we can denote as J n to the hopping rate of the n-th neighbour at position, e.g., (n, 0, 0). In Fig. 14(a) we plot the result of applying Eq. 79 with these positions for a system with N = 200 lattice sites and four different values of |V 0 /E R | = 4, 8, 12, 16. As expected, the larger V 0 , the smaller J 1 since the Wannier functions in each site get more localized and their overlap decrease. The longer range hoppings decay quickly with the distance due to the exponential localization of the Wannier wavefunctions.
Finally, in Fig. 14(b) we plot in solid lines the numerical density of states for the lattice depths considered in the panels (a-b), together with the expected one from the nearest neighbour model (in dashed black lines). For the shallowest lattice considered, in blue, we observe how an asymmetry between the two kinks in the middle of the band appear. However, as the potential depth increases one quickly recovers the nearest-neighbour density of states. Remarkably, the nonanalyticities giving rise to non-Markovian phenomena seem to survive even in the cases where the model deviates from the ideal nearest neighbour description.
To calculate the band structure of this potential it is convenient to express the position/momenta in the eigenvalue equation in terms of the primitive vectors of real/reciprocal space. For example, the Laplacian, ∇ 2 , written in terms of n = ∑ j n j c j is given by: (83) The potential in primitive coordinates is simply obtained by: V BCC (x, y, z) = V BCC (π(n 1 + n 2 − n 3 ), π(n 1 − n 2 + n 3 ), π(−n 1 + n 2 + n 3 )). With this change of variables, one can calculate E k using the prescriptions explained in Section IX A. Since the potential is not separable and the calculations are more demanding we calculate the band structure for smaller system sizes than the CS lattice, and then extrapolate the results to larger lattices to obtain a smooth density of states.
We numerically calculated E k using the potential V BCC (R) and realize that it deviates from the energy the energy dispersion of the nearest neighbour model (not shown). The underlying reason is that the potential in the direction of nextnearest neighbour is shallower than in the nearest neighbour direction, such that it still has a big effect in spite of the larger distance. To solve this problem we propose to use another set of 3 retro-reflected lasers with different frequency but virtually indistinguishable wavelength, as explained in Section IX A, and with negative V 0 , such that they cancel the V CS (R) contribution of V BCC (R). In that case the potential finally reads: Another possibility to obtain V BCC,2 (R) consists in summing three different sets of lasers, with the same amplitude E 0 and slightly different frequencies, with the following propagation/polarization vectors: The first set is composed by p 1/3 = π λ (1, ±1, 0) = −p 2/4 , with polarizationŝ ε 1,2/3,4 = 1 √ 2 (1, ∓1, 0). The second set is composed by p 5/7 = π λ (1, 0, ±1) = −p 6/8 , with polarizationsε 5,6/7,8 = 1 2 (1, 0, ∓1). The third set is finally composed by p 9/11 = π λ (0, 1, ±1) = −p 10/12 , with polarizationsε 9,10/11,12 = 1 √ 2 (0, 1, ∓1). The time-averaged potential resulting from the sum of the interference of each set leads to: + cos(x − z) + cos(y + z) + cos(y − z) = = cos(x) cos(y) + cos(x) cos(z) + cos(z) cos(y) .
Irrespective of the method used to obtain V BCC,2 (R), we calculate the lowest energy band of the atoms hopping in this potential and summarize the results in Figs. 15(b-c). In Fig. 15(b) we plot the scaling of the n-th nearest neighbour hopping for several lattice depths V 0 given in the legend, where we observe that indeed |J 1 | |J n>1 | for large lattice depths. Finally, in Fig. 15(c) we plot the numerical density of states for the lattice depths considered in the previous panels. For shallow lattices, the longer range hoppings renormalize the divergence of the density of states in the middle of the band. For the deeper lattices, the results converge to the ones expected from the nearest neighbour description.
Like we did for the BCC potential, to calculate the band structure it is convenient to express the position/momenta in the eigenvalue equation in terms of the primitive vectors. For example, the Laplacian, ∇ 2 , in terms of n = ∑ j n j c j is given by: The potential written in primitive coordinates can be obtained replacing: x = π(n 1 + n 2 ), y = π(n 1 + n 3 ), z = π(n 2 + n 3 ). With this change of variables, we calculate E k using the prescriptions explained in Section IX A. The results are summarized in Fig. 16, where we plot both the scaling of the n-th neighbour hopping for several V 0 /E R in Fig. 16(b), and the extrapolated numerical density of states in Fig. 16(c). We observe that the numerical density of states approaches the nearest neighbour one like we desired.

E. Diamond lattices
We already explained in Section IV that the diamond lattice is formed by two interspersed FCC lattices displaced by a vector f = 1 4 , 1 4 , 1 4 . Using this intuition, one could build an optical potential by doubling the configuration of the FCC  Table: number of n-th neighbours, distance and an example position of one of the n-th neighbour.
(b) Scaling of the dominant nearest neighbour contribution, J AB 1 , for several lattice depth potentials V 0 /E R . (c) Bath density of states for the V 0 /E R considered in panel (b). In dashed black lines we plot the expected result for the nearest neighbour model in the thermodynamic limit. Numerical simulations done for a bath of linear size 2 5 and momentum cut-off q max = 7. lattice to make a second displaced potential using the trick of making two slightly detuned sets of lasers. The final optical potential will be given by the sum of the two displaced potentials: The primitive vectors for the real/reciprocal space are those of the FCC potential, as well as the nearest neighbours within the same sublattice. Besides them, there are extra neighbours emerging from the coupling to the other sublattice sites, depicted and summarized in Fig. 17(a). Furthermore, another difference of this situation with re-spect to the case of a simple Bravais lattice is that we need to calculate not only the lowest, but also the first excited band, that we denote as E 1,k and E 2,k to characterize the behaviour of the system. In this type of models with a superlattice the Wannier function are not unambiguously defined [69], however, the maximally localized ones lead to the following formulas for the tunneling within the AA (BB) or AB lattice sites [70]: For the sake of concreteness, we will only plot here the scaling of the most relevant nearest neighbour coupling, J AB 1 , for several V 0 /E R in Fig. 17(b), and the corresponding density of states in Fig. 17(c), calculated extrapolating the results of a small system size simulation like we did in the previous cases. Differently from what happened in the other types of baths, with the potential configuration we propose there is no range of V 0 /E R that approximates D diam (ω) in all spectral regions. While the lowest band is reasonably well approximated for all the lattice depth range considered, the upper one deviates significantly for the idealized nearest neighbour model, being more similar around V 0 /E R ∼ 3. However, it is worth emphasizing that the regions around ω ≈ 0 does indeed show a singular bandgap, which points to the possibility of observing the most distinctive dynamics of this type of reservoirs.

X. CONCLUSIONS & OUTLOOK
Summing up, we have systematically characterized the quantum dynamics of QEs coupled to several 3D structured reservoirs with different qualitative features in their density of states. Through exact calculations, we predict the emergence of several phenomena beyond the traditionally considered band-edge related effects, such as: • Long-lived reversible dynamics for single QEs spectrally tuned within band frequencies in CS lattices.
• Directional emission and the emergence of perfect subradiant states in BCC lattices.
• Robust 3D anisotropic bound states which survive irrespective of the QE spectral detuning/coupling for FCC lattice. This is in stark contrast to standard 3D bound states appearing in isotropic band-edges, which for a fixed g/J merge into the scattering spectrum for a critical ∆.
• Through these robust bound states the FCC bath also mediates QE dipole-dipole interactions only in certain directions. The spatial decay of the interactions in these directions can be numerically fitted to: ∼ e −d √ δ n / √ δ n, very different from the isotropic Yukawa type interactions of standard 3D reservoirs.
• Sub-exponential relaxation of single QEs tuned at the singular bandgap of diamond lattices. Furthermore, when many QEs are interacting with the bath at this frequency one observes reversible exchange of excitations with frequency scaling as 1/n, being n the distance between emitters.
Furthermore, we also propose a way how to observe these effects with cold atoms in state-dependent optical lattices based on the proposal of Ref. [40]. More concretely, we i) provide the laser configurations to design optical potentials with the geometries of the bath considered, and ii) characterize their lowest energy band (or the two lowest in the case of diamond geometries). In particular, we study the associated timescales of the hopping models and to which extent their dynamics reproduce the physics of the nearest neighbour idealized models considered along the manuscript, taking their density of states as a figure of merit for the comparison. We conclude that the phenomena predicted along this manuscript is within the reach of this platform, that together with the recent experimental developments [42], foresee the observation of non-trivial dynamics in near future experiments.
The work constitutes a solid basis for future investigation of 3D structured quantum optical systems and opens several research directions. From the fundamental point of view, natural extensions of this work are the study of the emergent phenomenology in the many excitation regime [71], exploring the interplay between driving and the anisotropic collective dissipation appearing in these systems which may lead to the emergence of many-body entangled steady states [72,73], or the characterization of the phases emerging in the effective spin models. Another interesting direction consists in combining the ideas developed along this manuscript with the recent developments in confined photons in subwavelength atomic lattices [19][20][21][22]. From the more applied perspective, it will be interesting to find simpler laser configurations which give rise to the 3D models considered and make a more thorough characterization on the emergent dynamics with these potentials.