Super- and subradiance of clock atoms in multimode optical waveguides

Over the last decades, high precision spectroscopy as one of the leading fields in quantum metrology [1] has brought forward ground-breaking results such as the world's best atomic clocks [2–4] alongside all kinds of state-of-the-art devices as magnetometers, accelerometers or other quantum sensing devices [1]. To improve measurements one seeks to minimize quantum projection noise diminishing as $1/\sqrt{N}$ with an increasing number N of emitters. However, when many atoms are confined in a small volume this leads to a compromise between reducing projection noise and introducing detrimental interaction effects owing to larger atom densities and the resulting increased dipole–dipole couplings. One close to ideal approach for clocks is to place the atoms into a magic wavelength optical lattice, cancelling the differential lightshifts between the ground and excited levels, with one particle per site [5], where interactions occur via the remaining resonant dipole–dipole coupling only [6]. This leads to unprecedented precision and accuracy at the cost of a high preparation effort.

The platform of optical fibers has attracted a lot of interest in the last years with tapered nano-fibers used for trapping and observation of single emitters from cold atoms [7, 8] to nanoparticles [9] with applications as a single particle optical switch [10] or an optical isolator [11]. Similarly, atoms have been trapped inside hollow core photonic crystal fibers for spectroscopy [12–15] or nonlinear optics [16] and even single molecules were coupled to dielectric nano- waveguides [17]. In parallel, a great deal of theoretical investigations of the effective collective atom-field dynamics, the resulting light forces and possible applications in these systems have been carried out ranging from self organization [18–21], to nonlinear scattering [22, 23]. Applications of the collective radiative behavior into the light modes also point towards improved photonic memories [24]. Other recent studies include [17, 25–27].

Along these lines, the idea to use optical fibers and in particular hollow core photonic crystal fibers [28] for high precision spectroscopy was put forward by Okaba et al in 2014 [13]. This should allow for a collision-free setup of a large number of atoms trapped in a one-dimensional optical lattice with a filling factor of ≤1, which can be addressed in a uniform way via fiber guided modes. Here, confinement comes at the expense of some remaining atom-wall interactions and the sometimes unwanted enhanced long range dipole–dipole coupling. In their study Okaba et al saw clear indications of collective radiative emission as well as energy shifts and studied the effect of tight confinement and atomic interactions with the fiber wall. Nevertheless, even in their first pioneering study they were already able to reduce the relative uncertainty in measuring the ${}^{1}{S}_{0}{\to }^{3}{P}_{1}$ transition in ⁸⁸Sr to the order of 10⁻¹³.

For long lived atomic transitions as used in optical clocks, free space spontaneous emission is only a small hinderance to precision measurements. However, at higher densities (leading to a larger optical depth) even in free space scenarios, collective radiative effects can be increased drastically [29–31]. Such a fast collective decay (dubbed superradiance) can substantially decrease the sensitivity of a Ramsey type measurement. However, collective subradiant states coexist with the superradiant ones and the choice of favorable lattice and excitation geometries can lead to improved precision via subradiant state protection [32, 33].

In this paper we study superradiant/subradiant effects occurring in a 1D optical fiber geometry, where in addition to the collective coupling to the 3D free space vacuum modes, the guided discrete modes of the fiber have the potential to introduce an additional virtually infinite-range dipole–dipole coupling. This is particularly true in single transverse mode geometries [18]. Recently, closely related long range superradiance effects were observed in evanescent wave nano-fibers [34]. Moreover, we generalize the treatment to multimode fibers with a particular emphasis on photonic bandgap hollow core geometries as used in [13]. We consider a regular 1D chain of clock atoms coupled to several modes of a multimode optical fiber including free space spontaneous emission with inherent dipole–dipole coupling. As a reference point towards new physics we first study the enhanced spontaneous emission for a single mode fiber using various numerical approaches from the full master equation (ME) down to a simple mean field (MF) approach. In a further step we generalize these ideas to several fiber modes. While for small systems sizes the ME can be solved directly, we have to resort to enhanced mean-field methods for larger particle numbers.

We consider N identical two-level emitters with transition frequency ω₀ positioned at points $\{{\vec{r}}_{i}{\}}_{i=1}^{N}$ along a chain inside an (infinitely extended) optical fiber (see figure 1) of cross-section A = πa² (where a is an effective fiber radius). The atoms interact both with the free space modes as well as with the fiber guided modes. We will work in the approximation that the two processes can be separated, which is valid as long as the non-guided fiber modes are not that substantially different from the free-space modes, as is the case in the experiment by Okaba et al [13], where the fiber's diameter is large compared to the interatomic distance. In situations, where the unguided local field in the fiber needs to be evaluated explicitly, a Green's tensor formalism could be employed as described in [35].

Zoom In Zoom Out Reset image size

After the elimination of the radiation field, this results in dipole–dipole contributions stemming from the free-space modes (${{\rm{\Omega }}}_{{ij}}^{3{\rm{D}}}$) and guided modes (${{\rm{\Omega }}}_{{ij}}^{1{\rm{D}}}$). Additionally, incoherent mutual processes resulting from free space vacuum fluctuations are quantified by ${{\rm{\Gamma }}}_{{ij}}^{3{\rm{D}}}$ while the guided modes lead to mutual decay rates ${{\rm{\Gamma }}}_{{ij}}^{1{\rm{D}}}$.

2.1. Effective ME

The full Hamiltonian for the system of atoms and radiation modes can be decomposed into contributions originating from coupling to 3D free space and 1D fiber modes and we have

$$\begin{eqnarray}&&{H}_{\mathrm{full}}={H}_{{\rm{A}}}+{H}_{{\rm{F}}}^{3{\rm{D}}}+{H}_{{\rm{F}}}^{1{\rm{D}}}+{H}_{\mathrm{int}}^{3{\rm{D}}}+{H}_{\mathrm{int}}^{1{\rm{D}}},\end{eqnarray} \tag{ 1 }$$

with ${H}_{{\rm{A}}}={\omega }_{0}{\sum }_{i}{\sigma }_{i}^{+}{\sigma }_{i}^{-}$, where ${\sigma }_{i}^{\pm }$ denote the raising and lowering operators of the ith emitter and we assume identical transition frequencies ${\omega }_{0}={{ck}}_{0}=2\pi c/{\lambda }_{0}$ for all atoms. The free Hamiltonian of the radiation modes is split as ${H}_{{\rm{F}}}^{3{\rm{D}}}={\sum }_{\vec{k},\lambda }{\omega }_{k}{a}_{\vec{k},\lambda }^{\dagger }{a}_{\vec{k},\lambda }$ with the frequency $\omega =c\left|\vec{k}\right|$, the wave vector $\vec{k}$ and the polarization λ, and ${H}_{{\rm{F}}}^{1{\rm{D}}}={\sum }_{\eta }\omega {b}_{\eta }^{\dagger }{b}_{\eta }$ where $\eta =\left(\omega ,\lambda ,\nu ,f\right)$ with the frequency ω, polarization λ, the mode index ν and the propagation direction f = ±1. We have assumed complete independence of the free-space and guided modes such that $\left[{a}_{\vec{k},\lambda }^{\dagger },{b}_{\eta }\right]=\left[{a}_{\vec{k},\lambda },{b}_{\eta }^{\dagger }\right]=0$. Furthermore, we have $\left[{a}_{\vec{k},\lambda },{a}_{{\vec{k}}^{{\prime} },{\lambda }^{{\prime} }}^{\dagger }\right]={\delta }_{\vec{k}{\vec{k}}^{{\prime} }}{\delta }_{\lambda {\lambda }^{{\prime} }}$ and $\left[{b}_{\eta },{b}_{{\eta }^{{\prime} }}^{\dagger }\right]={\delta }_{\eta {\eta }^{{\prime} }}$. The two interaction terms in dipole and rotating wave approximation are given in the appendix.

The elimination of the radiation modes is done in the standard quantum optics procedure [36] by formally integrating the equations of motion for ${a}_{\vec{k},\lambda }$ and b_η and inserting the solution into an integro-differential equation for an arbitrary operator of the emitters. Under the Markov approximation (see [37, 38]) one can show that the dynamics of the reduced system of N atoms is properly described by the effective ME

$$\begin{eqnarray}&&\dot{\rho }={\rm{i}}\left[\rho ,H\right]+{ \mathcal L }[\rho ],\end{eqnarray} \tag{ 2 }$$

where the eliminated fields result in an effective dipole–dipole interaction via free space as well as via the guided modes, i.e.

$$\begin{eqnarray}&&H={\omega }_{0}\displaystyle \sum _{i=1}^{N}{\sigma }_{i}^{+}{\sigma }_{i}^{-}+\displaystyle \sum _{i\ne j}\left({{\rm{\Omega }}}_{{ij}}^{3{\rm{D}}}+{{\rm{\Omega }}}_{{ij}}^{1D}\right){\sigma }_{i}^{+}{\sigma }_{j}^{-}.\end{eqnarray} \tag{ 3 }$$

The two contributions are fundamentally different. The free space coupling is finite range with considerable strengths only around the distances of the order of a wavelength

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{ij}}^{3{\rm{D}}}\left(\xi ={k}_{0}{r}_{{ij}}\right)=-\displaystyle \frac{3{\rm{\Gamma }}}{4}\left[\displaystyle \frac{\cos \xi }{\xi }-\displaystyle \frac{\sin \xi }{{\xi }^{2}}+\displaystyle \frac{\cos \xi }{{\xi }^{3}}\right],\end{eqnarray} \tag{ 4a }$$

while the guided modes yield a sum of infinite range interactions, periodic along the fiber's axis (as in [39])

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{ij}}^{1{\rm{D}}}\left(\xi ={k}_{0}{z}_{{ij}}\right)=\displaystyle \frac{{{\rm{\Gamma }}}^{1{\rm{D}}}}{2}\ \displaystyle \sum _{\nu }{\chi }^{\nu }\sin \left({\beta }^{\nu }\xi \right),\end{eqnarray} \tag{ 4b }$$

where χ^ν expresses the coupling strength of each mode and β^ν is the propagation constant as defined in equation (A.16) in the appendix discussing the 1D field contributions. The 3D field-induced couplings are derived in [37].

The coherent part is accompanied by incoherent processes responsible for dissipative dynamics such as super- and subradiance. We have

$$\begin{eqnarray}&&{ \mathcal L }\left[\rho \right]=\displaystyle \frac{1}{2}\displaystyle \sum _{i,j}\left({{\rm{\Gamma }}}_{{ij}}^{3{\rm{D}}}+{{\rm{\Gamma }}}_{{ij}}^{1{\rm{D}}}\right)\left[2{\sigma }_{i}^{-}\rho {\sigma }_{j}^{+}-{\sigma }_{i}^{+}{\sigma }_{j}^{-}\rho -\rho {\sigma }_{i}^{+}{\sigma }_{j}^{-}\right].\end{eqnarray} \tag{ 5 }$$

As above, a finite range free-space contribution,

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{ij}}^{3{\rm{D}}}\left(\xi ={k}_{0}{r}_{{ij}}\right)=\displaystyle \frac{3{\rm{\Gamma }}}{2}\left[\displaystyle \frac{\sin \xi }{\xi }+\displaystyle \frac{\cos \xi }{{\xi }^{2}}-\displaystyle \frac{\sin \xi }{{\xi }^{3}}\right],\end{eqnarray} \tag{ 6a }$$

competes with the sum of infinite range fiber mediated decay rates,

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{ij}}^{1{\rm{D}}}\left(\xi ={k}_{0}{z}_{{ij}}\right)={{\rm{\Gamma }}}^{1{\rm{D}}}\ \displaystyle \sum _{\nu }{\chi }^{\nu }\cos \left({\beta }^{\nu }\xi \right).\end{eqnarray} \tag{ 6b }$$

The single particle emission rates are

$$\begin{eqnarray}&&{\rm{\Gamma }}=\displaystyle \frac{{k}_{0}^{3}{\mu }^{2}}{3\pi {\epsilon }_{0}}\qquad {{\rm{\Gamma }}}^{1{\rm{D}}}=\displaystyle \frac{{k}_{0}{\mu }^{2}}{{\epsilon }_{0}A}=\alpha {\rm{\Gamma }},\end{eqnarray} \tag{ 7 }$$

where $\alpha := {{\rm{\Gamma }}}^{1{\rm{D}}}/{\rm{\Gamma }}=3\pi /\left({k}_{0}^{2}A\right)$ and $\mu =\left|\vec{\mu }\right|$ is the atomic transition dipole moment, which we assume to be orthogonal to the propagation direction in the fiber, i.e. $\vec{\mu }\parallel {\vec{e}}_{x}$. The prefactor in the sum over all the guided modes of the 1D confined field reads

$$\begin{eqnarray}&&{\chi }^{\nu }={\left.\displaystyle \frac{{\rm{d}}{\beta }^{\nu }}{{\rm{d}}{k}}\right|}_{k={k}_{0}}\cdot \displaystyle \frac{A}{{C}^{\nu }}{\left|{E}_{x}^{\nu }\right|}^{2}\end{eqnarray} \tag{ 8 }$$

and incorporates the field E_x experienced by the dipoles, the normalization of the field ${\int }_{0}^{\infty }{\rm{d}}r{\int }_{0}^{2\pi }{\rm{d}}\varphi {\left|{{nE}}^{\nu }\right|}^{2}={C}^{\nu }$ for each mode and the (scaled) group velocity. The field itself is obtained by solving Maxwell's equations in cylindrical coordinates for a cylindrical fiber of cross section A = πa² and refractive index ${n}_{1}\gt 1$ cladded by vacuum n₂ = 1 (see appendix) in a similar manner to [40].

2.2. Numerical methods

The presented results mainly stem from numerical investigations based on the reduced ME of the N atoms under different initial conditions. As full simulations in the Hilbert space of N two-level systems are computationally costly, we make use of different approximation methods factorizing in terms of individual operators or pairs of operator correlations. In increasing orders of computational time, the methods consist of: (i) a MF model, (ii) an augmented MF model developed in [41] including particle–particle correlations, which we dub 'MF plus correlations' (MPC) and (iii) an exact integration of the ME as a basis of comparison [42]. For the exact ME solution the number of equations that need to be integrated scales as 4^N. The MF approach reduces the effort in the integration to 3N equations, i.e. to the expectation values  $\left\langle {\sigma }_{j}^{x}\right\rangle$, $\left\langle {\sigma }_{j}^{y}\right\rangle$ and $\left\langle {\sigma }_{j}^{z}\right\rangle$ for every atom (j = 1, ...N). In the MPC approach we include second order correlations and therefore our number of equations will scale as $3N\left(2N+1\right)$. More precisely, for MF we choose

$$\begin{eqnarray}&&\rho =\mathop{\displaystyle \bigotimes }\limits_{k}{\rho }^{(k)},\end{eqnarray} \tag{ 9 }$$

which effectively factorizes all correlations of the form $\left\langle {AB}\right\rangle \approx \left\langle A\right\rangle \left\langle B\right\rangle$, while the second order MPC approach is based upon the ansatz

$$\begin{eqnarray}&&\rho =\mathop{\displaystyle \bigotimes }\limits_{i}{\rho }^{(i)}+\displaystyle \sum _{j\lt k}\left({\rho }^{(j,k)}\otimes \mathop{\displaystyle \bigotimes }\limits_{i\ne j,k}{\rho }^{(i)}\right)\end{eqnarray} \tag{ 10 }$$

and we factorize $\left\langle {ABC}\right\rangle \approx -2\left\langle A\right\rangle \left\langle B\right\rangle \left\langle C\right\rangle +\left\langle A\right\rangle \left\langle {BC}\right\rangle +\left\langle B\right\rangle \left\langle {AC}\right\rangle +\left\langle C\right\rangle \left\langle {AB}\right\rangle$. In order to investigate the validity of our approximations we compare the two inexact methods with the solution of the full ME for reduced ensembles of not more than 10 atoms. Where the dynamics are well contained within the approximations, we can perform numerical simulations up to a few hundred atoms (for MPC) or up to 10⁴ atoms (for MF).

2.3. Initial states selected for the collective dynamics

Setups where a suitable choice of the inter-particle separation can lead to the cancellation of differential light shifts, as described e.g. in [43, 44], are of particular interest. We therefore focus on such a magic wavelength configuration with d/λ₀ ≈ 0.59, as used in optical lattice experiments with trapped Strontium atoms driven on the ¹S₀ to ³P₀ transition [13]. These distances feature relatively minute dipole–dipole energy shifts as well. We then fix the effective fiber propagation constant to β = 1.2 and set the ratio α (largely defined by the fiber's diameter, which can be modified to tailor the relative role of the free-space and confined field collective interactions) to 0.75. The dynamics is characterized by following the time evolution of the total excitation operator

$$\begin{eqnarray}&&{S}^{z}=\displaystyle \sum _{j}{\sigma }_{j}^{+}{\sigma }_{j}^{-},\end{eqnarray} \tag{ 11 }$$

from the initial value of N.

Generally, we find quantitative agreement between the exact ME solution and the MPC solution for short times and at least a qualitative one for longer times (see figure 2(a)). However, while two particle correlations between all pairs of emitters, long- and short-range for guided as well as free space modes, as included in the MPC method, are sufficient to capture the initial superradiant behavior, larger time subradiance implies more than two particle correlations. In contrast, the simple MF result is often far off and usually only useful for independently decaying ensembles; this is because the MF cannot capture collective effects such as entanglement or multipartite correlations. Especially when starting from a perfectly inverted state, the MF approximation fails to correctly account for the buildup of correlations and coherence. The MF equations are setup in such a way that they feature a driving term for the individual expectation values $\left\langle {\sigma }_{z}^{i}\right\rangle$ consisting of a sum over the other dipole's expectation values $\left\langle {\sigma }_{x}^{j}\right\rangle$ and $\left\langle {\sigma }_{y}^{j}\right\rangle$. When starting from a perfectly inverted state the latter two are zero for all times and thus the equation of motion for $\left\langle {\sigma }_{z}^{i}\right\rangle$ will result in uncoupled dynamics, as no seed for the build-up of superradiance is present.

Zoom In Zoom Out Reset image size

In the MPC treatment a perfectly inverted state will feature three-particle correlations of the form $\left\langle {\sigma }_{z}^{i}{\sigma }_{z}^{j}{\sigma }_{z}^{k}\right\rangle$, which, when factorized in to two-particle correlations and single-particle expectation values, can provide for this initial seed and thus result in superradiant dynamics for a perfectly inverted initial state. The full set of equations can be found in [41].

However, assuming an initial macroscopic non-zero dipole (which is normally associated with the initial spontaneous build-up of macroscopic coherence in superradiance), the MF approximation is closer to the real dynamics. We can achieve this regime by assuming imperfect excitation pulses close but not exactly at the π value. This could be the result of non-negligible dipole–dipole interactions and collective decay during the excitation.

Numerically we simulate the excitation under the following Hamiltonian (in an interaction picture removing the laser frequency) for the duration τ of the pulse

$$\begin{eqnarray}&&{H}_{\mathrm{exc}}=H+{H}_{{\rm{L}}}=\displaystyle \sum _{i\ne j}{{\rm{\Omega }}}_{{ij}}{\sigma }_{i}^{+}{\sigma }_{j}^{-}+{{\rm{\Omega }}}_{{\rm{R}}}\displaystyle \sum _{i}({\sigma }_{i}^{+}+{\sigma }_{i}^{-}),\end{eqnarray} \tag{ 12 }$$

supplemented with the collective Lindblad term describing decay. We have assumed resonant laser excitation with the Rabi frequency ${{\rm{\Omega }}}_{{\rm{R}}}$. For an intended π pulse we set the driving time τ by demanding that ${{\rm{\Omega }}}_{{\rm{R}}}\tau =\pi /2$. The interatomic separation is fixed at the magic wavelength trapping distance d/λ₀ = 0.59. Note, that at this distance, the one-dimensional collective decay (of infinite range) due to the fiber dominates over the free space collective decay (falling off with increasing distance). For the dynamics in figure 2(b), we consider a square pulse with Rabi frequency ${{\rm{\Omega }}}_{{\rm{R}}}=100{{\rm{\Gamma }}}^{1{\rm{D}}}$, for which the population of the excited state corresponds to a rotation with an angle θ = 0.95π on the collective Bloch sphere. Preserving some finite coherence in the spins as opposed to a perfect inversion. We compare independent decay for 1000 spins prepared in the perfectly inverted state with the situation where the ensemble is prepared as described above via the application of a excitation laser. The perfect π pulse leads to initial conditions $\langle {\sigma }_{j}^{x}\rangle (t=0)=0$ and $\langle {\sigma }_{j}^{y}\rangle (t=0)=0$ which evolve into $\langle {\sigma }_{j}^{x}\rangle (t)=\langle {\sigma }_{j}^{y}\rangle (t)=0$ for all times. This implies that the atoms are decaying independently. In reality, spontaneous emission events lead to the emergence of a collective macroscopic dipole which speeds up the collective decay. As figure 2(b) shows, this can be simulated by assuming the imperfect preparation scheme which creates an initial macroscopic dipole leading to faster decay than that of an independent ensemble.

Let us now analyze the effect of the competition between free radiation modes induced decay versus collective decay aided by guided modes. For increasing values of the ratio α, figure 3 quantifies the additional superradiance induced by the fiber. The parameters chosen here are: ${\lambda }_{0}=689\,\mathrm{nm}$, χ¹ = 1 and the propagation constant is β¹ = 1.2. While the free space superradiance is practically negligible at magic wavelength separations, the long-range fiber mediated interactions, even for modest α, can lead to considerable superradiant behavior. For example, for the parameters of [13], where λ₀ = 689 nm and a = 20 μm, the resulting ratio $\alpha \sim {10}^{-4}$, requires a large number of atoms only, i.e. more than 10⁴, to show a sizable effect of the guided modes' superradiance. This aspect is more clearly analyzed in figure 3(b) with the direct conclusion that the more emitters are present in the ensemble, the more dominant the 1D decay is over the 3D free space decay. The figure shows the time at which half the population has decayed as a function of the number of emitters and the ratio between 1D and 3D decay when starting from the fully inverted state.

Zoom In Zoom Out Reset image size

In standard Ramsey interferometry, ensembles of atoms are prepared in states exhibiting maximal dipole moments along some direction, typically via the application of π/2 pulses. As described in section 2.3, at least two distinct procedures can be distinguished: (i) transversal versus (ii) longitudinal excitation. While transversal excitation imprints the same phase on all atoms, longitudinal excitation imprints different phases along the preferential direction defined by the fiber axis. We perform a detailed analysis of the numerical methods used for the two cases which we illustrate in figure 4. For both situations, as concluded before, the MPC solution is a good estimate of the real dynamics while the MF solution depends strongly on the particularities of the initial state. In the transverse case, as shown in figure 4, the dynamics are mostly subradiant. The effects largely stem from the fiber-induced collective decay.

Zoom In Zoom Out Reset image size

Starting in the transversally pumped Ramsey state, which features a maximal dipole moment, the MF solution is much closer to the exact result but will still not completely concur with the ME results on a long time scale. This initial state has a substantial overlap with states not radiating into the fiber which adds a long-lived component in the inversion that is again captured by the MPC corrections. In the other extremal case of a fiber excited Ramsey state, where we have an initial non-zero value for the Pauli operator in the x-direction all models capture superradiance along the fiber quite well, in contrast to the independent decay. Therefore, a MF treatment represents a suitable choice to capture the system's dynamics when compared to the full ME solution. We observe that in the experimentally most relevant case of a Ramsey sequence driven through the fiber, superradiance is strongly present limiting measurement time and thus precision.

The diameter of optical fibers used in current experiments in conjunction with quantum emitters can range from ∼0.5 to ∼40 μm. While for sub-micron fibers the single mode approximation is well justified, for large diameter fibers one necessarily needs to consider more than one allowed mode. Of course, in the limit of very large diameters one would end up in the free space regime. This is equally true for evanescent fields of nanofibers as well as for the guided field inside a hollow core fiber. As the transverse mode functions are orthogonal, the Lindblad decay terms are independent and additional modes can be incorporated into our theory directly. Hence, in analogy to equations (4b) and (6b) the dipole–dipole interaction as well as the collective dissipation through the 1D guided modes is given by

$$\begin{eqnarray*}&&{{\rm{\Omega }}}_{{ij}}^{1{\rm{D}}}=\displaystyle \frac{{{\rm{\Gamma }}}^{1{\rm{D}}}}{2}\ \displaystyle \sum _{\nu }{\chi }^{\nu }\sin \left({\beta }^{\nu }\xi \right)\qquad {{\rm{\Gamma }}}_{{ij}}^{1{\rm{D}}}={{\rm{\Gamma }}}^{1{\rm{D}}}\ \displaystyle \sum _{\nu }{\chi }^{\nu }\cos \left({\beta }^{\nu }\xi \right),\end{eqnarray*}$$

where ν is again the mode index and the propagation parameter β^ν and χ^ν are particular to every mode and have to be determined numerically as a function of the fiber properties (see appendix).

For simplicity, we choose a model of an infinitely extended cylindrical fiber where analytical expressions exist; in principle, one can use any suitable model calculation (e.g. using programs as COMSOL) in order to obtain the field in arbitrary optical structures, such as hollow core photonic crystals, tapered waveguides or nano fibers as done in [45]. The only thing that will change with respect to our formalism are the values of χ^ν and β^ν. As the collective couplings are a superposition of trigonometric functions of different periodicity since β^ν is different for each mode, by using different fiber diameters it is thus possible to tailor the interaction range. This is illustrated in figure 5, where we observe regions of more and less pronounced collective couplings. Consequently, instead of treating the multimode problem as a detrimental effect, we propose here that one can carefully include a number of modes to design interactions of tailorable range. The quite pronounced exponential decay of the couplings with distance can be achieved by allowing each mode to have a certain (Lorentzian) bandwidth within the interval of allowed propagation constants as opposed to only one discrete contribution. The bandwidths are chosen in such a way that there is an overlap between the modes of a few percent [21].

Zoom In Zoom Out Reset image size

For a reliable prediction of collective dynamics of quantum emitters periodically arranged within the field of a fiber, one has to include infinite range interactions via the effective 1D geometry of the system as well as finite range contributions mediated by the free 3D electromagnetic vacuum. In principle, no matter how small the atom-fiber coupling is, the 1D contribution will always dominate and introduce superradiant decay for large enough atom numbers. The combination of both effects introduces extra line shifts, which are small but important for clock applications. Interestingly, shifts and collective decay are reduced in multimode fibers, where only atoms within a finite range will participate so that divergences of shifts and decay rates in the limit of large particle numbers are prevented. In some regimes where the build-up of particle–particle correlations can be neglected, MF calculations can give reasonable quantitative estimates of the collective system dynamics. However, two-particle correlations can quite reliably be accounted for in the MPC method derived in [41] and allow for a general simulation of the correct dynamics for hundreds of particles. Yet, some aspects of subradiance seem to necessarily include more than particle–particle quantum correlations falling outside the validity of the MPC approach. As a bottom line, a hollow core multimode fiber system based on a fiber guided magic lattice trap could thus be a favorable way for a technically simplified but still accurate and precise implementation of an optical clock. It requires, however, a careful design of the excitation process and lattice geometry to keep collective shifts and superradiance limited. In this respect a transformation to an active reference oscillator system looks very promising, if a clever way of pumping and collective cooling [46] via higher order fiber modes could be engineered.

Calculations were performed using the Quantum optics and Collective Spins Julia packages developed by Sebastian Krämer [47]. We thank M Moreno-Cardoner, A Asenjo-Garcia and D Plankensteiner for insightful discussions. We acknowledge support by the Austrian Science Fund via project P29318-N27 (LO) and the SFB FoQus project F4013 (H R). CG thanks the University of Innsbruck for hospitality and the Max Planck Society for financial support.

Our system consists of a cylindrical infinitely long vacuum clad fiber of radius a and refractive index n₁ embedded in vacuum (n₂ = 1) with N two-level atoms (ω₀, μ) fixed at positions $\{{z}_{j}{\}}_{j=1}^{N}$ along the fiber.

A.1. Electric field in and outside the fiber

For the guided field in the fiber we solve Maxwell's equations, which we sketch briefly. For the electric field along the fiber axis we have the wave equation in radial coordinates with the propagation parameter β

$$\begin{eqnarray}&&\square {E}_{z}=0.\end{eqnarray} \tag{ A.1 }$$

With the ansatz

$$\begin{eqnarray}&&{E}_{z}={Af}(r){{\rm{e}}}^{{\rm{i}}\nu \varphi }\cdot {{\rm{e}}}^{{\rm{i}}\left(\omega t-\beta z\right)}\end{eqnarray} \tag{ A.2 }$$

this reduces to

$$\begin{eqnarray}&&\left({\partial }_{r}^{2}+\displaystyle \frac{1}{r}{\partial }_{r}+\displaystyle \frac{1}{{r}^{2}}{\partial }_{\varphi }^{2}+\left({n}^{2}{k}^{2}-{\beta }^{2}\right)\right){E}_{z}=0\end{eqnarray} \tag{ A.3 }$$

with the refractive index n (n₁, n₂) and ${k}^{2}={\omega }^{2}/{c}^{2}={\epsilon }_{0}{\mu }_{0}{\omega }^{2}$. The same equation holds for the magnetic field H_z and from that we can calculate the other two electric field components as

$$\begin{eqnarray}&&{E}_{r}=\displaystyle \frac{-{\rm{i}}}{{\kappa }^{2}}\left(\beta \,{\partial }_{r}{E}_{z}+\displaystyle \frac{\omega {\mu }_{0}}{r}\,{\partial }_{\varphi }{H}_{z}\right),\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}&&{E}_{\varphi }=\displaystyle \frac{-{\rm{i}}}{{\kappa }^{2}}\left(\displaystyle \frac{\beta }{r}\,{\partial }_{\varphi }{E}_{z}-\omega {\mu }_{0}\,{\partial }_{r}{H}_{z}\right).\end{eqnarray} \tag{ A.5 }$$

For the radial component f(r) we obtain the differential equation for the Bessel functions

$$\begin{eqnarray}&&{\partial }_{r}^{2}f+\displaystyle \frac{1}{r}{\partial }_{r}f+\left({\kappa }^{2}-\displaystyle \frac{{\nu }^{2}}{{r}^{2}}\right)f=0.\end{eqnarray} \tag{ A.6 }$$

This results in the fields

$$\begin{eqnarray}f(r)=\left\{\begin{array}{ll}{{AJ}}_{\nu }({hr}) & 0\lt r\lt a\\ {{CK}}_{\nu }({qr}) & a\lt r\lt \infty \end{array}\right.,\end{eqnarray} \tag{ A.7 }$$

where $h=\sqrt{{n}_{1}^{2}{k}^{2}-{\beta }^{2}}$ and $q=\sqrt{{\beta }^{2}-{n}_{2}^{2}{k}^{2}}$. J_ν denotes the Bessel function of the first kind, while K_ν is the modified Bessel function of the second kind. Similar solutions exist for H_z (constants B and D). Thus, we can construct the electric fields for $r\lt a$

$$\begin{eqnarray}&&{E}_{z}={{AJ}}_{\nu }({hr}){{\rm{e}}}^{{\rm{i}}\nu \varphi },\end{eqnarray} \tag{ A.8 }$$

$$\begin{eqnarray}&&{E}_{r}=\displaystyle \frac{-{\rm{i}}}{{h}^{2}}\left(\beta {{hAJ}}_{\nu }^{{\prime} }({hr})+{\rm{i}}\omega {\mu }_{0}\displaystyle \frac{\nu }{r}{{BJ}}_{\nu }({hr})\right){{\rm{e}}}^{{\rm{i}}\nu \varphi },\end{eqnarray} \tag{ A.9 }$$

$$\begin{eqnarray}&&{E}_{\varphi }=\displaystyle \frac{-{\rm{i}}}{{h}^{2}}\left({\rm{i}}\beta \displaystyle \frac{\nu }{r}{{AJ}}_{\nu }({hr})-h\omega {\mu }_{0}{{BJ}}_{\nu }^{{\prime} }({hr})\right){{\rm{e}}}^{{\rm{i}}\nu \varphi }\end{eqnarray} \tag{ A.10 }$$

and for r > a

$$\begin{eqnarray}&&{E}_{z}={{CK}}_{\nu }({qr}){{\rm{e}}}^{{\rm{i}}\nu \varphi },\end{eqnarray} \tag{ A.11 }$$

$$\begin{eqnarray}&&{E}_{r}=\displaystyle \frac{-{\rm{i}}}{{q}^{2}}\left(\beta {{qCK}}_{\nu }^{{\prime} }({qr})+{\rm{i}}\omega {\mu }_{0}\displaystyle \frac{\nu }{r}{{DK}}_{\nu }({qr})\right){{\rm{e}}}^{{\rm{i}}\nu \varphi },\end{eqnarray} \tag{ A.12 }$$

$$\begin{eqnarray}&&{E}_{\varphi }=\displaystyle \frac{-{\rm{i}}}{{q}^{2}}\left({\rm{i}}\beta \displaystyle \frac{\nu }{r}{{CK}}_{\nu }({qr})-q\omega {\mu }_{0}{{DK}}_{\nu }^{{\prime} }({qr})\right){{\rm{e}}}^{{\rm{i}}\nu \varphi }.\end{eqnarray} \tag{ A.13 }$$

From the boundary condition for the tangential field components at r = a we find

$$\begin{eqnarray}&&C=\displaystyle \frac{{J}_{\nu }({ha})}{{K}_{\nu }({qa})}A\qquad D=\displaystyle \frac{{J}_{\nu }({ha})}{{K}_{\nu }({qa})}B,\end{eqnarray} \tag{ A.14 }$$

$$\begin{eqnarray}&&B=\displaystyle \frac{{\rm{i}}}{\nu }\displaystyle \frac{{ahq}\left({n}_{1}^{2}{{qJ}}_{\nu }^{{\prime} }({ha}){K}_{\nu }({qa})+{{\rm{i}}{n}}_{2}^{2}{{hJ}}_{\nu }({ha}){K}_{\nu }^{{\prime} }({qa})\right)}{\omega ({n}_{1}^{2}-{n}_{2}^{2}){\mu }_{0}\beta {J}_{\nu }({ha}){K}_{\nu }({qa})}A.\end{eqnarray} \tag{ A.15 }$$

Furthermore we obtain an eigenvalue equation for the propagation parameter β from the condition that the linear system for A, B, C, D has to have a solution, i.e.

$$\begin{eqnarray}&&\begin{array}{l}\left(\displaystyle \frac{{J}_{\nu }^{{\prime} }({ha})}{{{haJ}}_{\nu }({ha})}+\displaystyle \frac{{K}_{\nu }^{{\prime} }({qa})}{{{qaK}}_{\nu }({qa})}\right)\left(\displaystyle \frac{{n}_{1}^{2}{J}_{\nu }^{{\prime} }({ha})}{{{haJ}}_{\nu }({ha})}+\displaystyle \frac{{n}_{2}^{2}{K}_{\nu }^{{\prime} }({qa})}{{{qaK}}_{\nu }({qa})}\right)\\ \qquad -\,{\nu }^{2}{\left(\displaystyle \frac{1}{{\left({ha}\right)}^{2}}+\displaystyle \frac{1}{{\left({qa}\right)}^{2}}\right)}^{2}\displaystyle \frac{{\beta }^{2}}{{k}^{2}}=0.\end{array}\end{eqnarray} \tag{ A.16 }$$

This leaves us with one free parameter A which we fix by the normalization condition

$$\begin{eqnarray}&&{\int }_{0}^{2\pi }{\rm{d}}\varphi {\int }_{0}^{\infty }r\,{\rm{d}}r\,{n}^{2}{\left|E\right|}^{2}=1,\end{eqnarray} \tag{ A.17 }$$

where n = n₁ for r < a and n = n₂ for r ≥ a.

A.2. Interaction with atoms

We will now derive a ME for the atoms where coherent and incoherent atom–atom interactions stem from the collective coupling to both free and guided radiation modes. We start from the Hamiltonian

$$\begin{eqnarray}&&{H}_{\mathrm{full}}={H}_{{\rm{A}}}+{H}_{{\rm{F}}}^{3{\rm{D}}}+{H}_{{\rm{F}}}^{1{\rm{D}}}+{H}_{\mathrm{int}}^{3{\rm{D}}}+{H}_{\mathrm{int}}^{1{\rm{D}}},\end{eqnarray} \tag{ A.18 }$$

with ${H}_{{\rm{A}}}={\omega }_{0}{\sum }_{i}{\sigma }_{i}^{+}{\sigma }_{i}^{-}$, where ${\sigma }_{i}^{\pm }$ denote the raising and lowering operators of the ith emitter and

$$\begin{eqnarray}&&{H}_{{\rm{F}}}^{3{\rm{D}}}=\displaystyle \sum _{\vec{k},\lambda }{\omega }_{k}{a}_{\vec{k},\lambda }^{\dagger }{a}_{\vec{k},\lambda }\quad {H}_{{\rm{F}}}^{1{\rm{D}}}=\displaystyle \sum _{\eta }\omega {b}_{\eta }^{\dagger }{b}_{\eta },\end{eqnarray} \tag{ A.19 }$$

where $\eta =\left(\omega ,\lambda ,\nu ,f\right)$ with the frequency ω, polarization λ, the mode index ν and the propagation direction f = ±1. We have assumed complete independence of the free-space and guided modes such that $\left[{a}_{\vec{k},\lambda }^{\dagger },{b}_{\eta }\right]=\left[{a}_{\vec{k},\lambda },{b}_{\eta }^{\dagger }\right]=0$. The two interaction terms in dipole and rotating wave approximation read

$$\begin{eqnarray}&&{H}_{\mathrm{int}}^{3{\rm{D}}}=-\displaystyle \sum _{i=1}^{N}{\mathop{\overrightarrow{d}}\limits^{}}_{i}\cdot {\mathop{\overrightarrow{E}}\limits^{}}^{3{\rm{D}}}({\mathop{\overrightarrow{r}}\limits^{}}_{i})\approx -{\rm{i}}\displaystyle \sum _{i=1}^{N}\displaystyle \sum _{\mathop{\overrightarrow{k}}\limits^{},\lambda }{g}_{\mathop{\overrightarrow{k}}\limits^{}\lambda }\ {\mathop{\overrightarrow{e}}\limits^{}}_{\mathop{\overrightarrow{k}}\limits^{},\lambda }\cdot \overrightarrow{\mathop{\mu }\limits^{}}\ \left({a}_{\mathop{k}\limits^{},\lambda }{{\rm{e}}}^{{\rm{i}}\mathop{\overrightarrow{k}}\limits^{}{\mathop{\overrightarrow{r}}\limits^{}}_{i}}{\sigma }_{i}^{+}-{\rm{h}}{\rm{.}}{\rm{c}}{\rm{.}}\right),\end{eqnarray} \tag{ A.20 }$$

$$\begin{eqnarray}&&{H}_{\mathrm{int}}^{1{\rm{D}}}\approx -{\rm{i}}\displaystyle \sum _{i=1}^{N}{\int }_{0}^{\infty }\displaystyle \sum _{\eta }{g}_{\eta }\ \mathop{\overrightarrow{\mu }}\limits^{}\cdot \left({\mathop{\overrightarrow{e}}\limits^{}}_{\eta }{b}_{\eta }{{\rm{e}}}^{{\rm{i}}{f}{\beta }^{\nu }{{kz}}_{i}}{\sigma }_{j}^{+}-{\rm{h}}.{\rm{c}}.\right),\end{eqnarray} \tag{ A.21 }$$

where ${g}_{\vec{k}\lambda }=\sqrt{{\omega }_{k}/2{\epsilon }_{0}V}$ and

$$\begin{eqnarray}&&{g}_{\eta }=\sqrt{\displaystyle \frac{\omega }{2{\epsilon }_{0}{v}_{g}}}\left[\vec{\mu }\cdot E({r}_{j},{\varphi }_{j})\right],\end{eqnarray} \tag{ A.22 }$$

where v_g = dω/dβ is the group velocity in accordance with [40]. As the 3D and the 1D parts are linearly superimposed and thus can be treated independently, we will now focus on the 1D part only whereas the 3D derivation can be found e.g. in [37]. Furthermore we will restrict ourselves to one mode, i.e. neglect the sum over ν for now. From this we can calculate the time evolution of our field mode by solving

$$\begin{eqnarray}&&{\dot{b}}_{\eta }={\rm{i}}\left[H,{b}_{\eta }\right]=-{\rm{i}}\omega {b}_{\eta }+\displaystyle \frac{1}{\sqrt{2\pi }}\displaystyle \sum _{j=1}^{N}{g}_{\eta }^{* }\ {{\rm{e}}}^{-{\rm{i}}\beta (\omega ){z}_{j}}\ {\sigma }_{j}^{-},\end{eqnarray} \tag{ A.23 }$$

which yields

$$\begin{eqnarray}&&{b}_{\eta }(t)={b}_{\eta }(0)\ {{\rm{e}}}^{-{\rm{i}}\omega t}+\displaystyle \frac{1}{\sqrt{2\pi }}{\int }_{0}^{t}{\rm{d}}{t}^{{\prime} }\ {{\rm{e}}}^{-{\rm{i}}\omega (t-{t}^{{\prime} })}\,{g}_{\eta }^{* }\ \displaystyle \sum _{j=1}^{N}{{\rm{e}}}^{-{\rm{i}}\beta (\omega ){z}_{j}}{\sigma }_{j}^{-}\left({t}^{{\prime} }\right).\end{eqnarray} \tag{ A.24 }$$

We continue by writing down the equation of motion for an arbitrary system operator, i.e. any operator acting on our ensemble of two-level emitters. We substitute for b_η, replace ${\sum }_{\eta }\to \int {\rm{d}}\omega {\sum }_{\lambda ,\nu ,f}$ and have

$$\begin{eqnarray}\begin{array}{rcl}\dot{O} & = & {\rm{i}}\left[{H}_{{\rm{A}}},O\right]+\displaystyle \frac{1}{\sqrt{2\pi }}{\int }_{0}^{\infty }{\rm{d}}\omega \displaystyle \sum _{\lambda ,f}\displaystyle \sum _{j}\\ & & \left({g}_{\eta }{b}_{\eta }{{\rm{e}}}^{{\rm{i}}{f}\beta (\omega ){z}_{j}}\left[{\sigma }_{j}^{+},O\right]\right.-\left.{g}_{\eta }^{* }{b}_{\eta }^{\dagger }{{\rm{e}}}^{-{\rm{i}}{f}\beta (\omega ){z}_{j}}\left[{\sigma }_{j}^{-},O\right]\right).\end{array}\end{eqnarray} \tag{ A.25 }$$

Let us now focus on the interaction part and replace ${b}_{\eta }(t)$ and ${b}_{\eta }^{\dagger }$ with the inhomogeneous part from the expression from above (the homogeneous part is the in-field and we will neglect it too)

$$\begin{eqnarray}\begin{array}{rcl}\dot{O} & = & \displaystyle \frac{1}{2\pi }{\int }_{0}^{\infty }{\rm{d}}\omega \displaystyle \sum _{\lambda ,f}\displaystyle \sum _{j,k}\\ & & \left({\int }_{0}^{t}{\rm{d}}{t}^{{\prime} }{{\rm{e}}}^{-{\rm{i}}\omega (t-{t}^{{\prime} })}{\left|{g}_{\eta }\right|}^{2}{{\rm{e}}}^{{\rm{i}}{f}\beta (\omega )({z}_{j}-{z}_{k})}{\sigma }_{k}^{-}({t}^{{\prime} })\left[{\sigma }_{j}^{+},O\right]\right.\\ & & -\,\left.{\int }_{0}^{t}{\rm{d}}{t}^{{\prime} }{{\rm{e}}}^{{\rm{i}}\omega (t-{t}^{{\prime} })}{\left|{g}_{\eta }\right|}^{2}{{\rm{e}}}^{-{\rm{i}}{f}\beta (\omega )({z}_{j}-{z}_{k})}\left[{\sigma }_{j}^{-},O\right]{\sigma }_{k}^{+}({t}^{{\prime} })\right)\end{array}\end{eqnarray} \tag{ A.26 }$$

At this point we perform the Markov approximation (similar for ${\sigma }_{j}^{-}$)

$$\begin{eqnarray}&&{\sigma }_{j}^{+}\left({t}^{{\prime} }\right)\to {{\rm{e}}}^{-{\omega }_{0}\left(t-{t}^{{\prime} }\right)}{\sigma }_{j}^{+}(t)\end{eqnarray} \tag{ A.27 }$$

and rewrite the expression

$$\begin{eqnarray}\begin{array}{rcl}\dot{O} & = & \displaystyle \frac{1}{2\pi }{\int }_{0}^{\infty }{\rm{d}}\omega \displaystyle \sum _{\lambda ,f}\displaystyle \sum _{j,k}\\ & & \left({\int }_{0}^{t}{\rm{d}}{t}^{{\prime} }{{\rm{e}}}^{-{\rm{i}}(\omega -{\omega }_{0})(t-{t}^{{\prime} })}{\left|{g}_{\eta }\right|}^{2}{{\rm{e}}}^{{\rm{i}}{f}\beta (\omega )({z}_{j}-{z}_{k})}{\sigma }_{k}^{-}\left[{\sigma }_{j}^{+},O\right]\right.\\ & & -\,\left.{\int }_{0}^{t}{\rm{d}}{t}^{{\prime} }{{\rm{e}}}^{{\rm{i}}(\omega -{\omega }_{0})(t-{t}^{{\prime} })}{\left|{g}_{\eta }\right|}^{2}{{\rm{e}}}^{-{\rm{i}}{f}\beta (\omega )({z}_{j}-{z}_{k})}\left[{\sigma }_{j}^{-},O\right]{\sigma }_{k}^{+}\right).\end{array}\end{eqnarray} \tag{ A.28 }$$

For the integral over t' we use the Sokhotski–Plemelj theorem as

$$\begin{eqnarray}&&{\int }_{0}^{t}{\rm{d}}{t}^{{\prime} }{{\rm{e}}}^{\mp {\rm{i}}(\omega -{\omega }_{0})(t-{t}^{{\prime} })}=\mp {\rm{i}}{ \mathcal P }\displaystyle \frac{1}{\omega -{\omega }_{0}}+\pi \delta \left(\omega -{\omega }_{0}\right).\end{eqnarray} \tag{ A.29 }$$

Let us pull out ${\sum }_{j,k}$ perform the sum over the propagation direction f, set $\vec{\mu }=\mu \cdot {e}_{x}$ (which kills the sum over λ) and proceed

$$\begin{eqnarray}\begin{array}{rcl}\dot{O} & = & \displaystyle \sum _{j,k}\displaystyle \frac{1}{\pi }{\int }_{0}^{\infty }{\rm{d}}\omega \\ & & \left(\left(-{\rm{i}}{ \mathcal P }\displaystyle \frac{1}{\omega -{\omega }_{0}}+\pi \delta \left(\omega -{\omega }_{0}\right)\right){\left|{g}_{\omega }\right|}^{2}\cos \left(\beta {z}_{{jk}}\right){\sigma }_{k}^{-}\left[{\sigma }_{j}^{+},O\right]\right.\\ & & -\,\left.\left({\rm{i}}{ \mathcal P }\displaystyle \frac{1}{\omega -{\omega }_{0}}+\pi \delta \left(\omega -{\omega }_{0}\right)\right){\left|{g}_{\omega }\right|}^{2}\cos \left(\beta {z}_{{jk}}\right)\left[{\sigma }_{j}^{-},O\right]{\sigma }_{k}^{+}\right).\end{array}\end{eqnarray} \tag{ A.30 }$$

In this form we can roll the equation back to the time evolution of the density operator and obtain a Lindblad-type ME reading

$$\begin{eqnarray}&&\dot{\rho }={\rm{i}}\left[\rho ,\displaystyle \sum _{j\ne k}{{\rm{\Omega }}}_{{jk}}^{1{\rm{D}}}{\sigma }_{j}^{+}{\sigma }_{k}^{-}\right]+\displaystyle \sum _{{jk}}\displaystyle \frac{{{\rm{\Gamma }}}_{{jk}}^{1{\rm{D}}}}{2}\left(2{\sigma }_{j}^{-}\rho {\sigma }_{k}^{+}-{\sigma }_{j}^{+}{\sigma }_{k}^{-}\rho -\rho {\sigma }_{j}^{+}{\sigma }_{k}^{-}\right),\end{eqnarray} \tag{ A.31 }$$

where we have neglected the infinite Lamb shift for j = k. The two coupling parameters are calculated as

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Gamma }}}_{{jk}}^{1{\rm{D}}}}{2}=\displaystyle \frac{1}{\pi }{\int }_{0}^{\infty }{\rm{d}}\omega \pi \delta \left(\omega -{\omega }_{0}\right){\left|{g}_{\omega }\right|}^{2}\cos \left(\beta {z}_{{jk}}\right)=\displaystyle \frac{{\omega }_{0}{\mu }^{2}{\left[{E}_{x}\right|}^{2}}{2{\epsilon }_{0}{v}_{g}}\cos \left(\beta {z}_{{jk}}\right)\end{eqnarray} \tag{ A.32 }$$

and for the dipole–dipole shift we employ the residue theorem, create a pole at ${\omega }_{0}+{\rm{i}}\epsilon$ with $\epsilon \to 0$, close the integration path in the upper half plane and take only the real part of the integral, such that

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{jk}}^{1{\rm{D}}}=-\mathop{\mathrm{lim}}\limits_{\epsilon \to 0}{\mathfrak{R}}\displaystyle \frac{1}{\pi }\oint {\rm{d}}\omega \displaystyle \frac{{\left|{g}_{\omega }\right|}^{2}}{\omega -\left({\omega }_{0}+{\rm{i}}\epsilon \right)}\displaystyle \frac{1}{2}\left({{\rm{e}}}^{{\rm{i}}\beta {z}_{{jk}}}+{{\rm{e}}}^{-{\rm{i}}\beta {z}_{{jk}}}\right)\end{eqnarray} \tag{ A.33 }$$

$$\begin{eqnarray}\begin{array}{rcl} & = & -\displaystyle \frac{1}{2\pi }{\mathfrak{R}}\,2\pi {\rm{i}}{\left|{g}_{{\omega }_{0}}\right|}^{2}{{\rm{e}}}^{{\rm{i}}\beta {z}_{{jk}}}\\ & = & \displaystyle \frac{{\omega }_{0}{\mu }^{2}{\left[{E}_{x}\right|}^{2}}{2{\epsilon }_{0}{v}_{g}}\sin \left(\beta {z}_{{jk}}\right).\end{array}\end{eqnarray} \tag{ A.34 }$$

At this point it is important to realize, that all functional dependencies, i.e. E_x, v_g and β are evaluated at the transition frequency ω₀ and thus become a mere number.

The couplings as a function of the interatomic distance ξ = k₀d are depicted in figure A1. For the 1D coupling we show the spatial behavior of one mode. Taking several modes into account leads to a superposition of trigonometric functions with different periodicities and amplitudes as is discussed in the paper.

Zoom In Zoom Out Reset image size

A.3. Multiple modes

In order to allow for multiple modes in our model, we modify the interaction Hamiltonian to accommodate for several modes, i.e.

$$\begin{eqnarray}&&{H}_{\mathrm{int}}=-\displaystyle \frac{{\rm{i}}}{\sqrt{2\pi }}{\int }_{0}^{\infty }{\rm{d}}\omega \displaystyle \sum _{f,\lambda ,j}\displaystyle \sum _{\nu }\left({g}_{\eta }{b}_{\eta }{{\rm{e}}}^{{\rm{i}}{f}{\beta }^{\nu }(\omega ){z}_{j}}{\sigma }_{j}^{+}-{\rm{h}}.{\rm{c}}.\right).\end{eqnarray} \tag{ A.35 }$$

Doing this with the photonic operators for different modes commuting we get slightly modified coupling constants in a very similar derivation, which read

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{jk}}=\displaystyle \sum _{\nu }\displaystyle \frac{{\omega }_{0}{\mu }^{2}{\left[{E}_{x}^{\nu }\right|}^{2}}{{\epsilon }_{0}{v}_{g}^{\nu }}\cos \left({\beta }^{\nu }{z}_{{jk}}\right),\end{eqnarray} \tag{ A.36 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{jk}}=\displaystyle \sum _{\nu }\displaystyle \frac{{\omega }_{0}{\mu }^{2}{\left[{E}_{x}^{\nu }\right|}^{2}}{2{\epsilon }_{0}{v}_{g}^{\nu }}\sin \left({\beta }^{\nu }{z}_{{jk}}\right).\end{eqnarray} \tag{ A.37 }$$

For an impression of how the different modes will propagate, figure A2 depicts the dispersion relation as a function of the ratio between the fiber's radius and the wavelength of the propagating light. In a thin nano fiber, where $a/{\lambda }_{0}\lt 0.5$, only a few modes will be allowed, up to the extremal case, where only the fundamental mode is permitted. The wider the fiber becomes in relation to the light's wavelength, the more modes will be able to propagate. For $a/{\lambda }_{0}\approx 3$ we observe about 15 different modes. This implies that the tailoring of the interaction cannot only be done by placing the atoms at suitable distances to each other, but also by choosing a desired fiber radius.

Zoom In Zoom Out Reset image size

A.4. Numerical parameters

We now define a one-dimensional decay rate independent of the particular mode or the particular position of the atom as a prefactor

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{1{\rm{D}}}=\displaystyle \frac{{\omega }_{0}{\mu }^{2}}{{\epsilon }_{0}{cA}},\end{eqnarray} \tag{ A.38 }$$

where A is the cross section of the fiber A = a²π and is compensated for by the normalization of the electric field. We can then write

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{jk}}^{1{\rm{D}}}={{\rm{\Gamma }}}_{1{\rm{D}}}\displaystyle \sum _{\nu }{\chi }^{\nu }\,\cos \left({\bar{\beta }}^{\nu }{k}_{0}{r}_{{jk}}\right)\end{eqnarray} \tag{ A.39 }$$

and with $\beta =k\bar{\beta }$, where $\bar{\beta }$ is the dimensionless scaled propagation parameter

$$\begin{eqnarray}&&{\chi }^{\nu }=\left({\bar{\beta }}^{\nu }+{k}_{0}{\bar{\beta }}^{{\prime} }\right)\displaystyle \frac{A}{{C}^{\nu }}{\left[{e}_{x}^{\nu }\right|}^{2}.\end{eqnarray} \tag{ A.40 }$$

The parameters are evaluated numerically for each mode ν.