Light propagation in systems involving two-dimensional atomic lattices

TL;DR

This paper investigates how light interacts with two-dimensional atomic lattices, revealing collective effects, resonance shifts, and potential for emulating 1D waveguides, with implications for controlling light transmission in atomic systems.

Contribution

It provides a detailed analysis of light propagation in 2D atomic lattices, including effects of dipole interactions and potential to mimic 1D waveguides, extending understanding of light-matter interactions in such structures.

Findings

Lattice atoms exhibit multiple resonance frequencies due to dipole interactions.

Stacks of 2D lattices can emulate 1D waveguides with lossless light propagation.

Magnetic fields can induce transparency or opacity in the lattice at specific frequencies.

Abstract

We study the optical response of a 2D square lattice of atoms using classical electrodynamics. Due to dipole-dipole interactions, the lattice atoms polarize as if the lattice were an atom with up to three resonance frequencies, with cooperatively shifted resonances and altered transition linewidths. We show that when the distance between two 2D lattices is large enough and Bragg reflections are absent, the lattices interact among themselves as if they radiated a plane wave whose amplitude is in accordance with the radiation from a dipole moment continuously distributed in the lattice plane. We employ these results to study light propagation in stacks of 2D lattices, drawing on simple qualitative pictures of the response of a 2D lattice and light propagation in 1D waveguides. We show that a stack of 2D lattices may emulate regularly spaced atoms in a lossless 1D waveguide, and argue that in a suitable geometry the resonance shifts characteristic of 1D and 2D lattice structures may completely cancel to eliminate density dependent resonance shifts of atoms bound to a 3D lattice. A generalization to the case of anisotropic polarizability, such as in the presence of a magnetic field, reveals light frequencies induced by the magnetic field for which the lattice is either completely transparent, or completely opaque.

Equations

$$k=c=\hslash =\frac{1}{4\pi {ϵ}_0}=1.$$

$${\mathsf{G}}_{ij}(\mathbf{r}-{\mathbf{r}}^{\prime })=$$

$${\hat{\mathbf{e}}}_i\cdot \left\{(\hat{\mathbf{n}}\times {\hat{\mathbf{e}}}_j)\times \hat{\mathbf{n}}+[3\hat{\mathbf{n}}(\hat{\mathbf{n}}\cdot {\hat{\mathbf{e}}}_j)-{\hat{\mathbf{e}}}_j]\left(\frac{1}{r^2}-\frac{i}{r}\right)\right\}\frac{e^{ir}}{r}.$$

$$\delta =\frac{\omega -{\omega }_0}{\gamma }.$$

$$\alpha =-\frac{3}{2(\delta +i)}.$$

$$\mathbf{E}(0)=\sum _{\mathbf{n}\ne 0}\mathsf{G}({-\mathbf{R}}_{\mathbf{n}})\mathbf{d}.$$

$$E(0)$$

$$S(a)$$

$$S(a)=\underset{\eta =0+}{\lim }\left\{\underset{M\to \infty }{\lim }\sum _{\mathbf{n}\ne 0}^{|\mathbf{n}|\le M}\frac{e^{(i-\eta )\rho }\left[{\rho }^2-i\rho +1\right]}{2{\rho }^3}\right\}.$$

$$T(z,a)$$

$$=\sum _{\mathbf{n}}\frac{\left[r^4-i{r}^3+r^2\left(z^2+1\right)+3ir{z}^2-3z^2\right]e^{i(r-|z|)}}{2r^5};$$

$$r=\sqrt{(na{)}^2+z^2}.$$

$$\bar{T}(a)=\frac{2\pi i}{a^2},$$

$$E=E_0+S(a)d,d=\alpha E;$$

$$d=-\frac{3E_0}{2[\delta +\frac{3}{2}S(a)+i]}=-\frac{3E_0}{2(\Delta +i\Gamma )},$$

$$\Delta =\delta +\frac{3}{2}\Re [S(a)],\Gamma =1+\frac{3}{2}\Im [S(a)]$$

$$\bar{T}(a)=it(a),t(a)=\frac{2\pi }{a^2}\in \mathbb{R}.$$

$$E_R$$

$$E_T$$

$$\Gamma -\frac{3}{2}t(a)=0.$$

$$\Im [S(a)]=\frac{2\pi }{a^2}-\frac{2}{3}.$$

$$R,T={\left|\frac{E_{R,T}}{E_0}\right|}^2.$$

$$d=-\frac{D^2}{\Delta +i\Gamma }{E}_0,$$

$$E_R=itd=-\frac{it{D}^2}{\Delta +i\Gamma }{E}_0$$

$$E_T=E_0+E_R=\frac{\Delta +i(\Gamma -t{D}^2)}{\Delta +i\Gamma }{E}_0.$$

$$R+T=1+\frac{2t{D}^2(t{D}^2-\Gamma )}{{\Delta }^2+{\Gamma }^2}.$$

$$E_n=e^{i|z_n-z_m|}T(|z_n-z_m|,a)d_m\equiv G(n,m)E_m.$$

$$G(n,m)=-\frac{3T(|z_n-z_m|,a)}{2[\delta +\frac{3}{2}S(a)+i]}{e}^{i|z_n-z_m|}.$$

$$E_n=E_0{e}^{i{z}_n}+\sum _{m\ne n}G(n,m)E_m.$$

$$E(z)=E_0{e}^{iz}+\sum _n{G}_0(z,m)E_m,$$

$$G_0(z,m)=-\frac{3\bar{T}(a)}{2[\delta +\frac{3}{2}S(a)+i]}{e}^{i|z-z_m|}.$$

$$G(n,m)=-i\frac{\Gamma \zeta (n,m)}{\Delta +i\Gamma }{e}^{i\Phi (n,m)}{e}^{i|z_n-z_m|}.$$

$$\zeta (n,m)e^{i\Phi (n,m)}=\frac{T(|z_n-z_m|,a)}{\bar{T}(a)}.$$

Full text

I Introduction

Advances in experimental Keaveney et al. (2012); Pellegrino et al. (2014); Kwong et al. (2015); Bromley et al. (2016); Guerin et al. (2016); Jennewein et al. (2016); Roof et al. (2016); Corman et al. (2017) and numerical Javanainen et al. (1999); Chomaz et al. (2012); T. Bienaimé and Kaiser (2013); Javanainen et al. (2014); Bromley et al. (2016); Sutherland and Robicheaux (2016); Zhu et al. (2016); Javanainen et al. (2017) techniques have revitalized classical electrodynamics of material samples as a topic of frontline research. The idea is that in cold, dense atomic samples dipole-dipole interactions between the atoms mediated by light may lead to cooperative behavior of light-matter systems, and even result in a strongly correlated sample. It transpires that the standard classical electrodynamics of polarizable media Jackson (1999); Born and Wolf (1999) is a mean-field approximation, and as such could be inaccurate in cold dense samples Javanainen et al. (2014); Javanainen and Ruostekoski (2016).

Aside from 3D gaseous samples, cooperation of atoms in a 2D lattice has attracted recent interest Jenkins and Ruostekoski (2012); Bettles et al. (2016); Yoo and Paik (2016); Facchinetti et al. (2016); Facchinetti and Ruostekoski (2018); Shahmoon et al. (2017); Perczel et al. (2017a, b); Asenjo-Garcia et al. (2017); Bettles et al. (2017); Yoo (2018). In fact, studies of metamaterials with 2D arrays of dipoles that mix electric and magnetic dipole moments have made similar theoretical points earlier Jenkins and Ruostekoski (2013), and seem to be ahed of atomic lattices even in the experiments Jenkins et al. (2017), but we focus on the simpler case of ordinary atoms. For the present purposes the observation that a lattice with atoms may reflect back all of the incident light Jenkins and Ruostekoski (2012); Bettles et al. (2016); Shahmoon et al. (2017) is of particular interest. Going down in the number of dimensions, effectively 1D waveguides holding atoms are also a subject of long-standing theoretical and experimental activity Deutsch et al. (1995); Le Kien et al. (2004); Vetsch et al. (2010); Tiecke et al. (2014); Goban et al. (2014); Douglas et al. (2015); Goban et al. (2015); Sørensen et al. (2016); Corzo et al. (2016); Ruostekoski and Javanainen (2016, 2017). Light-induced correlations in 1D Ruostekoski and Javanainen (2016, 2017) has been one particular focus. We will use these results to investigate light propagation in 3D structures formed by stacking 2D lattices.

Given the recent activity in the 2D systems and the overlap of our paper with the existing literature, we open with a broad discussion of how the tenor of the present paper is different. First, as far as the 2D lattices are concerned, we investigate their response and light propagation directly, and mostly exactly. We do not engage in numerical simulations of necessarily finite 2D lattices, nor do we compute cooperative light-atom modes in either finite or infinite 2D lattices. Second, our approach is unapologetically classical electrodynamics. We even describe Zeeman splitting of the atomic levels by using an anisotropic polarizability for the atom. Finally, as a key technical item, in infinite lattices one encounters infinite sums that do not converge absolutely, and present numerical problems that have been addressed in various ways in the past Shahmoon et al. (2017); Perczel et al. (2017b); Bettles et al. (2017). Here we apply the usual exponential convergence factor, and carry out the sums in real space using numerical-analysis techniques to accelerate the convergence.

While the correspondence to the layout of the present paper is not one-to-one, our narrative runs as follows: We start with one 2D square lattice of atoms. In this case the response of the lattice to light can be studied almost analytically, except for certain infinite sums Yoo and Paik (2016). The values of the sums notwithstanding, thanks to the Lorentzian form of the response of an individual atom, a lattice responds to the incident light like an atom with modified resonance frequency and linewidth; or in the case of non-normal incidence of light, like an atom with up to three resonance frequencies.

Since stacking of the 2D lattices is a recurring topic here, we study the transfer of radiation from one lattice to another. While also presenting explicit counterexamples, we find that a simple approximation whereby the lattice radiates as if the atomic dipoles were smeared continuously across a plane is often a good description of the radiation from a 2D lattice.

As a side effect of our studies of light propagation between 2D lattices, we revisit the case of complete reflection from a 2D lattice. Cooperation between the atoms is naturally inseparable from the rest of the physics of the system, but the cooperation per se is not responsible for total reflection: The combination of a single Lorentzian resonance and energy conservation suffices for total reflection. Bragg scattering of light from the 2D lattice may lead to seeming violation of energy conservation, so we also discuss Bragg scattering.

We proceed to combine considerations from different dimensionalities. We point out that a stack of 2D lattices can emulate a perfect 1D waveguide with no losses of light. In another example, by viewing a 3D lattice as as stack of 2D lattices, we show how cooperative shifts in 2D and shifts characteristic of 1D waveguides may combine to eliminate completely the density dependent resonance shifts analogous to collisional shifts. Such shifts occur Chang et al. (2004) even if each atom is bound to a lattice site.

In our final example, basically an alternative take of Refs. Facchinetti et al. (2016); Facchinetti and Ruostekoski (2018), we add to our approach a magnetic field and the attendant anisotropic response of the atoms to the driving field. We demonstrate with numerical and analytical arguments that the magnetic field may induce both perfect reflection and perfect transmission of light through a 2D lattice.

A retrospective discussion in Sec. VII closes the present paper.

II Setup

II.1 Electrodynamics primer

We take the light to be classical. This can be rigorously justified for an atom with a $J=0\rightarrow J^{\prime}=1$ transition in the limit of asymptotically low light intensity, and when photon recoil effects are negligible Ruostekoski and Javanainen (1997). The latter is the case if the lattice binds the atoms to regions small compared to the wavelength of light. Of course, the confinement is never perfect, and even if nothing else, the atoms are always subject to quantum mechanical zero-point fluctuations. These could be taken into account, say, by direct numerical simulations Ruostekoski and Javanainen (2017), but here we assume fixed lattice positions of the atoms.

We use custom units Javanainen et al. (2017) chosen in such a way that the numerical values for the wave number of light and for certain natural constants are

[TABLE]

As usual, time dependent quantities are written in terms of the slowly-varying positive-frequency quantities, factoring out the time dependence at the frequency $\omega$ of the driving light. The dipole propagator ${\sf G}$ that gives the electric field at ${\bf r}$ from a dipole with the amplitude ${\bf d}$ at ${\bf r}^{\prime}$ as ${\bf E}({\bf r})={\sf G}({\bf r}-{\bf r}^{\prime})\,{\bf d}$ is $3\times 3$ matrix, a tensor. In our units it has the cartesian components

[TABLE]

Here $\hat{\bf e}_{i}$ are the cartesian unit vectors, and $r$ and $\hat{\bf n}$ are the distance from the source point to the field point and the unit vector directed from the source point to the field point. This tensor is a function of the difference between the coordinates, and satisfies ${\sf G}_{ij}({\bf r}-{\bf r}^{\prime})={\sf G}_{ji}({\bf r}-{\bf r}^{\prime})={\sf G}_{ij}({\bf r}^{\prime}-{\bf r})$.

The final salient part of the electrodynamics has to do with the polarization of the atoms. We scale the detuning, the difference $\omega-\omega_{0}$ between the light frequency and the atomic transition frequency, to the HFHM linewidth $\gamma$ of the optical transition, defining

[TABLE]

For a $J=0\rightarrow J^{\prime}=1$ transition in an atom at low light intensity (and also for a classical isotropic charged harmonic oscillator) the induced dipole moment and the electric field are parallel. With the present units and conventions their quotient, the polarizability, reads

[TABLE]

II.2 Self-sum

We study a square lattice. For definite notation, we say that the squares are aligned with the $x$ and $y$ axes, and $z$ is the direction perpendicular to the lattice. The lattice spacing is denoted by $a$; $a=2\pi$ would be one wavelength. The lattice sites are specified as ${\bf R}_{\bf n}=a{\bf n}=a(n_{x}\hat{\bf e}_{x}+n_{y}\hat{\bf e}_{y})$ with integers $n_{x}$ and $n_{y}$. We denote $n=|{\bf n}|$; ${\bf n}=0$ means the lattice site at the origin. Assume identical dipoles ${\bf d}$ at all lattice sites. Ignoring the self-field of the atom at the origin, the field at the origin is

[TABLE]

By lattice translation invariance, this is actually the field at all sites of the lattice. If the dipole is circularly polarized in the plane of the atoms, so obviously is the field, and the now-scalars dipole and electric field are related by

[TABLE]

Since the result holds for both left- and right-circularly polarized quantities, it actually holds for any polarization of the dipole and the electric field in the plane of the lattice.

The issue here is that the sum is not absolutely convergent. We cure this with the usual convergence factor $e^{-\eta\rho}$, and let $\eta\rightarrow 0+$ at the end of the calculation. The other problem from the point of view of numerics is the infinite sum. Our task therefore is to carry out numerically the operations

[TABLE]

We call $S(a)$ the self-sum. We have developed code that computes $S(a)$ to a preset numerical precision. The technical idea is to accelerate the convergence to the limit $\eta\rightarrow 0$ using Richardson extrapolation Press et al. (2007). The details are discussed in Appendix A.1.

It turns out that the expression (7) does not necessarily converge in the limit $\eta~{}\rightarrow~{}0$. When lattice spacing is increased, this happens for the first time at $a=2\pi$, when the lattice spacing equals the wavelength of the driving light. The reason is Bragg scattering. For $a=2\pi$ light scattered from all atoms interferes constructively in the $x$ and $y$ directions, which in the limit of an infinitely large lattice produces an infinite field strength. Such in-plane Bragg scattering happens every time there are parallel lines of lattice points with the shortest distance between the lines equal to a wavelength. Elementary solid state physics adapted from three to two dimensions shows that this happens at lattice spacings $a=2\pi\sqrt{m_{1}^{2}+m_{2}^{2}}$, where $m_{1}$ and $m_{2}$ are integers, and not both equal to zero. A comprehensive discussion of Bragg scattering is deferred to Appendix B.

Figure 1 shows the self-sum, real (dashed blue line) and imaginary (solid red line) parts as a function of the lattice spacing scaled to $2\pi$ 111There is a functionally equivalent figure in the Supplemental Material of Ref. Shahmoon et al. (2017). The lowest Bragg reflections are expected at $a/2\pi=1$, 1.41, 2, and 2.24, and those are indeed the exceptional values of the lattice spacing in the figure. The real and imaginary parts are peculiar in that they abruptly jump to infinity when a Bragg-reflection lattice spacing is crossed from above (real part) or below (imaginary part), and diverge smoothly on the other side. Although it is not very obvious from the figure, at small $a$ the real and imaginary parts, respectively, diverge as $1/a^{3}$ and $1/a^{2}$.

II.3 Transfer sum

Let us now consider the radiation out of the plane of the lattice. Specifically, imagine an identical second lattice displaced by the distance $z>0$ from the original lattice. By virtue of lattice translation invariance, the total field from the dipoles of the original lattice is the same at all sites of the second lattice. The field equals $E(0)=T(a,z)d$, where the transfer sum $T$ is

[TABLE]

The transfer sum explicitly includes the center site ${\bf n}=0$, whereas in the self-sum the site ${\bf n}=0$ is included only indirectly in the radiative damping rate and the Lamb shift of the atom. Second, the distances between the source and target points $r$ depend on the translation $z$ between the two lattices. Third, the original lattice is expected to radiate a field with a propagation phase $e^{iz}$ toward the increasing coordinate $z$ and $e^{-iz}$ toward the decreasing coordinate. We have explicitly canceled the free-propagation phase from the transfer sum. Because of the $|z|$ in the exponential, the expression (8) is actually valid for both signs of the translation $z$. Finally, the sum (8) is no more absolutely convergent than (6). This is handled by adding a long-distance convergence factor $e^{-\eta\rho}$ just like in Eq. (7). The numerical issues are essentially the same as those encountered with Eq. (7).

There is also an analytical argument to be made. Namely, for any $z\neq 0$ and in the limit $a\rightarrow 0$, the terms in the transfer sum vary less and less with the changes of the lattice site index $\bf n$ between neighboring sites, and the sum may be approximated by the integral over $\bf n$. In the integral it is possible to do the limit $\eta\rightarrow 0$ analytically as well. All of this makes an interesting exercise in Mathematica, and the result is

[TABLE]

independently of the value of $z$. A priori, one expects that this form of the transfer sum is accurate whenever $|z|\gg a$.

A moment’s reflection shows that $\bar{T}(a)d$ is the electric field that would ensue if there were a continuously distributed dipole moment in the plane, with the surface density $d/a^{2}$. Deviations of $T(z,a)$ from $\bar{T}(a)$ can be attributed to the discreteness of the dipoles in the lattice. The lumped nature of the dipoles shows qualitatively in two aspects. First, the transfer sum also has Bragg reflection singularities for the same values of $a$ as the self-sum. Second, in the limit $z\rightarrow 0$ with a fixed $a$, the transfer sum obviously converges to the same value as the self-sum, except for the added contribution from the dipole at ${\bf n}=0$. The field of this dipole at the distance $|z|$ diverges with $z\rightarrow 0$, so that $T(z,a)$ diverges with $z\rightarrow 0$ for any fixed $a>0$.

We illustrate the behavior of the transfer sum in Fig. 2. The range of either argument $z$ or $a$ does not reach down to zero or up to $2\pi$, so there are no outright divergences. However, the onset of the divergence with $z\rightarrow 0$ for $z<a$ is clearly visible, and hints of the exceptional behavior for $a=2\pi$ are also discernible in the 3D plot of the absolute value of $T(z,a)/\bar{T}(a)$ (top panel). For the most part, though, the value of $T(z,a)/\bar{T}(a)$ is close to one. This is demonstrated in the contour plot of the real part (bottom panel) that only shows the region of $z$ and $a$ with $0.9<\Re[T(z,a)/\bar{T}(a)]<1.1$. For $z\gtrsim a$, $\bar{T}(a)=2\pi i/a^{2}$ is a passable approximation for the transfer sum $T(z,a)$. In fact, except for problems associated with Bragg reflections, at large distances $z$ between the lattices, $\bar{T}(a)$ becomes an increasingly accurate approximation to the transfer sum $T(z,a)$. Moreover, at large distances the radiation from the lattice obviously makes a plane wave. $\bar{T}(a)d$ is a good approximation for the electric field not only at the sites of a second lattice, but in the far field also in an entire plane parallel to the lattice.

II.4 Total reflection from a lattice

We now combine what we have so far for an investigation of light transmission and reflection for one 2D lattice. Take an incoming field with the amplitude $E_{0}$ coming perpendicularly to the lattice, then the field $E$ and the dipole moment $d$ at each lattice site satisfy

[TABLE]

see Eq. (3) for the polarizability $\alpha$. This allows us to solve the dipole moment as a function of the incoming field as

[TABLE]

where

[TABLE]

are the effective detuning and transition linewidth. The dipole moment is as for a single atom, except that the cooperative response of the lattice has shifted the resonance by -\hbox{\frac{3}{2}}\Re[S(a)], converting the detuning to the effective detuning $\Delta$, and the linewidth is changed to $\Gamma$.

Next consider the reflected and transmitted fields far away from the lattice, plane waves such that the long-distance transmission amplitude $\bar{T}(a)$ applies. We write

[TABLE]

The total reflected and transmitted amplitudes are then

[TABLE]

Since $\Gamma$ is positive, the transmitted field equals zero if and only if both the effective detuning $\Delta$ equals zero, and the imaginary part of the numerator in $E_{T}$ is also zero,

[TABLE]

This can be cast in the form

[TABLE]

As can be verified by comparing against numerical results, this actually is the dependence of the imaginary part of the self-sum on the lattice spacing in the interval $a\in(0,2\pi)$. Others have shown the same using analytical arguments Shahmoon et al. (2017).

Thus, while $a\in(0,2\pi)$ holds true, there is always a detuning $\delta$ such that $\Delta=0$, and no light gets transmitted through the lattice. As may be seen from Fig. 1, there are two values of the lattice constant in the interval $a\in(0,2\pi)$ with $\Re[S(a)]=0$, and hence for on-resonance excitation with $\delta=0$ also $\Delta=0$ holds true. These are the lattice spacings, approximately $0.2\times 2\pi$ and $0.8\times 2\pi$, for which the lattice does not transmit any resonant light Bettles et al. (2016); Shahmoon et al. (2017).

Our development also provides unique insights into the zero transmission, and the full reflection that accompanies it. To begin with, let us define the power transmission and reflection coefficients

[TABLE]

Naively, $R+T$, the sum of reflection and transmission coefficients, equals one if energy is conserved, as it should be since in steady state all electromagnetic energy coming to each atom also gets radiated away. In Fig. 3 we show the sum $R+T$ (solid black line) and the transmission coefficient $T$ (dashed red line) as a function of the lattice spacing. The detuning $\delta$ is always chosen so that the effective detuning equals $\Delta=0$. The difference of the two curves, of course, is the reflection coefficient $R$. The sum $R+T$ is not always equal to one, seemingly indicating that energy is not conserved.

The explanation is that $R$ and $T$ only incorporate the energy that propagates as a plane wave in the direction perpendicular to the lattice. However, when nontrivial Bragg scattering is possible for $a\geq 2\pi$, there are always Bragg scattered waves present that transport energy to other directions, and these are not included in the sum $R+T$. Compare Figs. 1 and 3: there are discontinuities in the curves at exactly those $a$ for which a new Bragg order emerges.

Now, assume that the dipole moment associated with the lattice shows a generic single-Lorentzian resonance of the form

[TABLE]

where $D^{2}$ is a positive constant, $\Delta$ is a measure of detuning from resonance, and $\Gamma$ is the linewidth. Furthermore, assume that in the far field the radiation is out of phase with the dipole moment by $\pi/2$, so the radiated field is of the form

[TABLE]

with some positive $t$. On the side of the reflection, this is the total field. On the side of the transmission, one adds the incoming field,

[TABLE]

The sum of reflection and transmission coefficients can be manipulated into the form

[TABLE]

This can be identically equal to one only if $\Gamma=tD^{2}$. But this is also precisely the condition that there exists a detuning such that the transmission and reflection coefficients are equal to zero and one, respectively. The other assumptions except a single Lorentzian resonance are generic to atoms interacting with light, and even have obvious counterparts in the general theory of linear response.

In the case of a single Lorentzian resonance energy conservation per se dictates that there must be a detuning for which total reflection prevails for a 2D lattice on normal incidence. Energy conservation works in this way because for plane wave excitation, and in the absence of Bragg scattering, the interference from the lattice atoms forces the radiation from the lattice to propagate either in the direction of the incoming light, or in the opposite direction. Cooperation between the atoms is not directly needed, although it evidently must be part of a consistent theoretical description. A single Lorentzian resonance is a sufficient condition for total reflection, but it will turn out below, Sec. VI, that it is not a necessary condition.

Related arguments have been made before Facchinetti et al. (2016); Facchinetti and Ruostekoski (2018), but our logic is different. We place the Lorentzian form of the resonance at front and center, which is the cue for our further developments about non-normal incidence in Sec. IV.

III Stack of lattices

We next consider the situation when $N$ 2D square lattices with the lattice constant $a$ are stacked so that the $x$ and $y$ coordinates of the atoms are the same for all lattices. The driving light propagates in the perpendicular $z$ direction. The positions of the lattices are denoted by $z_{n}$. For a regularly spaced stack with the spacing $\Delta z$ we could write $z_{n}=(n-1)\Delta z$, $n=1,2,\ldots,N$, but almost nothing in our formal development depends on the regular spacing.

We know that a 2D lattice can completely reflect back all of the incoming light, just like a single atom in a 1D channel for light Javanainen et al. (1999); Ruostekoski and Javanainen (2016, 2017). One would thus think that the existing understanding of atoms in an effectively 1D waveguide could be transferred to a stack of lattices. This notion is by and large correct. The basic difference is that in a 1D waveguide there is only one independent value of the electric field (for a given polarization mode) at every 1D coordinate, call it $z$. The response of an atom will impose a condition tying the incoming and scattered lights on both sides of the atom, and one can set up a transfer-matrix analysis Ruostekoski and Javanainen (2017). Unfortunately it is not so for a stack of 2D lattices. The electric field falling on the sites of a lattice does not uniquely determine the electric field in the whole plane of the lattice, and therefore neither the field that would excite the next lattice. This will directly and indirectly cause differences between a 1D waveguide and a stack of 2D lattices. Our question is, how much of a difference?

Now, given the electric field $E_{m}$ at the $m^{\rm th}$ lattice, the lattice atoms develop a dipole moment $d_{m}$ as before, Eq. (11). The lattice then produces a scattered field whose amplitude on the atoms of a lattice at $z_{n}$ equals

[TABLE]

Here we have put back the propagation phases of light that were transformed away in the definition of the transfer sum. The propagator for the light is

[TABLE]

Given an incoming plane wave with the amplitude $E_{0}$, the fields at the sites of the lattices are related by

[TABLE]

This is an inhomogeneous set of linear equations for the fields $E_{n}$, which can be solved. Given the fields, we can compute the dipoles $d_{n}$, and from them the total transmitted field. Far behind the stack the long-distance form of the transfer sum $\bar{T}(a)$ is operative, and at large $z$ we have the transmitted field

[TABLE]

with

[TABLE]

Transmission and reflection coefficients of the stack of lattices may be computed from the far field (25).

Suppose the lattice spacing satisfies $a\in(0,2\pi)$, so that (15) holds true. We may write Eq. (23) in the form

[TABLE]

Here $\zeta>0$ and $\Phi\in\mathbb{R}$ are the modulus and the phase of the complex ratio that characterizes the difference between the transfer sum and its long-distance form,

[TABLE]

The counterpart of the propagator (27) for atoms in a 1D waveguide Ruostekoski and Javanainen (2016, 2017) is, in the present units,

[TABLE]

Here $\zeta$ is the fraction of the light energy emitted by an atom that ends up back in the waveguide, as opposed to leaking to free space. The cooperative resonance frequency and linewidth of the lattice, instead of the resonance quantities of a single atom, are the obvious differences between a lattice and an atom in a 1D waveguide. Moreover, in a 1D waveguide there would be no adjustments $\Phi$ to the propagation phases, and, in lieu of the factors $\zeta(n,m)$ that depend on the positions of both lattices, a (homogeneous) 1D waveguide would have a constant $\zeta$. For an ideal waveguide, $\zeta=1$.

Since $T(|z_{n}-z_{m}|,a)\rightarrow\bar{T}(a)$ when either $a\rightarrow 0$ or $|z_{n}-z_{m}|\rightarrow\infty$, both of these limits produces a simulacrum of a perfect 1D waveguide. As far as light transmission through a 1D waveguide goes, the distances between the successive atoms only matters modulo half of a wavelength, $\pi$ in our present units Ruostekoski and Javanainen (2016, 2017). Thus, by placing the lattice planes far enough apart, it is in principle possible to emulate the physics of an arbitrary string of atoms in a perfect 1D waveguide to an arbitrary precision using a stack of 2D lattices.

The final question is about the quality of the emulation. In Fig. 4 we plot the transmission as a function of the scaled detuning $\Delta/\Gamma$ for a stack of $N=4$ lattices. Figures 5 and 6 of Ref. Ruostekoski and Javanainen (2017) shows analogous graphs for a 1D waveguide, though none of them for the exact same parameter values.

The scaling $\Delta/\Gamma$ removes the cooperative resonance shift and broadening of each 2D lattice from consideration. The dashed red line is for the geometric constants $a=\Delta z=\pi$, half of a wavelength, and shows the superradiant resonance with the width $4\,\Gamma$. The solid line is for the constants $a=\Delta z=\pi/2$. It displays a sharp stop band with two narrow transmission features due to Fano resonances in the cooperative response of the lattices in the stack. As far as we can tell numerically, these resonances still have exactly unit transmission, but the positions are slightly shifted from what one would expects on the basis of atoms in an ideal 1D waveguide: $\Delta/\Gamma\simeq\pm 1.45$ versus $\delta=\pm\sqrt{2}\simeq\pm 1.41$. To an unaided eye, the graphs in Fig. 4 would be indistinguishable from the corresponding graphs for an ideal 1D waveguide.

Close to $a=2\pi$, starting at about 15% below a wavelength, and on both sides of $a=2\pi$, we see substantial deviations from the picture of an ideal 1D waveguide, even after scaling away the cooperative effects of each individual 2D lattice by plotting the results as a function of $\Delta/\Gamma$. An example is provided in Fig. 5. This shows the transmission for a stack of four lattices at the fixed spacing $\Delta z=\pi/2$, for different lattice constants. At $a=\pi$ , solid black line, the result is similar to the result shown in Fig. 4 for $a=\Delta z=\pi/2$, but the other two graphs for $a=1.8\,\pi$ (dashed red line) and $a=2.2\,\pi$ (dotted blue line) deviate from the corresponding 1D fare substantially.

We seek further illumination from time evolution. Let us denote the total field strength at lattice $n$ by $E_{n}$, then the linear dipole at site $n$ satisfies the time evolution equation

[TABLE]

The steady-state solution gives $d_{n}=\alpha E_{n}$ with the polarizability (3), as it should. We also state the assumption that the propagation time of light across the sample is much smaller than the response time of the atoms, something like the inverse of the linewidth. Then the field $E_{n}$ satisfies

[TABLE]

as an instantaneous equation. This gives a coupled equation for the dipoles

[TABLE]

The inhomogeneous version of Eq. (32) with $E_{0}=0$ gives the cooperative radiative modes of the stack of lattices. Specifically, with the ansatz $d_{m}(t)=\tilde{d}_{m}e^{-i\omega t}\equiv e^{(iu-v)t}$, we have the eigenvalue equation

[TABLE]

The effective detuning and linewidth for each mode, $u$ and $v$, incorporate the cooperative effects from the radiative coupling between the lattices. As usual, the corresponding eigenvectors made of $\tilde{d}_{n}$ determine the coupling of the eigenmodes to the driving plane wave.

Once more, in the case $a\in(0,2\pi)$ and under the generally quite good an approximation $T(|z_{n}-z_{m}|,z)=T(a)$, the time dependent collective modes and collective behavior of the stack of lattices would be the same as the behavior of the atoms in a lossless 1D waveguide, save for the lattice shift and line broadening in the detuning $\Delta$ and linewidth $\Gamma$.

However, $T(z,a)\simeq\bar{T}(a)$ breaks down around $a\simeq 2\pi$. What happens upon approach, and crossing, the limit $a=2\pi$ is shown in Figs. 6. This is for a stack of $N=4$ lattices, for the constant spacing between the lattices $\Delta z=\pi$. The corresponding ideal 1D waveguide would have four eigenmodes with no resonance shifts, $u=0$, whereas one eigenmode would have a superradiant decay rate $v=4$ and the other three would have the extreme subradiant decay rate $v=0$ Ruostekoski and Javanainen (2017). The broad resonance in Fig. 4 (dashed red line) originates from the coupling of the driving light to the superradiant mode.

Figure 6 shows the decay rates (top) and resonance shifts (bottom) of the eigenmodes as a function of the lattice constant $a$. The detuning is chosen so that $\Delta=0$. The results are presented in units $\Gamma$, and the effects of the 2D lattice per se have again been renormalized away. The way we have defined the detuning, the cooperative resonance shift actually equals $-u$, as given in the figure. Below but up to quite close to $a=2\pi$, the analog of the 1D waveguide prevails. The behavior for $a>2\pi$ is different, but the damping rates (in the units $\Gamma$, which is singular in itself) appear to interpolate continuous across the Bragg reflection singularity at $a=2\pi$.

In summary, a stack of 2D lattices generally speaking behaves closely like a 1D ideal waveguide with embedded atoms, although these quasiatoms have the resonance frequency and the radiative linewidth of the 2D lattice. The notable exceptions we have demonstrated occur at lattice spacings such that 2D Bragg scattering is, or is about to become, a consideration.

IV Arbitrary angle of incidence

Up to now, the light has always propagated in the direction perpendicular to the lattice. We now take up the case of an arbitrary direction of the incoming light.

IV.1 Formulation

Henceforth we express the direction of the incoming plane wave using spherical polar coordinates. With our choice of the units the length of the wave vector of the incoming light equals one, so we write

[TABLE]

$\theta=0$ means perpendicular incidence, and we assume that the incoming light propagates from “left” to “right” so that $0\leq\theta<\pi/2$. The complex polarization vector of the incoming beam is denoted by $\hat{\bf e}$. For a valid incoming plane wave it must be that ${\bf k}\cdot\hat{\bf e}=0$. Here the dot $\cdot$ stands for the inner product as for real vectors.

The obvious effect of the tilt is to cause a rolling phase for both the electric field and the dipoles according to the projection of the wave vector on the plane of the atoms, ${\bf k}_{\parallel}=\sin\theta\cos\phi\,\hat{\bf e}_{x}+\sin\theta\sin\phi\,\hat{\bf e}_{y}$, but a single vectorial amplitude still characterizes both the dipole moment and the electric field at the positions of the atoms. We write the dipolar field, excluding the self-field, at the center point of the lattice in terms of a self-sum tensor,

[TABLE]

The matrix $\sf S$ is symmetric. Moreover, symmetry considerations suggest, and explicit calculations confirm, that the $z$ components of the light field and of the dipoles do not couple to the $x$ and $y$ components, giving ${\sf S}_{xz}={\sf S}_{yz}=0$.

Likewise, the field transferred to the center site of a lattice at the distance $z$ away is characterized by what we call (somewhat confusingly, given the standard terminology) the transfer matrix $\sf T$,

[TABLE]

The purpose of the exponential prefactor is to cancel the free-propagation phase over the distance between the planes; $k_{\perp}=\sqrt{1-{\bf k}_{\parallel}^{2}}=\cos\theta$. The matrix $\sf T$ is symmetric as well. The behavior of the tensor ${\sf T}$ under the reflection $z\rightarrow-z$ warrants caution: the $z$ component of the electric field radiated by the components of the dipoles that lie in the plane of the atoms have opposite signs on the opposite sides of the plane, but the field from the $z$ components of the dipoles is the same at both $z$ and $-z$. This symmetry implies that ${\sf T}_{xz}(z)=-{\sf T}_{xz}(-z)$ and ${\sf T}_{yz}(z)=-{\sf T}_{yz}(-z)$, but these and their matrix-symmetry counterparts are the only elements of ${\sf T}$ that change upon the reflection $z\rightarrow-z$.

It turns out that the transfer matrix ${\sf T}$ also has a large-distance/dense-lattice limit analogously to the transfer sum $T$. By comparing with numerical computations we have found that for positive $z$ the limit is

[TABLE]

with

[TABLE]

see Appendix A.2.

Particularly interesting from our viewpoint is that the large-distance limit of the transfer matrix may also be written as

[TABLE]

where ${\sf P}_{\perp}({\bf k})$ stands for the projector to the subspace orthogonal to the vector ${\bf k}$. This is how the formalism enforces the condition that the reradiated field is transverse. When the lattice is viewed at the angle $\theta$ away from perpendicular, the apparent area density of the dipoles is ${\bf d}/(a^{2}\cos\theta)$. This neatly explains the $\cos\theta$ in the denominator of (46).

Bragg reflections are a possibility for oblique incidence, even more so than for normal incidence. This is discussed in detail in Appendix B. However, unless Bragg reflections are specifically mentioned, we continue to assume that they are absent.

IV.2 Single plane of atoms

Special cases proliferate when we allow non-perpendicular incidence. We avoid detailed descriptions of quantitative results, and concentrate on prominent qualitative features of reflection and transmission from a single plane of atoms. Our observations are extracted from, or supported by, numerical calculations. Although we do not elaborate on this, the matrix $\sf S$ inherits symmetries from the underlying square lattice that help to discover the results. It is not an accident that we occasionally call the matrices $\sf S$ and $\sf T$ “tensors.”

We have the incoming field ${\bf E}_{0}$ and the total field ${\bf E}$ at the center site related by

[TABLE]

where $\sf A$ is the polarizability tensor. The response of an atom with the $J=0\rightarrow J^{\prime}=1$ transition in the absence of an external magnetic field is isotropic, and $\sf A$ simply equals polarizability times the unit matrix. However, we retain the option of anisotropic atomic response in our formalism for applications below, in Sec. VI. It is a simple matter to solve the total field and the atomic response,

[TABLE]

We write the reflected and transmitted fields on the left ($z<0$) and on the right ($z>0$) in the form

[TABLE]

In the rest of this section we always use the large-$z$ limit $\bar{\sf T}$ for the transfer matrix. The reflected field propagates in the direction

[TABLE]

Correspondingly, the $z\rightarrow-\infty$ limit of the transfer matrix $\bar{\sf T}_{R}$ projects to the subspace orthogonal to ${\bf k}_{R}$.

So far we take the polarizability tensor $\sf A$ to be diagonal, corresponding to the usual isotropic polarizability $\alpha$ as in (3). The induced dipole is then

[TABLE]

The matrix $\sf S$ has three eigenvalues, and (without proof) we expect it to have three linearly independent eigenvectors. The component of the electric field along the eigenvector corresponding to the eigenvalue $s_{i}$ drives a resonance with the effective detuning and linewidth

[TABLE]

Now, take any incoming beam of light characterized by the polarization vector $\hat{\bf e}$. The polarization vector may be decomposed into a linear superposition of the eigenvectors of $\sf S$, and there are in general three separate resonances with separate resonance frequencies and linewidths that show up in the response; one for each eigenvector. The matrix $\sf S$ has no elements that couple the $z$ polarization to any polarization in the plane, so at this stage the $z$ components of the dipoles and of the electric field, if any, get processed separately. In the subsequent reradiation stage governed by the transfer matrices $\bar{\sf T}$ and $\bar{\sf T}_{R}$ the eigenvectors are projected to the two-dimensional subspaces of transverse polarizations allowed for the exiting plane waves, but in the most general case the outgoing light still carries along behavior from all three resonances.

For normal incidence, $\theta=0$, there is no difference between the $x$ and $y$ directions, and the incoming field must be polarized in the plane of the atoms. Therefore ${\sf S}_{xx}={\sf S}_{yy}$ and ${\sf S}_{xy}~{}=~{}0$ must hold true, and the dipolar resonance in the $z$ direction cannot be excited. Since the $x-y$ part of the matrix ${\sf S}$ is diagonal, all polarizations in the $x-y$ plane are eigenvectors with the same eigenvalue. We are back to the single-resonance case we have discussed already.

For non-normal incidence there in general is a $z$ component in the polarization of the incoming light, which is subject to its own separate resonance. For a general propagation direction there are also two elliptically polarized eigenmodes of $\sf S$ in the $xy$ plane, which are linearly independent though not necessarily orthogonal in the sense of the proper inner product for complex vectors. We expect two resonances in the $x-y$ plane as a result. The situation simplifies if the projection of the propagation vector ${\bf k}_{\parallel}$ to the lattice plane is along a direction of high symmetry. If ${\bf k}_{\parallel}$ lies in a direction of a lattice axis ($\phi=0,\,\pi/2,\,\ldots$), there are two resonances linearly polarized along the axes. If ${\bf k}_{\parallel}$ bisects an angle between the lattice axes ($\phi=\pi/4,\,3\pi/4,\,\ldots$), there are likewise two eigenvectors of $\sf S$ in the plane polarized parallel and perpendicular to $\bar{\bf k}_{\parallel}$.

Taking into account the transversality of light, we have two guaranteed ways to produce just a single resonance in a square lattice of atoms with oblique incidence: Choose the propagation vector of light in such a way that its component in the lattice plane either lies in the direction of a lattice axis or bisects the axes, and choose the linear polarization of the light in such a way that it lies in the lattice plane. This simplification is based on the symmetries. There may be other symmetries or just plain accidental special cases of a single resonance, but we have not run into any.

At this point the question of whether it is possible to attain zero transmission or perfect reflection for non-normal incidence practically asks itself. First, even in the normal-incidence case perfect reflection can only happen on exact lattice-shifted resonance. We therefore hypothesize that perfect reflection from a single lattice is not possible if more than one resonance gets excited, as the frequency of the light could not be the right one for all resonances at the same time. Our experience supports this hypothesis — the difference from perfect reflection may be small, but we have always been able to detect it. Moreover, we have studied numerically transmission of the eigenmodes of $\sf S$ regardless of whether they are proportional to legal polarizations. We surmise that complete reflection is only possible if the eigenmode, in fact, is a valid transverse polarization for a plane wave, and gets projected untruncated by the transfer matrices $\sf T$ and ${\sf T}_{R}$. These considerations together suggest that the special cases of the preceding paragraph are precisely the cases when one can reach complete reflection. A magnetic field that causes additional mixing of in-plane and out-of-plane polarizations (Sec. VI) changes things, but otherwise our numerical experiments support this hypothesis.

These concepts are illustrated in Fig. 7 showing the transmission and reflection coefficients along with the sum ($R+T=1$) of the two for non-perpendicular incidence. The lattice spacing here is $a=\pi$, half of the wavelength, the light comes in at an oblique angle of $0.4\,\pi$ ($\theta=0$ means perpendicular incidence and $\theta=0.5\,\pi$ would be light coming in along the lattice plane), and the projection of the wave vector onto the lattice plane makes the angle $\phi=0.125\,\pi$ with the $x$ axis, i.e., bisects the angle between the $x$ axis and the direction that, in turn, bisects the angle between the $x$ and $y$ axes. The polarization of the light is in the plane of the atoms. The two resonances with the corresponding positions and widths $-0.325+0.389\,i$ and $0.399+3.00\,i$ are obvious in the figure, as is the fact that the reflection coefficient never reaches one.

Our computations show that, as long as there are no Bragg reflections, the sum of transmission and reflection coefficients again equals one. It should be the case by virtue of energy conservation, but we find it remarkable how much of the mathematics of the tensors $\sf S$ and $\sf T$ must align just right to produce this result. If one grants energy conservation, the perfect reflection in the case of the polarization of the light that is a joint eigenvector of $\sf S$ and $\sf T$ is presumably a simple extension of the argument in the 1D case for a single Lorentzian resonance. One could then even find the corresponding damping rates $\Gamma_{i}$ as a function of the lattice spacing $a$ analytically, analogously to the 1D case.

V Canceling light shifts in a lattice

It is a commonly held notion in high-precision spectroscopy that binding the atoms to a lattice will eliminate collisions, and thereby the associated collision shifts. This idea need not be an unqualified success Chang et al. (2004). We have already demonstrated cooperative shifts of the resonance originating from the dipole-dipole interactions of the atoms in a 2D lattice in Fig. 1. The same dipole-dipole interactions are largely responsible for collision shifts Lewis (1980). Moreover, the analysis of regularly spaced atoms in a 1D waveguide shows possibly very large shifts of the resonances when the spacing between the atoms is close to half of the wavelength Ruostekoski and Javanainen (2017), which would be the lattice spacing if one could use light close to resonance to make an optical lattice to hold the atoms. This is because the atoms effectively delimit cavities, which pull the resonance. At present we view a 3D lattice as a 1D stack of 2D lattices, and argue that these 1D and 2D shifts may cancel if the geometry of the experiment is chosen properly.

We first revisit the relevant material in 1D Ruostekoski and Javanainen (2017), so far in fully dimensional quantities. The lattice has the spacing $d$. In addition there is the retention ratio $\zeta$, as also described in Sec. III. In the limit of a large number $N$ of atoms and a small retention ratio $\zeta$ such that $N\zeta\lesssim 1$ (whereupon the optical thickness of the sample is at most comparable to unity), the resonance of a 1D waveguide experiences a cavity shift

[TABLE]

It may be very large for an inopportune lattice spacing.

Next return to 2D and 3D lattices. If every site is occupied, on resonance even a very small sample is opaque, and in addition to the cooperative line shifts there could be large line broadenings as a result of the optical thickness. Having only a small fraction of the lattice sites occupied at random might be an advantage.

Unfortunately, exact modeling of the random atomic positions most likely requires direct numerical simulations. We attempt to get away with an ad-hoc model. Suppose the filling factor of the lattice is $\zeta$. If we, on the other hand, multiply the polarizability by an adjustable constant $\zeta$, $0<\zeta\leq 1$, the effect is to reduce the induced dipole moments by the factor $\zeta$, as if the area density and volume density of the atoms were both reduced by the factor $\zeta$. We equate these two roles of the constant $\zeta$, which is a mean-field type approximation.

As the final ingredient of the model, we use the large-distance limit of the transfer sum $T$ or the transfer matrix $\sf T$. This effectively means that, in transmission of light, a lattice behaves as if the dipole moment were distributed continuously over the lattice plane. This is also a mean-field approximation.

Now consider a stack of 2D lattices with the lattice constant $a$ and lattice spacing $d$ a normal incidence. When the polarizability is reduced by the factor $\zeta$, this is just as in a 1D waveguide holding atoms with the fraction $\zeta$ of the radiated energy retained; compare with Sec. III. Reverting to the units of the present paper, such atoms in a 2D lattice produce an atom-like response with the linewidth \Gamma=1+\hbox{\frac{3}{2}}\zeta\,\Im[S(a)] and lattice-induced resonance shift s_{L}=-\hbox{\frac{3}{2}}\zeta\,\Re[S(a)]. Let us assume that the filling factor is small enough that $\zeta\,\Im[S(a)]$ is small compared to unity, so the linewidth stays approximately unchanged. By comparing with Eq. (54), we have the total shift at normal incidence

[TABLE]

Next move on to non-normal incidence, applying the development in Sec. IV. If the purpose of putting the atoms in a lattice is to improve precision of spectroscopy, one would not want to introduce two or even three separate resonances with their own shifts. We assume that the non-normal incidence is in a configuration for which only one resonance with the complex eigenvalue $s$ gets excited. Besides, incidence at the angle $\theta$ reduces the effective wave number for propagation from lattice plane to lattice plane by the factor $\cos\theta$. We thus have an expression for the resonance shift

[TABLE]

For a numerical illustration we consider a simple cubic lattice with $d=a$. The lattice spacing equal to half of the magic wavelength $813\,{\rm nm}$ in strontium with the transition wavelength $698\,{\rm nm}$, which gives $a=3.66$. We direct the light so that the projection of the propagation direction bisects the angle between lattice axes, $\phi=\frac{1}{4}\,\pi$. The total resonance shift $\delta_{T}$ as a function of the angle of incidence $\theta$ is plotted in Fig. 6. There are two angles where the shift gets very large owing to the cotangent behavior of a 1D lattice, but also two angles at which the shift cancels. In the other case of a single resonance with ${\bf k}_{\parallel}$ along a lattice axis, an angle of incidence with no density shift exists only after (2D) Bragg scattering would have set in already.

Our argument could be questioned on various grounds. In particular, we have earlier produced an example Javanainen and Ruostekoski (2016) in which mean-field theory for atom-light interaction fails by a factor of ten. We break up this issue in two parts: Can a partly filled lattice plane be thought of as kind of a superatom with mean-field parameters? Is mean-field theory sufficient to describe propagation of light between the superatoms?

Examples we have worked out numerically in Appendix C suggests that the mean-field approximation for a 2D lattice is semiquantitatively valid. Two factors may help here. First, in a typical application the lattice spacing is about half of a wavelength, and the density expressed as a dimensionless number is at most $\sim(1/2)^{3}\sim 0.1$. One can detect deviations from mean-field theory even at such “low” densities  Javanainen et al. (2014, 2017); Javanainen and Ruostekoski (2016), but no qualitative failure. Second, the lattice sets a minimum distance between the atoms and thereby limits the contributions from the $1/r^{3}$ and $1/r^{2}$ parts of the dipole-dipole interaction.

In a 3D cubic lattice the self-sum and transfer sum are very similar in structure. Encouraged by the success of the mean-field approximation within a plane, we therefore surmise that the mean-field theory for light propagation between the planes also remains valid even if there are empty sites. Overall, the errors coming with the mean-field approach should not be large enough to invalidate our conclusions.

While putting atoms in a lattice will not automatically get rid of dipole-dipole interactions and the ensuing analog of collision shifts, a judicious choice of the experimental parameters, including the orientation of the lattice, could eliminate density dependent resonance shifts when the atoms are held in an optical lattice.

VI Magnetic field

In the present section we consider the effects of the magnetic field on the response of the lattice. The studies in this direction were initiated in Refs. Facchinetti et al. (2016); Facchinetti and Ruostekoski (2018), which propose to use 2D atom lattices for magnetometry. Many of the results we bring up were already introduced in these papers, although the reasons for the results in our treatment are occasionally different in subtle ways.

The technical item of the present section is to give up on isotropic polarizability. The polarizability tensor $\sf A$ may be obtained from quantum mechanics for an arbitrary atomic level scheme and direction of the magnetic field, a computation that is implemented numerically as part of the “Software Atom” Javanainen (2017). For a $J=0\rightarrow J^{\prime}=1$ transition it may even be obtained from a classical analysis of the motion of a charge bounded to a harmonic oscillator potential, and acted on by both a sinusoidally oscillating electric field and a static magnetic field. We outline this argument in Appendix D. The numerical parameter that encapsulates the strength of the magnetic field, $\omega_{B}$, equals the Zeeman splitting in the excited level in units of the linewidth. It is defined here based on the frequency difference between the states $m^{\prime}=1$ and $m^{\prime}=0$, and is negative if the Lande factor of the level $J^{\prime}=1$ is positive.

Now that the polarizability is a nontrivial $3\times 3$ matrix, we need to use the full formalism for the reflected and transmitted fields as in Eqs. (48)-(50), and in the most general case we may have to consider an arbitrary direction and polarization of the incident field as well. We have found in numerical examples with full sets of basically random parameters that, in the absence of Bragg scattering, energy is again conserved, $R+T=1$.

We discuss as an example the special case when the light is incident perpendicularly to the atomic plane and is linearly polarized in the $x$ direction, while the magnetic field is in the $y$ direction. The motion of the charges induced by the magnetic field is in the $z$ direction, perpendicular to the atomic plane and in the direction of propagation of light. The dipole moment will have a component in the $z$ direction as well.

Facchinetti et al. Facchinetti et al. (2016); Facchinetti and Ruostekoski (2018) have studied this situation using the eigenvalues and eigenvectors of the joint atom-field system along similar lines as in our discussion around Eq. (32). They note the existence of a subradiant mode in which the atoms store the electromagnetic field for extended periods of time. The proximate reason is that the atoms get polarized in the direction of propagation of the light, whereupon they do not radiate. The authors also introduce a two-mode model that, with proper extrapolation to the limit of infinite lattice, is mathematically equivalent to our method.

As an example, consider Fig. 9. We set the lattice spacing $a=1.1\,\pi$ and the Zeeman splitting $\omega_{B}=-0.15$ (solid red line) or $\omega_{B}=-1.1$ (blue dashed line), and compute the reflection coefficient numerically. This figure is our counterpart of the right panel of Fig. 4 in Facchinetti et al. (2016). In particular, the dip in both curves indicates excitation of the subradiant cooperative mode. Unlike in the figures in Refs. Facchinetti et al. (2016); Facchinetti and Ruostekoski (2018) showing results for finite-size lattices, in our results for an infinite lattice the dips in both curves reach right down to zero. This indicates no reflection, and perfect transmission. Also, as illustrated with the help of the dotted black line at $R=1$, the reflection coefficient always (for $\omega_{B}\neq 0$) touches unity at two detunings; one on each side of the detuning giving the complete transmission.

We may easily gain semi-analytical insights into the behaviors of reflection and transmission. To begin with, neither the polarizability tensor $\sf A$ nor the self-sum tensor $\sf S$ couple the $y$ components of vectorial quantities to the $x$ and $z$ components, so we only need a two-dimensional description. The required tensors are

[TABLE]

where the elements of the self-sum tensor $\sf S$ need to be computed numerically, and depend only on the lattice spacing $a$. We have the complete solution, electric field amplitudes and dipole moments, from Eqs. (48)-(50).

Before proceeding to the results we note that these solutions are the same as the steady-state solutions to the two-mode model of Refs. Facchinetti et al. (2016); Facchinetti and Ruostekoski (2018), except that in these papers the counterparts of the coefficients $S_{xx}$ and $S_{zz}$ are primarily extracted from the from the eigenvalues of the time evolution of the finite lattice system. The authors then proceed to extrapolation to an infinite lattice and present results either the same or similar to ours, although there are differences in both the logic and the use of approximations (of which we have none). References Facchinetti et al. (2016); Facchinetti and Ruostekoski (2018) show no complicated polarizability tensors, and it might come as a surprise that the models are materially equivalent. The explanation is that the steady state is already built into the polarizability tensor. It is not valid while electric field and polarization depend on time, whereas both time dependent relations between electric field and polarization and steady state could be solved from the differential equations in Refs. Facchinetti et al. (2016); Facchinetti and Ruostekoski (2018). On the other hand, since the system is linear, frequency dependence of the steady state also uniquely determines the time dependent response, so we could go backward and find the time dependence if needed.

Putting in the incoming field ${\bf E}_{0}=\hat{\bf e}_{x}$, the radiated (far field) component of the electric field $r=\hat{\bf e}_{x}\cdot{\bf E}_{T}$ is directly the reflection amplitude. We have

[TABLE]

where $\bar{T}(a)$ is given by Eq. (13) as before. The reflection coefficient $R$, as in Fig. 9, equals $R=|r|^{2}$.

The familiar structure of polarizability is emerging here, compare with Eq. (11). In the absence of the magnetic field, $\omega_{B}=0$, the result agrees with the one we already have for perpendicular incidence. The reflection amplitude can equal zero only if the fractional expression in the denominator diverges at some detuning $\delta$, which, in turn, is only possible if $\Im[S_{zz}]=-\frac{2}{3}$. This, in fact, was the case for all values of $a\in(0,2\pi)$ we tried. For a given lattice spacing $a$, there is a detuning equal to -\hbox{\frac{3}{2}}\Re[{\sf S}_{zz}] at which the term proportional to $\omega_{B}^{2}$ diverges, which gives $r=0$, and the lattice appears totally transparent. This holds no matter what the nonzero value of the magnetic field is.

There is a physical explanation for the special value of the imaginary part of ${\sf S}_{zz}$. In the case of perpendicular incidence the self-sum tensor ${\sf S}$ has an eigenmode in the $z$ direction that cannot radiate. By energy conservation, this means that the mode is not damped either. The linewidth analogous to the linewidth for the transverse components of polarization would therefore be $\Gamma_{z}=1+\frac{3}{2}\Im[S_{zz}]=0$, giving $\Im[S_{zz}]=-\frac{2}{3}$. Energy conservation together with the transverse nature of the electric field dictates the presence of subradiance, not cooperation between the atoms per se.

As we have noted already, the imaginary part of ${\sf S}_{xx}$ is also known. Let us write

[TABLE]

obviously with $R_{x}=\Re[S_{xx}]$, $R_{z}=\Re[S_{zz}]$. The reflection amplitude is bound by one in absolute value, and the maximum absolute value 1 is attained at two values of the detuning $\delta$ that depend on the real parts of the coefficients of the matrix $\sf S$ and on the magnetic field:

[TABLE]

We then have both perfect reflection, and zero transmission.

We can calculate $R=|r|^{2}$, $T=|1+r|^{2}$, and $R+T$ from Eqs. (58) and (59), and ask if the result equals one for all detunings $\delta$ and magnetic-field splittings $\omega_{B}$. As another peculiar point of mathematics, this turns out to be the case as soon as the imaginary parts of ${\sf S}_{xx}$ and ${\sf S}_{zz}$, as well as $\bar{T}(a)$, all have their known values. $R_{x}$ and $R_{z}$ have some definite values for each lattice constant, but energy conservation and the main qualitative features of the optical response do not depend on these values.

The $z$ component of polarization may be found similarly to the radiated field. At the detuning when only the $z$ component remains, $\delta=-\frac{3}{2}\,R_{z}$, we have the dipole moment amplitude per atom

[TABLE]

This can be large, and especially so if the resonances associated with the in-plane and perpendicular dipoles coincide, $R_{x}=R_{z}$. We draw these resonance shifts in Fig. 11 as a function of the lattice spacing, $R_{x}$ as a dashed blue line (this is part of the dashed blue line in Fig. 1) and $R_{z}$ as a solid red line. For the lattice spacing of $0.537\times 2\pi$ (compare with the lattice spacing $0.55\times 2\pi$ in Fig. 9) such a confluence of the resonances in fact occurs.

VII Discussion

We have investigated the optical response of atoms bound to a 2D lattice, allowing for non-normal incidence and anisotropic polarizability that occurs when a magnetic field is present. We have also explicitly looked into the issues that come up when one stacks 2D lattices to make 3D lattice structures. Our methods are a mix of analytical and numerical calculations, symbolic computation, and numerical analysis. Especially the algebra associated with Sec. VI is massively tedious, and might be difficult to get right without a tool such as Mathematica. Our amalgam of computer-based methods should be of some interest in its own right.

Of the physics principles that have come up, we emphasize energy conservation. In steady state an atom reradiates all of the energy that it removes from the electromagnetic field. With a plane wave coming in to a 2D lattice of atoms, in the absence of Bragg scattering, interference between the radiators dictates that only transmitted and reflected plane waves can be present in the far field. It is not a particular surprise that under such circumstances the sum of reflection and transmission coefficients equals one in all cases we have tried. What boggles the mind is how the mathematics and numerics with the self-sum tensor, transfer matrix, and polarizability tensor always manage to line up in this way.

We turn the puzzle into a virtue, however. The self-sum and the transfer sum are not absolutely convergent, and by doing them differently could have given very different results. We have tamed the sum with an exponential convergence factor, an old and widely used trick everywhere in theoretical physics. The resulting $R+T=1$ lends some credence to this process.

While studying light transmission in a stack of 2D lattices, we have pointed out the close analogy of the system with atoms inside a 1D waveguide for light. We have argued that it is possible to make the stack to emulate light propagation in the 1D situation, but also discuss the limitations of such a scheme.

The original motivation for the present study was Ref. Bettles et al. (2016) that headlines the observation from numerical simulations that a 2D lattice can be opaque. While one might be tempted to ascribe this behavior to cooperative effects, we have pointed out that the combination of energy conservation and single-Lorentzian form of the resonance line of the lattice suffices for this result. As a matter of principle, cooperativity does not enter the argument at all. Related discussions have appeared in the literature Jenkins and Ruostekoski (2012); Facchinetti et al. (2016); Facchinetti and Ruostekoski (2018), but we have also turned our argument into a tool to discuss the behavior of the lattice of atoms at non-normal incidence.

Indeed, building on the analyses of 1D vs. 2D lattices and non-normal incidence, we have produced an example that might be relevant in experiments. We have shown that, by picking a suitable propagation direction and polarization for the driving light, it may be possible to cancel density dependent resonant shifts for atoms bound to a 3D lattice. This point is not trivial: The atoms, being constrained, do not move and collide, but the dipole-dipole interactions that are a major contribution to collision shifts still remain. In general, density shifts are present even for trapped atoms.

Finally, we have developed quantitative theory of the optical response for a configuration with a magnetic field acting on a lattice of atoms. This scheme was developed in Refs. Facchinetti et al. (2016); Facchinetti and Ruostekoski (2018) starting from finite-size numerical simulations. Besides giving semi-analytical results for an infinite lattice, we show that in this particular configuration any nonzero magnetic field can induce perfect transmission of light, and a total of two tunings of the light exist when the lattice is perfectly reflecting. The transparency can again be traced to back to energy conservation, but the total reflection does not match the single-resonance scheme we have identified for the case with no magnetic field. At present we do not understand the reason for perfect reflection, let alone how it gets embedded into the mathematics.

Our explicit solutions are seemingly limited to narrow special cases. A discussion of the restrictions is warranted.

With some lesser details such as no recoil effects and absence of optical pumping discussed elsewhere Javanainen et al. (2017), the classical electrodynamics we have used is valid basically in the limit of low light intensity Ruostekoski and Javanainen (1997). It is probably impossible in practice to solve this problem in full quantum electrodynamics, and analyses that implicitly or explicitly assume that at any time at most a single photon is present invariably seem to come back to classical electrodynamics Prasad and Glauber (2000); Svidzinsky et al. (2010); Balik et al. (2013); T. Bienaimé and Kaiser (2013). Even if the literature about lattice systems of atoms interacting with light is thick with quantum systems, with a few notable exceptions (e.g., Tiecke et al. (2014); Douglas et al. (2015)) both the atoms and the light have so far mostly behaved essentially classically in both theories and experiments, and as such fit our framework.

The assumption about a 2D square lattice is benign; the same methodology works for any simple (Bravais) lattice. In the general case there will be three atomic resonances, and simplifications can be sought based on symmetries of the system. Lattices with an $n$-atom basis can be dealt with in an analogous way Perczel et al. (2017b), though the dimensions of the matrices would then be $3n\times 3n$. Possible resonance frequencies multiply accordingly, and the $3n$-dimensional eigenvectors need not be transparent. If the 2D Bravais lattices are not far enough apart, even for perpendicular incidence with $N$ lattices we have at least $N$ simultaneous equations to solve. That was, in fact, done for Figs. 4-6. If the angle of incidence is not perpendicular, the number of coupled linear equations necessarily triples, and for an $n$-atom basis the dimensionality of the linear algebra problem becomes $3nN$. And so on. These are straightforward technical complications, although the numerical effort can grow and the hopes for qualitative understanding can dwindle drastically.

An infinite lattice is not possible in practice, but it is a feature not a bug in our theory. One can easily solve light propagation in atomic systems in all relevant dimensions with direct numerical simulations. However, the demands on computational resources would usually dictate a limited number of the atoms, and a limited size of the sample. The optics of finite-size samples would then have to be disentangled from the simulation results Javanainen et al. (2017). We wish to unearth generic principles, and it would get difficult. In the same vein, we do not produce listings of quantitative results, but instead we have built a toolbox that anyone interested in specific questions in this area might find helpful.

Acknowledgements.

This work is supported in part by the National Science Foundation, Grant Nr. PHY-1401151.

Appendix A Numerical details

A.1 Richardson extrapolation

In the computations of the self-sum and the transfer sum we have to do two limits. We wish to have automated control of the limits, given a prescribed goal for precision. This is a fairly demanding problem in numerical analysis. Here we describe our solution.

By symmetry, the sum in (7) can be carried out in the first quadrant with only nonnegative integers, which saves 75% of the numerical effort.

Next, consider the sum in Eq. (7) for a given $\eta$. We would like to know how large a value of the upper limit $M=M(\eta)$ is needed to reach a prescribed relative error $\epsilon_{M}$ from the $M\rightarrow\infty$ limit of the sum. To this end, let us estimate the residual by replacing the discrete sum by an integral, using the leading $1/n$ dependence of the summand in the estimate. We find

[TABLE]

In the limit $M\rightarrow 0$ the integral does not represent the infinite sum particularly well for several obvious reasons, but it does give an idea of the scale of the sum

[TABLE]

The requirement that the relative error of the sum for the given upper limit $M$ equals $\epsilon_{M}$ gives the equation $|R_{S}(M)/R_{S}(0)|=\epsilon_{M}$, from which we may solve the upper limit $M$:

[TABLE]

We add another heuristic condition that the value $M(a,\eta)$ be at least 16.780913. The essentially random non-integer value is an attempt to reduce the probability that including or excluding a member $\bf n$ in the sum (7) depends on round-off errors.

One may study the precision by which the sum for any given upper limit $M$ approximates the infinite sum by making $M$ larger, and seeing how much the result changes. In this way we have found empirically that the expression (64) for the limit $M$ is numerically useful when one aims at a relative truncation error of the sum on the order of $\epsilon_{M}$. The computation time scales like $M^{2}$, and thus like $1/a^{2}$ for small $a$. We have devised various ad-hoc procedures to mitigate this problem, but most expediently we avoid small values of $a$ that apparently do not yield interesting results anyway.

The final limit $\eta\rightarrow 0+$ is more tricky. As may be seen from Eqs. (7) and (64), the number of the terms in the sum scales as $1/\eta^{2}$, and the numerical effort can get substantial in a brute-force attempt to make $\eta\rightarrow 0$.

We resort to Richardson extrapolation Press et al. (2007) instead. First we calculate a sequence of values for the sum, decreasing the convergence parameter $\eta$ by a factor $2$ every time:

[TABLE]

Assuming an expansion of the form

[TABLE]

for the sum as a function of the parameter $\eta$, we see that for small enough $\eta$ the ratio

[TABLE]

should be approximately 2 for all $n$. Empirically, it is so for all but some isolated values of $a$. So, next we form a sequence of numbers

[TABLE]

that has the property that the terms proportional to $\eta$ cancel, and the leading error in the sequence $S^{(1)}_{n}$ is proportional to $(\eta/2^{n})^{2}$. This would mean that the ratio as in Eq. (67) for the sequence $S_{n}^{(1)}$ should be approximately 4, which it empirically is. Continuing iteratively in this way, the leading error in the sequence

[TABLE]

should be proportional to $(\eta/2^{n})^{k+1}$. One can go in this way to higher orders $k$ to accelerate the convergence.

The main numerical issue is that the limit $\eta\rightarrow 0$ apparently genuinely diverges at the onset of new Bragg reflections, like $\eta^{-1/2}$ at $a=2\pi$. As far as we can tell, the sum converges for all values of $a$ that do not give a Bragg reflection in the plane of the lattice, but when $a$ approaches a new Bragg reflection, the values of the parameter $\eta$ for which the expansion (66) is useful become smaller and smaller. Correspondingly, Richardson extrapolation in the order $k$ may seem to converge, but gives a spurious result.

Our method to combat this problem runs as follows. We compute members of the sequence (65) one by one, typically starting with $\eta=0.08$ that has proven convenient for our purposes. Once we have three values, we do a Richardson extrapolation to the order that we have denoted by $k=2$, with an error proportional to $\eta^{3}$. We then keep adding terms to the sequence (65) and do the Richardson extrapolation with the last three members of the sequence to the order $k=2$, until the result is either deemed convergent or the process takes too much time and is terminated.

Specifically, if Richardson extrapolation works as advertised, the result with the leading error from using the last three terms of the sequence (65) up to the order $n$ is

[TABLE]

where $K$ is approximately a constant. Given this expression, we may then combine the last two extrapolation results into the formulas

[TABLE]

The former is in fact the Richardson extrapolation to the next order $k=3$ with the error $\propto\eta^{4}$, and the latter is essentially an estimate for the relative error of the extrapolation to the order $k=2$. We terminate the sequence (65) at the first sequence order $n$ such that, for the given relative error $\epsilon$, we have

[TABLE]

and declare the result (71). If convergence does not occur by some order $n$, typically 8, we call the result “error.” Of course, by increasing the limit order $n$, we could in principle force convergence whenever we are not at an exact Bragg diffraction value of $a$.

The final caveat is that the numbers $S^{(0)}_{n}$ obtained by truncating the sums at the limiting values $M(\eta)$ are already approximate, aiming at a relative error of $\epsilon_{M}$. If the noise in the sequence $S^{(0)}_{n}$ from this approximation is too large, Richardson extrapolation will not converge to the accuracy $\epsilon$. We safeguard against this failure by choosing $\epsilon_{M}=0.1\epsilon$. In our numerical surveys we typically use the relative-error parameter $\epsilon=10^{-3}$. Except very close to Bragg reflection points, this indeed gives relative errors smaller than or comparable to $\sim 10^{-3}$. Such errors are barely noticeable by eye in figures such as Fig. 3.

The same process may be used for the transfer sum $T(z,a)$. For the transfer sum the limit $\eta\rightarrow 0$ apparently also diverges at Bragg reflection points such as $a=2\pi$, and presents numerical problems nearby.

The basically same numerics also works for the tensor counterparts of the self-sum and transfer sum, $\sf S$ and $\sf T$, except that it deals with $3\times 3$ complex matrices not complex numbers. In the convergence criterion (73) we use the metric induced by the $\ell^{2}$ norm of a matrix so that, for instance, the “length” of the matrix $\sf S$ equals

[TABLE]

For numbers, $1\times 1$ matrices, this concept boils down to the absolute value as in (73), and in the matrix case it is not thrown off by very small entries in the matrices that occur because round-off errors cause nonzero values where the exact result would be zero.

We were able to run all of the results we needed with controlled precision as per above, on Mathematica for the scalar case and using mostly C++ for the matrix case. We can think of significant improvements to the numerical procedures, but the return would probably not justify the effort.

A.2 Transfer matrix

As to the long-distance analytical form of the transfer matrix, (41) and (45), we firstly without further ado replace the sum (36) with an integral over the site index ${\bf n}$, and express it in terms of polar coordinates in the plane. The integral over the angles is easy to carry out analytically using Mathematica. We choose the propagation vector in Eq. (34) with $\phi=0$, and are left with integrals such as

[TABLE]

where $J$ stands for the Bessel function. We do this integral numerically, again using Mathematica. For very small $\eta>0$ the value turns out the be independent of $z>0$, and by varying the angle $\theta$ we find that it equals $-\sin\theta$. This is the upper-right element of the matrix $M(\theta)$ in Eq. (41). The other elements come from similar arguments. Given the simple results it is obvious that a completely analytical procedure to find the elements of the matrix $M(\theta)$ should exist, but we have not bothered to find one. Finally, the matrix $R(\phi)$ of Eq. (45) represents rotation about the $z$ axis, and serves to extend the result to an arbitrary azimuthal angle $\phi$ under the evidently correct assumption that in such rotations of the direction of the incoming light the matrix $\sf T$ transforms as a tensor.

Appendix B Bragg scattering

Like in the main text, let us denote the wave vector of the incoming light by

[TABLE]

Its projection to the plane of the atoms is

[TABLE]

The incoming light paints a phase pattern on the atomic dipoles, which become proportional to $e^{i{\bf k}_{\parallel}\cdot{\bf R}_{\bf n}}$. The question is about the spatial pattern, especially in the far field, that the dipoles radiate in return. More specifically, there will be at least two obvious plane waves, the reflected field and the radiated part of the transmitted field. Our goal is to characterize the wave vectors $\bar{\bf k}$ for all possible radiated plane waves.

The principle that cracks the case is that a plane wave with the wave vector $\bar{\bf k}$, when incident on the atoms, would have to paint the same phase pattern $e^{i{\bf k}_{\parallel}\cdot{\bf R}_{\bf n}}$. Any vector $\bar{\bf k}$ whose projection to the plane satisfies $\bar{\bf k}_{\parallel}={\bf k}_{\parallel}$ fits the bill. However, not all of the $\bar{\bf k}$ with this property qualify. Namely, in steady state the frequency of the light is fixed, and therefore so is the absolute value of the vector $\bar{\bf k}$; $|\bar{\bf k}|=1$ in our units.

These considerations immediately leave two alternatives: a wave with $\bar{\bf k}={\bf k}$ that co-propagates with the incoming plane wave, and the reflected wave $\bar{\bf k}={\bf k}_{R}$ in which the $z$ component has flipped the sign. The transmitted wave is the superposition of the incoming field and the reradiated field with the wave number ${\bf k}$.

Nonetheless, we have not yet exhausted all possibilities. The dipoles reside at discrete positions, and any $\bar{\bf k}_{\parallel}$ that gives the same phase pattern as ${\bf k}_{\parallel}$ at the positions of the dipoles also qualify. This occurs precisely when these vectors differ by a reciprocal lattice vector,

[TABLE]

where $m_{x}$ and $m_{y}$ are arbitrary integers. But the radiation still has to occur at a wave vector whose absolute value equals $1$, so we have the condition ${\bar{\bf k}_{\parallel}}^{2}\leq 1$ to leave room for the $z$ components $\bar{k}_{z}=\pm\sqrt{1-{\bar{\bf k}_{\parallel}}^{2}}$. We have an inequality that may be rearranged to give

[TABLE]

Here $m_{x}=m_{y}=0$ always qualify, and give the forward and backward scattered waves. Any other combination of $m_{x}$ and $m_{y}$ that satisfy the inequality (79) signals the presence of a nontrivial Bragg-scattered wave. The condition for Bragg scattering depends on the azimuthal angle $\phi$. If this angle is varied freely, Bragg scattering for a given order $m_{x},\,m_{y}$ may be found for lattice spacings that satisfy

[TABLE]

Bragg scattering may thus occur for $a>\pi$, starting from a half-wavelength lattice spacing when the angle of incidence gets extremely oblique, $\theta\rightarrow\pi/2$. On the other hand, by having the incoming wave propagate at an angle $\phi=\pi/4$ so that the projection ${\bf k}_{\parallel}$ bisects the angle between lattice axes, one may delay the onset of Bragg scattering up to $a=\sqrt{2}\,\pi$.

2D Bragg scattering is also valid Bragg scattering in an evenly spaced 3D stack of 2D lattices. One may show easily that if there is no 2D Bragg scattering in a square lattice, there will be no nontrivial 3D Bragg scattering in a simple cubic lattice made by stacking 2D square lattices either. Other than warning that this may change if the spacing between the 2D lattices is larger than the lattice constant, we do not go into the details.

Appendix C The mean-field approximation

Our mean-field approximation for radiation in the lattice plane states that, if the fraction of lattice sites occupied is $\zeta$, on the average the field radiated by the other atoms on any given atom, the self-sum $S$, and the shift of the resonance $s=-\frac{3}{2}\Re[S]$, all get multiplied by $\zeta$. In this Appendix C we demonstrate with specific examples that the approximation is reasonable.

The mean-field approximation is not trivial because of multiple scattering and the ensuing cooperative phenomena: An atom radiates a field, which strikes another atom, which radiates a field, which comes back to the original atom, and so on. Here we study cooperative effects in the spirit of our earlier simulations of light propagation in an atomic sample Javanainen et al. (2014, 2017), by solving a set of linear equations that takes into account both the incoming light and the secondary fields from all atoms. The solution will give, among other things, the total field from all other atoms falling on any given atom.

In our examples we take the sample to be a circular cut from the lattice with the spacing $a=\pi$ (half of a wavelength). We include the sites that are at most the distance $R$ from the origin. Each site is occupied by an atom with the probability $\zeta$. We take a circularly polarized plane wave with the amplitude $E_{0}$ coming in on the lattice in a perpendicular direction, on resonance so that $\delta=0$.

If there are empty sites, even at the center of the sample the radiation from the atoms in general does not have the same circular polarization as the incoming light. However, on the average the polarization would be circular at the center, and in an infinite lattice at any site. We therefore project the electric field to the polarization of the incoming light, and characterize the secondary radiation from the other atoms on a given atom with a single complex amplitude. We compute this electric field $\bar{E}$ either at the center of the sample, or at a lattice site closest to the center that is occupied.

Modeling the secondary field as $\bar{E}=dS$, one can easily solve from Eqs. (10) and (11) the self-sum $S$, and obtain the shift of the resonance:

[TABLE]

We find it for a number of samples with the sites filled at random with the probability $\zeta$, and report the average as well as its standard deviation deduced from the samples. The fluctuations from sample to sample can be large, so the number of samples used in the averaging range up to a million. The computations are implemented by deriving suitable C++ classes from the ones used in our earlier 3D simulations, and the code works internally in 3D. Our standard sample radius $R=40\,a$ encloses $5025$ lattice sites, which means that the key operation is to solve up to a $15075\times 15075$ set of dense linear equations.

In the case with $R=40\,a$ and $\zeta=1$, all sites occupied, the result for the finite-size sample is $s=0.7958$. It should be compared with the result for an infinite lattice from Sec. II.2, $s=0.8006$. This level of the agreement, $0.6\%$, is probably fortuitous, but it lends credence to our methods.

Representative results are shown in Fig. 11, which plots the shift divided by the filling factor, ${s/\zeta}$, as a function of $\zeta$ over a range of two and a half orders of magnitude. The standard deviations are also given. Black filled circles are for the sample radius $R=40\,a$, and the field is computed at an occupied lattice site closest to the center. Data marked by red crosses are computed at the center of the circuclar cut of the lattice, whether there is an atom at the center or not. The points marked by blue open squares are for $R=80\,a$, and for an atom closest to the center.

In an infinite lattice, if it did not matter whether we are looking at a filled or non-filled site, and if the mean-field approximation were exact, the results would all come out at the same value of $s/\zeta$. The range of the values of $s/\zeta$ in the figure is on the order of a few tens of per cent. Variation of the magnitude we are seeing in Fig. 11 would not invalidate the conclusions of Sec. V. Besides, much of it evidently comes from factors other than the inaccuracy of the mean field approximation per se.

Appendix D Magnetic field and polarizability

In the present Appendix D we derive the polarizability tensor of an atom with a $J=0\rightarrow J=1$ transition in the presence of the magnetic field specified in terms of the polar angles $\theta$ and $\phi$ as

[TABLE]

classically. We think of the atom as a charged ($q$) damped ($\gamma$) harmonic oscillator (frequency $\omega_{0}$, mass $m$) under the influence of both the driving electric field {\bf E}(t)=\hbox{\frac{1}{2}}{\bf E}e^{-i\omega t}+\rm{c.c.} and the constant magnetic field B.

Newton’s second law gives the equation of motion for the position of the charge

[TABLE]

We look for the steady state with the Ansatz {\bf r}(t)=\hbox{\frac{1}{2}}{\bf r}_{0}e^{-i\omega t}. The true physical solution is ${\bf r}(t)+{\bf r}^{*}(t)$, so that ${\bf d}=q{\bf r}_{0}$ is the positive-frequency component of the dipole moment resulting from the driving field with the positive-frequency component ${\bf E}$. The dipole moment is a linear function of the electric field, ${\bf d}={\sf A}{\bf E}$, which defines the polarizability tensor $\sf A$.

The calculations are tedious, so we proceed with just an outline that we have filled in with Mathematica. Like in the undergraduate exercise of a driven damped harmonic oscillator, we find resonance denominators of the type $\omega^{2}-\omega_{0}^{2}$, yet we know from modeling of resonance behavior of atoms that it is the detuning $\Delta=\omega-\omega_{0}$ that matters most. We therefore develop a “resonance approximation” as follows:

Write the polarizability tensor in the form ${\sf A}={\cal A}/\omega$. Up to the point when all results are put together, we deal with the tensor $\cal A$. Compare this with the standard undergraduate exercise, in which one would make the resonance approximation by writing the resonance denominator as $\omega^{2}-\omega_{0}^{2}=(\omega-\omega_{0})(\omega+\omega_{0})\simeq 2\omega\Delta$. We are getting ahead of things here, but wish to point out that this scaling is particularly useful as in our custom units $\omega=1$.

• 2.

Write the tensor $\cal A$ in terms of $\omega_{0}$ and $\omega=\omega_{0}+\Delta$.

• 3.

Expand all components of the tensor $\cal A$ as partial fractions with the first power of $\Delta$ in the denominators, zeroth power in the numerators.

• 4.

Replace each partial fraction with the leading term in its expansion in the limit $\omega_{0}\rightarrow\infty$.

• 5.

Keep only the partial fractions that are of the dominant, zeroth, order in $\omega_{0}$.

This approximation produces terms with resonance denominators of the form $\Delta+i\gamma$ and $\Delta\pm\Omega_{B}+i\gamma$, where $\Omega_{B}={qB}/{2m}$ is half of the cyclotron frequency for the given magnetic field $B$. Finally, we compare these results in the case of zero magnetic field with the expression (3) for scalar polarizability, and write the corresponding expressions with the resonance denominators that include the magnetic field. The polarizability tensor is unwieldy, but here is the $xx$ component as an example:

[TABLE]

with $\omega_{B}=\Omega_{B}/\gamma$.

While the derivation we have outlined is classical and seems ad-hoc, the results agree with fully quantum mechanical numerical calculations as reported by the “Software Atom” Javanainen (2017) for the $J=0$ to $J^{\prime}=1$ transition. In the quantum case $\omega_{B}$ equals the frequency difference caused by the magnetic field between the states $m^{\prime}=1$ and $m^{\prime}=0$ of the excited-state manifold, expressed in units of the linewidth of the transition, and is negative if the Lande factor for the level $J^{\prime}=1$ is positive. Another point in favor of our version of the resonance approximation is that energy is conserved: We have checked with an explicit calculation that, in steady state and over a period of the oscillations, the energy that the driving electric field puts on the charge equals the energy that the oscillating dipole radiates.