Modified dipole-dipole interactions in the presence of a nanophotonic waveguide

When an emitter ensemble interacts with the electromagnetic field, dipole-dipole interactions are induced between the emitters. The magnitude and shape of these interactions are fully determined by the specific form of the electromagnetic field modes. If the emitters are placed in the vicinity of a nanophotonic waveguide, such as a cylindrical nanofiber, the complex functional form of these modes makes the analytical evaluation of the dipole-dipole interaction cumbersome and numerically costly. In this work, we provide a full detailed description of how to successfully calculate these interactions, outlining a method that can be easily extended to other environments and boundary conditions. Such exact evaluation is of importance as, due to the collective character of the interactions and dissipation in this kind of systems, any small modification of the interactions may lead to dramatic changes in experimental observables, particularly as the number of emitters increases. We illustrate this by calculating the transmission signal of the light guided by a cylindrical nanofiber in the presence of a nearby chain of emitters.


Introduction
An ensemble of emitters coupled to a common environment displays collective behaviour.This includes the enhanced and inhibited emission of photons from the ensemble (so-called super and subradiance, respectively) [1,2,3], and the emergence of induced dipole-dipole interactions among the emitters [4,5,6].The rates associated with these processes can be modified via the specific choice of the environment and its boundary conditions, which may be tailored in order to obtain desirable properties.For example, introducing a dielectric or metallic structure in close proximity to an atomic gas, modifies the local electromagnetic spectrum and hence the collective behaviour of the ensemble [7,8,9,10,11].This effect was noticed by Purcell in 1950 [12] and has later been verified in various experiments involving atoms or electrons in cavities, and molecules coupled to dielectric structures [13,14,15,16].Among these structures, so-called nanophotonic waveguides, such as single mode optical nanofibers [17,18] or integrated photonic nanostructures [19,20,21,22], particularly stand out since they provide strong and homogenous coupling between the emitters and long coherence times [23], and also due to the confined or guided field modes carried by these nanostructures, which can lead to propagation-directiondependent (chiral) emission [24].The translationally invariant nature of these guided modes gives rise to infinitely ranged coherent dipole-dipole interactions and incoherent couplings between the emitters.These all-to-all interactions have recently facilitated, for example, the observation of super-and subradiance [25,26,27], the study of non-trivial topology [28,29,30], the realization of long-lived photon storage and multiple photon bound states [31,32,33,34,35] and the investigation of new dynamical phases and phase transitions [36,37].The modes propagating in the space outside of the nanostructure are often referred to as radiation or unguided modes.Unfortunately, since these modes can propagate in any direction in space, their contribution to the coherent and incoherent interactions between the emitters is highly non-trivial, containing, for example, frequency integrals that either converge very slowly or not at all.Hence, they are often approximated by the free-space modes [26,38], an approximation that, however, inevitably breaks down when the emitters are close to the nanostructure.A number of works have nicely provided expressions for these contributions in terms of complex integration contours or calculated them numerically under strong approximations, like low-frequency cut-offs [31,39,40,41], but up to now how to compute these contributions exactly remains an open challenge.
In this paper, we address precisely this challenge by not only providing a comprehensive guide for the analytical calculation of guided and unguided mode contributions to the dipoledipole interactions and collective dissipation induced among emitters trapped near a single mode optical nanofiber, but also supplying -to the best of our knowledge for the first time-a method for the numerical calculation of the exact contribution of the radiation modes that overcomes convergence issues that typically make their accurate numerical determination extremely challenging.The method that we outline here can be generally used for the calculation of dipoledipole interactions induced by other nanostructures, given that the mode profile functions are known [11,42,10,32,43].We show that, due to the presence of the nanofiber, the radiation field in its vicinity can be significantly altered, giving rise to contributions to the dipole-dipole interactions that significantly differ from the free-field counterparts.Owing to the collective character of both coherent and incoherent interactions, when measuring experimentally observable properties such as the transmission signal of fiber-guided light, these differences are even more evident as the number of emitters is increased.

System and master equation
We consider an ensemble of N identical emitters placed near the surface of a nanophotonic waveguide (see Fig. 1a).Each emitter is modelled as a two-level system with ground and excited states |g⟩ and |e⟩, respectively, separated by an energy ℏω a .Within the dipole and Born-Markov approximations, the dynamics of the reduced density matrix ρ, which contains the internal degrees of freedom of the emitters, is determined by the Each emitter is a two-level system with an excited (|e⟩) and a ground (|g⟩) state separated by an energy ℏω a .b): Due to their coupling to the electromagnetic field modes, the emitters interact through the exchange of virtual photons.c): These interactions have contributions stemming from the guided (V gd ) and the radiation (V rd ) modes.
master equation [10,38] where σ α = |g α ⟩⟨e α | is the spin-1/2 ladder operator for the α-th emitter.The first term in the master equation describes dipole-dipole or exchange interactions, which are mediated by the exchange of virtual photons between the emitters (see Fig. 1b).The exchange rate between the α-th and β-th emitter is given by the magnitude of the coefficient V αβ .The second term of (1) describes incoherent photon emission or dissipation in the system, which in general possesses a collective character.Here, the diagonal terms of the dissipation matrix Γ αβ describe the single-emitter spontaneous decay rates, which for the α-th emitter is thus given by Γ αα .The sign and magnitude of the off-diagonal elements Γ αβ for α ̸ = β are responsible for the collective character of the emission, its landmark being the well-known superand subradiant photon emission [44,45].
Our aim is to obtain the coefficient matrices V αβ and Γ αβ , which in turn fully determine all properties in the system.These can be determined via the electromagnetic Green's tensor Ḡ(r α , r β , ω) of the environment evaluated at the emitter positions r α and r β as [10,11,46] where P denotes the principal value and d α is the transition dipole moment for the α-th emitter.The electromagnetic Green's tensor is the solution to the wave equation where k = ω/c and ϵ(r, ω) is the dielectric constant of the medium.The Green's tensor satisfies the Schwarz reflection principle Ḡ * (r, r ′ , ω) = Ḡ(r, r ′ , −ω * ) and the Onsager reciprocity, ḠT (r, r ′ , ω) = Ḡ(r ′ , r, ω).As discussed in Appendix A, the complex function ω 2 Ḡ(r, r ′ , ω) is by definition holomorphic in the complex upper half-plane and bounded at infinity.Hence, the dipole-dipole interaction may be written in terms of the real part of the Green's tensor where we have made use of the Kramers-Kronig relations.
As shown above, for the calculation of the dipole-dipole and dissipation matrices we need to be able to separate the real and imaginary components of the Green's tensor, respectively.As the Green's tensors we consider are in general written as integrals over the frequency, to do so we introduce a small imaginary shift to the frequency Ḡ(r, r ′ , ω) = lim ϵ→0 + Ḡ(r, r ′ , ω + iϵ) [47] and use the Sokhotski-Plemelj theorem applied on the real line to find the imaginary part and compute the real part using the Kramers-Kronig relation.However, as discussed in Appendix A, this is valid only when the function we consider is well-defined and holomorphic in the frequency complex upper half-plane (including the real line).Decomposing the Green's tensor into its longitudinal and transverse components one can see that, while the transverse component is indeed well-behaved, the longitudinal component is singular at ω → 0 [48].To illustrate this, let us consider the case where the medium is the vacuum (i.e.free space with no boundaries).Here, ϵ(r, ω) = 1 everywhere and the solution to the wave equation (2) is known analytically and given by 2  , where r = r − r ′ , r = |r|, and r = r/r.Decomposing the vacuum Green's tensor into its longitudinal and transverse components as Ḡ0 (ω) = Ḡ0∥ (ω) + Ḡ0⊥ (ω) [49,50] (see Appendix A), we observe that, while the transverse component is a holomorphic function in the complex upper halfplane, the longitudinal component is not, as it contains a second order pole at ω = 0.This in turn means that the usual form of the Kramers-Kronig relation cannot be used to find the real part of the full vacuum Green's tensor.It can rather only be validly applied to its tranverse part.With this, we find that where we have used the fact that the longitudinal part is purely real.We will consider this relation when treating the more complicated situation of a cylindrical nanofiber.

Green's tensor for a cylindrical nanofiber
The Green's tensor can be expressed in terms of the eigenmodes of the positive frequency part of the electric field as Here, E n (r, ω) are the eigenvectors of the her- with eigenvalues λ n .This representation allows us to easily split the contribution to the Green's tensor coming from the eigenmodes guided by the waveguide and the unguided (or radiative) ones as Ḡ(r, r ′ , ω) = Ḡgd (r, r ′ , ω) + Ḡrd (r, r ′ , ω) [51,52,53], which leads in turn to the decomposition of the dipole-dipole and collective dissipation matrices: V αβ = V gd αβ +V rd αβ and Γ αβ = Γ gd αβ +Γ rd αβ , see Fig. 1c.Let us analyze these contributions separately for the case of a cylindrical dielectric nanofiber of radius r f characterized by a refractive index n 1 placed in an infinite vacuum of refractive index n 2 = 1.To do so, we will use a cylindrical coordinate system (r, ϕ, z), centered in the fiber core with the z-direction being along the fiber (see Fig. 1a).

Guided modes
Using the expression for the guided modes of the electric field given in Appendix B, the guided contribution to the Green's tensor is given by where (µ) ≡ (β, l, f ) are the labels of the guided modes.Here, l = ±1 is the polarization of the mode, f = ±1 is the propagation direction for the guided modes in the fiber and e (µ) (r) are the guided profile functions (note, that we assume that the nanofiber supports only the fundamental HE 11 guided mode).In addition, β is the longitudinal propagation constant in the fiber, which for each value of ω is determined via the fiber eigenvalue equation (see Appendix B), q is a variable associated to the field outside the nanofiber and k ′ = ω ′ /c.The expression for the Green's tensor may be written in terms of the mode frequency ω by using the relation where β ′ (ω) = dβ dω .The guided contribution to the Green's tensor is purely transverse.We can therefore apply the Sokhotski-Plemelj theorem to find its imaginary part which gives a guided contribution to the collective decay matrix that reads The real part of the guided Green's tensor is found by using the usual Kramers-Kronig relation The principal value integral in this equation can be calculated analytically choosing an appropriate contour [26,38], and one finds that is the guided contribution to the dipole-dipole interaction.

Radiation modes
The radiation modes of the electric field can be both transverse and longitudinal, and the radiation contribution to the Green's tensor may be decomposed as Ḡrd (ω) = Ḡrd⊥ (ω) + Ḡrd∥ (ω).On the other hand, one can also separate the Green's tensor as a sum of the vacuum modes and the ones scattered by the nanofiber as Ḡrd (ω) = Ḡ0 (ω) + Ḡsc (ω).The scattered modes are, however, purely transverse, allowing to replace the longitudinal radiation Green's tensor with the well known longitudinal part of the vacuum Green's tensor and write Using the expression for the transverse radiation modes of the electric field given in Appendix C, the radiative contribution to the Green's tensor is given by where (ν) ≡ (β, q, m, l) are the labels of the radiation modes, with l = ±1 and m = 0, ±1, ±2, ... labelling the polarization and order of the mode.
In contrast to the guided mode case, β = k cos θ is here a continuous variable for each value of the frequency ω and the variable q = k sin θ is characteristic of the field outside the nanofiber.Finally, e (ν) (r) is the radiation profile function of the transverse electric field.The Green's tensor may be rewritten in terms of the mode frequency ω and the angle θ as Note, that, as explained in Appendix C, we have chosen specifically the eigenmode decomposition adopted in [53], which differs from, for example, the one adopted in [38] in the accompanying normalization factors.However, let us point out that the two formulations are equivalent and hence give the same results.
Similarly as for the guided modes, we use the Sokhotski-Plemelj theorem to find the imaginary part of the radiation Green's tensor using the fact that the transverse component is holomorphic in the upper half complex plane and well-behaved on the real frequency axis.The imaginary part becomes then such that the corresponding contribution to the dissipation matrix reads where now β a = ωa c cos θ and q a = ωa c sin θ are both functions of the angle θ.
We use the imaginary part of the Green's tensor evaluated at ω = ω a to illustrate in Fig. 2 how much the presence of the nanofiber disturbs the field around a point dipole compared to the situation where the dipole is in free space.In particular, we can clearly observe that the field is most dramatically modified when the dipole moment is perpendicular to the axis of the nanofiber (y and x directions, second and third row in Fig. 2).
The real part of the radiation Green's tensor can be found by applying the (modified) Kramers-Kronig relation to the transverse part of the radiation Green's tensor, ii } at a distance of 100 nm from the surface of a silica nanofiber with r f = 250 nm.The frequency ω a = 2πc/λ a corresponds to the D2 transition in Cs, with λ a = 852 nm.The refractive index of the silica nanofiber is modeled using a Sellmeier equation [54].
such that the radiation contribution to the dipoledipole interaction is given by Due to the complicated frequency dependence of the imaginary part of the radiation Green's tensor, an analytical treatment of the principal value integral in equation ( 4) is nontrivial and requires the introduction of branch cuts and the application of contour integration techniques [31].Due to these complications, in the literature typically this contribution is approximated by its free-field counterpart, which is known analytically [26,38].However, as illustrated by Fig. 2, the field might actually be strongly modified by the presence of the nanofiber, and these modifications, even when small, may lead to large differences in observables when treating systems with many emitters.For these reasons, in the following we describe how to efficiently calculate numerically the exact value of V rd αβ .

Numerical evaluation of the radiative dipole-dipole interactions
The first challenge one encounters to evaluate numerically the integral in ( 4) is the singularity in the integrand, happening at ω = ω a .This singularity can be avoided by using a fast Fourier transform, as introduced in Refs.[55,56], and using the Schwarz reflection principle, such that V rd αβ may be written as However, the main obstacle for the calculation of V rd αβ lies in the imaginary part of the Green's tensor, in particular in its behavior at large values of ω.To illustrate this, let us consider again the Green's tensor in vacuum (3).In the limit of large ω, i.e. when kr ≫ 1, the imaginary part of this tensor yields Im{ Ḡ0 (r, r ′ , ω)} ≈ 1 4πr (1 − 3 r2 ) cos kr kr Here we can see that the integrand in the frequency integral in (5) for large ω ≫ ω a is an oscillating function with period proportional to ω a λ a /r αβ , where λ a = 2πc ωa is the transition wavelength.The amplitude of the oscillations either decays very slowly, as 1/ω (e.g. in vacuum when the dipoles are aligned with the displacement vector, i.e. | dα • rαβ | = 1, see Fig. 3a), or is constant.In the former case, reaching convergence of the integral over frequency requires then an extremely high cut-off frequency ω c ≫ ω a λ a /r αβ , which grows the closer the emitters are to each other.Moreover, in general, even a large value Single cut-off Averaging Figure 3: Numerical integration.Imaginary part of the Green's tensor in vacuum (in arbitrary units) as a function of ω/ω a for two dipoles at a distance a from one another pointing a) parallel and b) perpendicular to the displacement vector between them for a/λ a = 1 and a/λ a = 1/2 (black and blue solid lines, respectively).The insets show the relative error in % when calculating the vacuum dipole-dipole interaction V 0 αβ as a function of a/λ a by using a direct numerical integration with a single cut-off frequency ω c /ω a = 10λ a /a (black dashed) and using the averaging method with a range of cut-offs ω c /ω a ∈ [2λ a /a, 10λ a /a].The sharp peaks occur at the zeros of V 0 αβ .
for ω c does not necessarily ensure convergence (see e.g.Fig. 3b for emitters in the free-field with | dα • rαβ | = 0, where the amplitude of the oscillations remains constant for large ω).
In addition to this issue, note that for the numerical calculation of the imaginary part of the Green's tensor in the presence of the cylindrical nanofiber, in principle an infinite sum of modes m should be considered.In practice, the sum is truncated at some finite value m c .However, for a fixed distance of the emitters with respect to the nanofiber surface, the amount of modes necessary to ensure convergence grows dramatically with increasing ω.This means that the computational cost of increasing the cut-off frequency is enormous.In order to circumvent this problem, we make use of the periodic nature of the integrand in the limit of large ω.This allows us to approximate the result of the frequency integral as the average of the outcomes for a range of cut-off frequencies comprising a few oscillations.This results in a much faster convergence of the integral when the integrand decays as 1/ω, as the maximum cut-off frequency necessary for convergence is relatively low.Most importantly, when the amplitude of the oscillations remains constant, this allows to obtain convergence which would remain otherwise elusive.
In order to illustrate the power of this approach, we have applied it to calculate numerically the dipole-dipole interactions in vacuum.We show in the insets of Fig. 3 the relative error of the numerical calculation compared to the exact one (known analytically in this case), using a "single cut-off" direct integration of (4) and an "averaging" method for the frequency integration in expression (5).It is evident that the second approach gives dramatically better results, particularly when the dipoles are perpendicular to their displacement vector.This encourages us to apply this numerical method to the nanofiber case, where the main features of the Green's tensor as a function of the frequency are similar to the vacuum ones.
We have made publicly available a code that follows this recipe to calculate both guided and radiation contributions to both dipole-dipole and incoherent interactions in [57].

Modified dipole-dipole interaction between two atoms near a nanofiber
We apply here the numerical method presented above to the simple case of two identical twolevel atoms placed at a fixed distance from the fiber surface x a , fixed azimuthal angle ϕ = 0 and fixed separation a. Figure 4 shows the radiation contribution to the dipole-dipole interaction V rd for distances x a = 50nm and x a = 100nm from the fiber surface compared to the vacuum counterpart V 0 as a function of a/λ a for three different dipole orientations.The dipole-dipole interaction in all cases is normalized by the singleatom spontaneous decay rate in vacuum, γ = |d| 2 ω 3 a /(3c 3 πϵ 0 ℏ).As expected from the results for the imaginary The radiation contribution to the dipole-dipole interaction V rd (in units of the single particle decay rate γ) for a distance of x a = 50 nm (red dotted) and x a = 100 nm (blue dashed) from the fiber surface as a function a/λ a .The black solid line is the analytic solution in the case of vacuum.Each panel displays the results for a different orientation of the dipoles with respect to the surface of the nanofiber: parallel, binormal and normal for the upper, middle and lower panel, respectively.The parameters are identical to those used in Fig. 2.
part of the Green's tensor (Fig. 2), the most important modifications of the dipole-dipole interactions are found when the atomic dipole moment orientation is perpendicular to both the axis and the surface of the nanofiber (bottom panel).Here, the differences between the calculated radiation components and their vacuum counterparts are not only large (e.g., 70% difference at a/λ a = 1), but also persist at large distances between the atoms i.e. at a ≥ λ a .
Finally, to further illustrate the power of our approach, let us consider a pair of dipoles separated by a = λ parallel to the displacement vector (Figs.3a and 4a).Using our numerical method, here we calculate the vacuum interaction to a pre-cision of 0.004% error using an upper cut-off frequency of ω c = 10ω a .To get the same precision using a standard numerical integration method, a cut-off frequency of ω c = 40ω a would be needed.
For ω = 10ω a we have to include around m c = 90 modes, while for ω = 40ω a around 260 modes are needed (needing about six times longer to calculate the imaginary part of the Green's tensor at that frequency point alone).

Effect of the collective interaction on the transmission of fiber-guided light
The effect of the modifications that we find in the dipole-dipole interaction between a pair atoms becomes greatly amplified when considering observables in a large ensemble of many emitters.
The key to understand this amplification is that the emitters couple collectively to its environment.This in turn means that each of the N 2 terms V αβ and Γ αβ of the dipole-dipole and dissipation matrices play a role in the (collective) coherent and incoherent dynamics of the system.Hence, small variations of the coefficients V αβ and Γ αβ in the master equation ( 1) have a large effect on observables such as the photon emission rate or the spectrum of the light that is absorbed and emitted from the ensemble.
We illustrate this now by investigating the transmission signal of fiber-guided light when a periodic chain of N atoms is placed in the vicinity of the fiber surface (see Fig. 5).Here, in the cylindrical coordinate system, the coordinates of atom α are given by r α = (x a , 0, (α − 1)a).The system is driven by a weak probing field of frequency ω p detuned from the atomic resonance frequency by the detuning ∆ = ω p − ω a and with Rabi frequency Ω.The light field enters the nanofiber from its left, is guided through the nanofiber, interacts with the atomic chain, and its transmitted signal is measured at a position z to the right of the chain in the fiber as [10,31] Here Êgd R is the right propagating component of the guided probe field, given by a sum of the input and scattering components as where γ gd = Γ gd αα is the single-atom decay rate of each atom into the guided modes.
In the single excitation sector, the atomic wave function may be written as In the low saturation regime (weak probe laser), c G (t) ≈ 1, and inserting the ansatz into the master equation ( 1), the time evolution of the probability amplitude coefficients c α e (t) can be found to be determined by where the components of η are η α = Ωe iβzα and Γ = Γ αα is the total single-atom decay rate (sum of the decay rates into guided and radiation modes).Finally, the effective Hamiltonian H eff given by has been introduced for α ̸ = β.
In Figure 5, we show the guided light transmission signal (6) as a function of ∆ in the stationary state [where ċe (t) = 0] when a periodic chain of N = 5, 10 and 20 atoms with nearest neighbor separation a/λ a = 0.1 is placed at x a = 50, 100 and 200 nm from the surface of the nanofiber.We compare in all cases the transmission with the result obtained under the approximation that the radiation component of the dipole-dipole interaction is given by the free-field V 0 αβ .One can easily observe here that indeed large deviations arise in the transmission spectrum.These differences become more evident the closer the atoms are to the nanofiber (where each element V rd αβ differs more from V 0 αβ , see Fig. 4) and the larger the number N of atoms in the chain (where more elements V rd αβ participate in the collective dynamics).In all cases one can observe, not only measurable shifts of the (subradiant) narrow and (superradiant) broad resonance peaks, but an overall deformation of the spectrum as N increases.Note, that in current experiments typically several hundreds of atoms can be trapped near the nanofiber [58], making the spectra differ even more substantially.

Conclusions and outlook
We have provided a detailed recipe for the analytical and numerical calculation of the interaction of an ensemble of emitters in the presence of a cylindrical nanofiber.In particular, we have shown that the dipole-dipole interactions mediated by the radiation modes outside the nanofiber differ significantly (up to a 70% in some of the cases considered) from the ones obtained in vacuum, specially when the emitters' transition dipole moments are oriented normal or binormal to the fiber surface.These differences affect substantially the collective properties of the system, such as the transmission signal of fiber-guided light, where large shifts in the resonances are observed as the number of emitters is increased.Note, however, that we have calculated the transmission in the single-excitation limit (linear optics regime).Beyond this limit, we expect that, while features related to superradiance will not dramatically affected by modifications in the dipole-dipole interactions [59], these will become crucial when studying the excitation and description of many-body long-lived (subradiant) states, e.g. for the realisation of quantum memories [60,61].
A The Green's tensor and the applicability of the Kramers-Kronig relation In this Appendix we discuss the applicability of the Kramers-Kronig relation for the Green's tensor.
The application of the Kramers-Kronig relation to a complex function χ(ω) requires that the function satisfies two conditions.The function has to be holomorphic in the upper half plane, i.e. the function contains no poles in the complex upper half plane and the function has to converge to a finite value as |ω| → ∞ in the upper half plane.
The Green's tensor converges as |ω| → ∞ in the upper half plane, but it can have poles on the real line at ω = 0. Hence, the Green's tensor does not in general satisfy the conditions for the Kramers-Kronig relation.However by multiplying the Green's tensor by ω 2 , the small and large frequency limits may be applied [62] where the components of M are M ij < ∞, which is exactly the Kramers-Kronig conditions.Therefore, the Kramers-Kronig relation may always be applied to the function ω 2 Ḡ(r, r ′ , ω), but not always to the Green's tensor itself.
The transverse Green's tensor is in general holomorphic in the complex upper half-plane including the real line and hence satisfies the conditions for the Kramers-Kronig relation, while the longitudinal part is in ω = 0.The transverse and longitudinal part of the Green's tensor may be found by means of the transverse and longitudinal projection operators [48,63] where the projection operators are defined as As an example, let us now decompose the vacuum Green's tensor into its longitudinal and transverse components as [49,50] Ḡ0∥ (r, ω) =(3 r2 − 1) where r = r − r ′ , r = |r| and r = r/r.Here, we can see that, while the transverse component of the Green's tensor is a holomorphic function in the complex upper half-plane, the longitudinal component is not, as it contains a second order pole at ω = 0.

B Guided modes
The guided modes of the nanofiber are denoted by (µ) = (β, l, f ), where β is the longitudinal propagation constant, l = ±1 is the polarization and f = ±1 is the propagation direction in the fiber.The guided modes of the electric field at a position r = (r, ϕ, z) may be expressed as where e (µ) (r) is the guided profile function of the fiber with components in cylindrical coordinates, for r < r f , given by e (µ) r (r) =iC Here, κ 2 = k 2 n 2 1 − β 2 and q 2 = β 2 − k 2 characterizes the field inside and outside the fiber respectively, k = ω/c is the free space propagation constant and n 1 is the refractive index of the fiber.The longitudinal propagation constant β of the fiber may be determined by solving the fiber eigenvalue equation The functions J m (x) and K m (x) are the Bessel function of first kind and modified Bessel function of second kind respectively, which has the following form in the asymptotic limit (|x| ≫ 1) The primed Bessel functions are their derivatives The parameter s is defined as where k 0 = ω 0 /c and ω 0 is the argument frequency of the guided Green's tensor.
The profile function of a guided mode is and is normalized

C Radiation modes
The radiation modes of the nanofiber are denoted by (ν) = (β, q, m, l), where β is the longitudinal propagation constant, m = 0, ±1, ±2, ... is the mode order and l = ±1 is the polarization.The characteristics of the field inside and outside the fiber are κ 2 = k 2 n 2 1 − β 2 and q 2 = k 2 − β 2 respectively.Unlike for the guided modes, the variables β and q for the radiation modes are continuous variables for each value of k and may be expressed in terms of an angle θ β = k cos θ, q = k sin θ.
The radiation modes of the electric field at position r = (r, ϕ, z) may be expressed as E (ν) (r) = e (ν) (r)e imϕ e iβz .
We will only consider the transverse modes of the radiation electric field, ∇ • [ϵ(r)E (ν) (r)] = 0.In this case, e (ν) (r) is the transverse radiation profile function of the fiber with components in cylindrical coordinates, for r < r f , given by e C j H (j) m (qr).
By choosing the constant B = ilηA, the normalization of the profile function leads to the equations and

Figure 1 :
Figure1: Dipole-dipole interactions between emitters next to a nanofiber.a): N emitters in vacuum are placed in the vicinity of a nanofiber with radius r f .Each emitter is a two-level system with an excited (|e⟩) and a ground (|g⟩) state separated by an energy ℏω a .b): Due to their coupling to the electromagnetic field modes, the emitters interact through the exchange of virtual photons.c): These interactions have contributions stemming from the guided (V gd ) and the radiation (V rd ) modes.

Figure 2 :
Figure 2: Modified radiation field near a nanofiber.Diagonal component of the imaginary part of the Green's tensor Im{G ii (r, r ′ , ω a )} with i = x, y, z (bottom, middle and upper row, respectively, all in arbitrary units) around a dipole placed at r ′ = 0 as a function of the coordinate r in the zx-plane (left and middle column) and in the yx-plane (right column).The cross in the middle of the dipole should be interpreted as the dipole pointing into the plane.The left column shows Im{G 0ii } in vacuum (free field), while the middle and right columns display Im{G rd ii } at a distance of 100 nm from the surface of a silica nanofiber with r f = 250 nm.The frequency ω a = 2πc/λ a corresponds to the D2 transition in Cs, with λ a = 852 nm.The refractive index of the silica nanofiber is modeled using a Sellmeier equation[54].

Figure 4 :
Figure4: Modified dipole-dipole interactions.The radiation contribution to the dipole-dipole interaction V rd (in units of the single particle decay rate γ) for a distance of x a = 50 nm (red dotted) and x a = 100 nm (blue dashed) from the fiber surface as a function a/λ a .The black solid line is the analytic solution in the case of vacuum.Each panel displays the results for a different orientation of the dipoles with respect to the surface of the nanofiber: parallel, binormal and normal for the upper, middle and lower panel, respectively.The parameters are identical to those used in Fig.2.

Figure 5 :
Figure5: transmission signal.Transmission signal T in the stationary state of fiber-guided light driving the atoms with Rabi frequency Ω and detuning ∆ for a chain of N atoms with lattice constant a = 0.1λ a placed at a distance x a from the surface of the nanofiber.The black solid line represents the transmission T 0 using the vacuum instead of the full radiation dipole-dipole interactions.

m
(x) are the Bessel function of first kind and Hankel functions of j-th kind respectively, which takes the following form in the asymptotic limit(|x| ≫ 1) m (x) ≈ 2 πx e i(x− mπ 2 − π 4 ) , H (2) m (x) ≈ 2 πx e −i(x− mπ 2 − π 4 ) .The prime denotes the derivative of the functions, J ′ m (x) = d dx J m (x) and H The coefficients C j and D j are related to the coefficients A and B C j = (−1) j iπq 2 r f 4 (AL j + iµ 0 cBV j ), D j = (−1) j−1 iπq 2 r f 4 (iϵ 0 cAV j − BM j ),