Photon-antibunching in the fluorescence of statistical ensembles of emitters at an optical nanofiber-tip

Single-photon sources are essential building blocks of devices in quantum information science [1]. An ideal source emits exactly one photon into a single optical mode at an arbitrarily chosen time, with each photon being indistinguishable from all others generated by the source. The demand that no more than a single photon is emitted at a time typically requires the use of single emitters, for instance single atoms, molecules, quantum dots, or single defects in crystal lattices [2–10]. The fabrication of such single emitter sources is however technically demanding [11], which motivates the search for new approaches that are simpler to build. The idea presented in this paper is to collect the fluorescence from a nanoscale volume that contains a statistically varying number of N emitters with an average on the order of $\left\langle N\right\rangle \sim 1$. This contains the possibility that more than one emitter is present within the volume. However, in this case the excitation and emission process can be dramatically altered due to interactions between the emitters, [12]. By using cold atoms as emitters the Doppler effect of moving atoms can be eliminated. We theoretically analyze the light statistics of the fluorescence of a statistical mixture of one and two atoms with and without interactions by simulating coherent excitation by a laser field emerging from the fiber tip, and incoherent emission using a quantum Monte Carlo approach. From the resulting time stamps at which photons are emitted the corresponding g⁽²⁾-function is calculated. The paper is organized as follows: section 2 defines the problem of how statistical ensembles of emitters can emit light with subpoissonian statistics and large brightness, which at first glance seem to be contradicting requirements. Section 3 contains the theoretical model used for calculating the fluorescence from two atoms interacting with a classical pump light field and by dipole–dipole interactions among them. In section 4 an analytic solution of the Hamilton-operator (neglecting the decay terms) shows that an increasing interaction strength leads to a retardation of the simultaneous excitation of both atoms which is eventually responsible for the evolution of a subpoissonian photon statistics. Section 5 contains results of the quantum Monte Carlo simulations with interaction strength δ₁₂ and collective correction γ₁₂ to the decay rate as free parameters. The corresponding g⁽²⁾ photon-autocorrelation functions exhibit photon-antibunching for increasing values of δ₁₂. Eventually, section 6 deals with the experimental situation of Rubidium atoms at the apex of a nanoscale fiber tip. The simulations take both the interatomic interaction and the interaction with the surface of the nanotip into account. The dipole–dipole interaction strength δ₁₂ and the collective correction γ₁₂ are determined using realistic parameters. The corresponding g⁽²⁾ function is calculated from a large number of random atomic positions.

The physical mechanism described in this paper can be applied to any statistical ensemble of emitters in a nanoscale volume. However, the specific system which we have in mind is that of an optical fiber which is tapered down at one of its ends to a nanotip in order to increase interactions with dipole emitters nearby the tip [13–16]. We assume the tip is immersed into a gas of cold atoms with number density ρ, and a classical light field emerging from the fiber end excites atoms in a small volume V in front of its apex, figure 1. Thus, an average number of

$$\begin{eqnarray}&&\mu =\rho V\end{eqnarray} \tag{ 1 }$$

atoms is within the volume. The exact number N of atoms in the volume is poissonianly distributed. i.e. the probability of finding N atoms in the volume is given by

$$\begin{eqnarray}&&{P}_{\mu }(N)=\displaystyle \frac{{\mu }^{N}}{N!}{{\rm{e}}}^{-\mu }.\end{eqnarray} \tag{ 2 }$$

We will first show how the light statistics is influenced by such a poissonianly distributed number of non-interacting emitters. The g⁽²⁾(τ)-function

$$\begin{eqnarray}&&{g}^{(2)}(\tau )=\displaystyle \frac{\left\langle {\hat{a}}^{\dagger }(t){\hat{a}}^{\dagger }(t+\tau )\hat{a}(t+\tau )\hat{a}(t)\right\rangle }{{\left\langle {\hat{a}}^{\dagger }(t)\hat{a}(t)\right\rangle }^{2}}\end{eqnarray} \tag{ 3 }$$

is proportional to the probability density of the emission of a photon at a time difference τ after the emission of a first photon. The total photon annihilation operater $\hat{a}$ refers to a single optical mode which interacts with the N atoms. We assume that the interatomic distance is well below the optical wavelength such that propogation effects between photons emitted from distant atoms can be neglected. The g⁽²⁾-function is scaled such that ${g}^{(2)}(\tau )\to 1$ for sufficiently large τ where no correlations between emission processes exist. Antibunching, which is the signature of single photon emission, corresponds to a value of g⁽²⁾(0) = 0. For a single atom the emission is antibunched due to the time it takes to re-excite the atom once it is projected to its ground state after the emission of a photon. For an arbitrary fixed number of N independent atoms the value of ${g}_{N}^{(2)}(0)$ scales with the atom number as

$$\begin{eqnarray}&&{g}_{N}^{(2)}(\tau =0)=\displaystyle \frac{N-1}{N}.\end{eqnarray} \tag{ 4 }$$

This result is caused by the effect that, whenever one of the N atoms emits a photon, any of the other N − 1 atoms can emit a photon, too. Finally, the contributions of different atom numbers to the g⁽²⁾-function are weighted with their statistical probability P_μ(N) of finding N atoms in the volume and with the rate R_N at which photons are emitted by the N atoms:

$$\begin{eqnarray}&&{g}_{\mu }^{(2)}(\tau =0)=\displaystyle \frac{1}{R}\displaystyle \sum _{N=1}^{\infty }{g}_{N}^{(2)}(\tau =0){P}_{\mu }(N){R}_{N},\end{eqnarray} \tag{ 5 }$$

with normalization $R={\sum }_{N}{P}_{\mu }(N){R}_{N}$. Note that g⁽²⁾—functions of different realizations can be averaged with their relative weight, when the timescale on which the realizations are stable is larger than the correlation time; in our case i.e. the timescale on which the position of an atom changes substantially. The condition is very well fulfilled for ultracold atoms at the nanotip. With a velocity of 1 mm s⁻¹ they move over a distance of 1 nm within a time of 1 μs, which is much larger than the correlation time which is on a timescale of $1/\left(2\pi {\gamma }_{0}\right)=27\,\mathrm{ns}$ for Rubidium with natural linewidth γ₀. For N independent emitters the rate increases as R_N = N R₁, i.e. the number of emitted photons increases linearly with the number of atoms. Thus, the g⁽²⁾-function quickly tends to a value of one for increasing mean atom number μ, as shown in figure 1. Nevertheless, antibunching can be observed also with statistical ensembles of emitters, if $\mu \ll 1$, where—in turn—the probability P_μ(0) that no atom is in the volume gets high. This effect reduces the brightness B of the source, similar to other non-deterministic sources like those based on parametric down-conversion [8]. We define the brightness as total photon emission rate normalized to the photon emission rate of a single atom:

$$\begin{eqnarray}&&B(\mu )=\displaystyle \frac{1}{{R}_{1}}\displaystyle \sum _{N=1}^{\infty }{P}_{\mu }(N){R}_{N}.\end{eqnarray} \tag{ 6 }$$

For independent atoms the brightness

$$\begin{eqnarray}&&B(\mu )=\displaystyle \sum _{N=1}^{\infty }{P}_{\mu }(N)N=\mu \end{eqnarray} \tag{ 7 }$$

is equal to the average atom number μ, see figure 1. Interactions are able to reduce the emission rate R_N for atom numbers N > 1. In an extreme case, interactions would completely suppress the excitation of more than one atom, and the emission rate for a single atom would be R₁ = γ₀, whereas R_N = 0 for N > 1. The corresponding value of the correlation function would be g⁽²⁾(0) = 0 for any μ with a maximum brightness of B ≈ 0.37 for μ = 1, see dashed lines in figure 1. The statistics of this single-photon source shows perfect antibunching while the brightness is comparable to a single-emitter based photon source. The following sections deal with the influence of dipole interactions on the light statistics. We focus specifically on the case of one and two atoms. This is justified by the fact that for increasing atom number in the volume the suppression will be even stronger due to the decreasing distance between the atoms. Moreover, for a mean atom number of μ = 1 the Poisson-probability that three atoms are within the volume is only P(3) = 0.06.

Zoom In Zoom Out Reset image size

The theoretical method we use in this paper is known as the quantum jump approach to dissipative dynamics [17–20]. The atomic state $\left|{{\rm{\Psi }}}_{j}\right\rangle$ of atom number j = 1, 2 is described in the matrix basis of the ground state $\left|g\right\rangle$ and the excited state $\left|e\right\rangle$

$$\begin{eqnarray}&&\left|g\right\rangle =\left(\begin{array}{c}1\\ 0\end{array}\right),\left|e\right\rangle =\left(\begin{array}{c}0\\ 1\end{array}\right).\end{eqnarray} \tag{ 8 }$$

The states are coupled by a laser field with detuning Δ = ω_L − ω₀ between the laser frequency ω_L and the atomic transition frequency ω₀ of the atoms. The Rabi frequency of the coupling is ${{\rm{\Omega }}}_{j}={\rm{d}}{{ \mathcal E }}_{j}/{\hslash }$, with dipole matrix element d of the transition, and electric field amplitude ${{ \mathcal E }}_{j}$ of the light at the position of the corresponding atom. In the rotating frame, the effective single-atom Hamiltonian reads

$$\begin{eqnarray}{H}_{j}={\hslash }\left(\begin{array}{cc}0 & {{\rm{\Omega }}}_{j}/2\\ {{\rm{\Omega }}}_{j}/2 & {\rm{\Delta }}\end{array}\right).\end{eqnarray} \tag{ 9 }$$

The total wavefunction $\left|{\rm{\Psi }}\right\rangle$ describing the state of both atoms including correlations is an element of the four-dimensional tensor product Hilbert space with the vector basis

$$\begin{eqnarray}\begin{array}{rcl}\left|{gg}\right\rangle & = & \left|g\right\rangle \otimes \left|g\right\rangle ,\left|{ge}\right\rangle =\left|g\right\rangle \otimes \left|e\right\rangle ,\\ \left|{eg}\right\rangle & = & \left|e\right\rangle \otimes \left|g\right\rangle ,\left|{ee}\right\rangle =\left|e\right\rangle \otimes \left|e\right\rangle ,\end{array}\end{eqnarray} \tag{ 10 }$$

where the operator left and right to ⨂ corresponds to atom no. 1 and 2, respectively. The Hamiltonian acting on the tensor space is given by the single atom Hamiltonians (9) via

$$\begin{eqnarray}&&{H}_{{\rm{a}}}={H}_{1}\otimes {\mathbb{1}}+{\mathbb{1}}\otimes {H}_{2},\end{eqnarray} \tag{ 11 }$$

with two-dimensional identity matrix ${\mathbb{1}}$. Thus, the Hamiltonian without interactions is

$$\begin{eqnarray}{H}_{{\rm{a}}}={\hslash }\left(\begin{array}{cccc}0 & \displaystyle \frac{{{\rm{\Omega }}}_{2}}{2} & \displaystyle \frac{{{\rm{\Omega }}}_{1}}{2} & 0\\ \displaystyle \frac{{{\rm{\Omega }}}_{2}}{2} & {\rm{\Delta }} & 0 & \displaystyle \frac{{{\rm{\Omega }}}_{1}}{2}\\ \displaystyle \frac{{{\rm{\Omega }}}_{1}}{2} & 0 & {\rm{\Delta }} & \displaystyle \frac{{{\rm{\Omega }}}_{2}}{2}\\ 0 & \displaystyle \frac{{{\rm{\Omega }}}_{1}}{2} & \displaystyle \frac{{{\rm{\Omega }}}_{2}}{2} & 2{\rm{\Delta }}\end{array}\right),\end{eqnarray} \tag{ 12 }$$

Dipole–dipole interactions between the two atoms are described by the interaction Hamiltonian [21]

$$\begin{eqnarray}\begin{array}{rcl}W & = & {\hslash }{\delta }_{12}\left({S}^{+}\otimes {S}^{-}+{S}^{-}\otimes {S}^{+}\right)\\ & = & {\hslash }\left(\begin{array}{cccc}0 & 0 & 0 & 0\\ 0 & 0 & {\delta }_{12} & 0\\ 0 & {\delta }_{12} & 0 & 0\\ 0 & 0 & 0 & 0\end{array}\right),\end{array}\end{eqnarray} \tag{ 13 }$$

with interaction potential strength δ₁₂ and operators

$$\begin{eqnarray}&&{S}^{+}=\left|e\right\rangle \left\langle g\right|,\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{S}^{-}=\left|g\right\rangle \left\langle e\right|.\end{eqnarray} \tag{ 15 }$$

The complete Hamiltonian is thus given by

$$\begin{eqnarray}&&{H}_{\mathrm{tot}}={H}_{{\rm{a}}}+W.\end{eqnarray} \tag{ 16 }$$

The evolution of the density matrix is determined by the master equation in its Lindblad form [21]

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}\rho }{{\rm{d}}{t}} & = & -\displaystyle \frac{{\rm{i}}}{{\hslash }}\left[{H}_{\mathrm{tot}},\rho \right]-\displaystyle \sum _{i,j=1}^{2}{\gamma }_{{ij}}\left(\displaystyle \frac{1}{2}{S}_{i}^{+}{S}_{j}^{-}\rho \right.\\ & & +\left.\displaystyle \frac{1}{2}\rho {S}_{i}^{+}{S}_{j}^{-}-{S}_{j}^{-}\rho {S}_{i}^{+}\right).\end{array}\end{eqnarray} \tag{ 17 }$$

Here, the matrix

$$\begin{eqnarray}{\gamma }_{{ij}}=\left(\begin{array}{cc}{\gamma }_{1} & {\gamma }_{12}\\ {\gamma }_{12} & {\gamma }_{2}\end{array}\right)\end{eqnarray} \tag{ 18 }$$

is defined with decay rates γ₁ and γ₂ of the individual atoms and collective correction γ₁₂ to the decay rate which will be introduced in section 6, and

$$\begin{eqnarray}&&{S}_{1}^{\pm }={S}^{\pm }\otimes {\mathbb{1}},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{S}_{2}^{\pm }={\mathbb{1}}\otimes {S}^{\pm }.\end{eqnarray} \tag{ 20 }$$

The quantum state trajectory is simulated in time steps with length dt. For each time step the probability that a jump m = 1, 2 occurs is given by

$$\begin{eqnarray}&&{p}_{m}=\left\langle {\rm{\Psi }}\right|{L}_{m}^{\dagger }{L}_{m}\left|{\rm{\Psi }}\right\rangle {\rm{d}}{t}.\end{eqnarray} \tag{ 21 }$$

The collective jump operators

$$\begin{eqnarray}&&{L}_{1}=\sqrt{{{\rm{\Lambda }}}_{1}}\left(\alpha {S}_{1}^{-}+\beta {S}_{2}^{-}\right),\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{L}_{2}=\sqrt{{{\rm{\Lambda }}}_{2}}\left(-\beta {S}_{1}^{-}+\alpha {S}_{2}^{-}\right).\end{eqnarray} \tag{ 23 }$$

are derived following [20] as eigenvectors of the matrix (18) with eigenvalues

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{1}=\bar{\gamma }+\sqrt{{\left({\gamma }_{12}\right)}^{2}+{{\rm{\Delta }}}_{\gamma }^{2}},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{2}=\bar{\gamma }-\sqrt{{\left({\gamma }_{12}\right)}^{2}+{{\rm{\Delta }}}_{\gamma }^{2}},\end{eqnarray} \tag{ 25 }$$

where

$$\begin{eqnarray}&&\bar{\gamma }=\displaystyle \frac{1}{2}\left({\gamma }_{1}+{\gamma }_{2}\right),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\gamma }=\displaystyle \frac{1}{2}\left({\gamma }_{1}-{\gamma }_{2}\right).\end{eqnarray} \tag{ 27 }$$

They are linear combinations of the single atom jump operators ${S}_{j}^{-}$ with coefficients

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{{\rm{s}}{\rm{ign}}({\gamma }_{12})\left({{\rm{\Delta }}}_{\gamma }+\sqrt{{\left({\gamma }_{12}\right)}^{2}+{{\rm{\Delta }}}_{\gamma }^{2}}\right)}{\sqrt{{\left({\gamma }_{12}\right)}^{2}+{\left({{\rm{\Delta }}}_{\gamma }+\sqrt{{\left({\gamma }_{12}\right)}^{2}+{{\rm{\Delta }}}_{\gamma }^{2}}\right)}^{2}}},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\beta =\displaystyle \frac{\left|{\gamma }_{12}\right|}{\sqrt{{\left({\gamma }_{12}\right)}^{2}+{\left({{\rm{\Delta }}}_{\gamma }+\sqrt{{\left({\gamma }_{12}\right)}^{2}+{{\rm{\Delta }}}_{\gamma }^{2}}\right)}^{2}}}.\end{eqnarray} \tag{ 29 }$$

For equal single atom decay rates γ₁ = γ₂ = γ, the collective jump operators are given by

$$\begin{eqnarray}&&{L}_{1}=\sqrt{\displaystyle \frac{\gamma +\left|{\gamma }_{12}\right|}{2}}\left({\rm{s}}{\rm{ign}}({\gamma }_{12}){S}_{1}^{-}+{S}_{2}^{-}\right),\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{L}_{2}=\sqrt{\displaystyle \frac{\gamma -\left|{\gamma }_{12}\right|}{2}}\left(-{S}_{1}^{-}+{\rm{s}}{\rm{ign}}({\gamma }_{12}){S}_{2}^{-}\right),\end{eqnarray} \tag{ 31 }$$

corresponding to superradiant and subradiant decay, respectively. In the case considered in section 6 where the atoms may have different distance from the nanotip, the single atom decay rates may also differ from each other. Finally, the decay terms can be included into an effective Hamiltonian

$$\begin{eqnarray}&&{H}_{\mathrm{eff}}={H}_{\mathrm{tot}}-{\rm{i}}{J},\end{eqnarray} \tag{ 32 }$$

with

$$\begin{eqnarray}&&J=\displaystyle \frac{1}{2}\displaystyle \sum _{m}{L}_{m}^{\dagger }{L}_{m}.\end{eqnarray} \tag{ 33 }$$

The calculated probabilities (21) for the individual jumps are compared with a random evenly distributed number $r\in \left[0,1\right]$. If $r\lt {p}_{0}\,=\,1-({p}_{1}+{p}_{2})$ no jump occurs, and the wavefunction evolves according to Hamiltonian dynamics with H_eff. As the effective Hamiltonian is not Hermitian (due to the decay terms), the wavefunction is renormalized to its absolute value after each step. The Hamiltonian dynamics can thus be summarized as

$$\begin{eqnarray}&&\left|{\rm{\Psi }}\right\rangle \to {\left|{\rm{\Psi }}\right\rangle }^{{\prime} }=\displaystyle \frac{\exp (-{{\rm{i}}{H}}_{{\rm{e}}{ff}}{\rm{d}}{t})\left|{\rm{\Psi }}\right\rangle }{{\left|\exp (-{{\rm{i}}{H}}_{{\rm{e}}{ff}}{\rm{d}}{t})\left|{\rm{\Psi }}\right\rangle \right|}^{2}}.\end{eqnarray} \tag{ 34 }$$

Note that the value of dt has to be chosen such that ${p}_{1}+{p}_{2}\ll 1$.

If, on the other hand, a jump occurs, i.e. for r > p₀, which jump operator L_m is applied is given by the value of r compared to the smallest integer m = 1, 2 such that ${\sum }_{n\,=\,0}^{m}{p}_{n}\gt r$. This condition is a clever way of determining the jump operator depending on its probability. After the jump, the wavefunction is renormalized like in the Hermitian case. Thus, the jump dynamics is summarized as

$$\begin{eqnarray}&&\left|{\rm{\Psi }}\right\rangle \to {\left|{\rm{\Psi }}\right\rangle }^{{\prime} }=\displaystyle \frac{{L}_{m}\left|{\rm{\Psi }}\right\rangle }{\left|{L}_{m}\left|{\rm{\Psi }}\right\rangle \right|}.\end{eqnarray} \tag{ 35 }$$

Before we apply this quantum jump approach to specific situations we present in the following section a general argument why dipole interactions can lead to the suppression of a simultaneous emission of photons from both atoms. For that reason we analytically solve the Schrödinger equation for two atoms including the dipole–dipole interactions.

The Schrödinger equation with the Hamiltonian H_tot from (16) including dipole–dipole interactions (but without the decay terms) can be solved analytically. We consider here the case where Δ = 0, and Ω₁ = Ω₂ = Ω in order to keep the calculation as simple as possible. Thus, we have

$$\begin{eqnarray}{H}_{\mathrm{tot}}={\hslash }\left(\begin{array}{cccc}0 & \displaystyle \frac{1}{2}{\rm{\Omega }} & \displaystyle \frac{1}{2}{\rm{\Omega }} & 0\\ \displaystyle \frac{1}{2}{\rm{\Omega }} & 0 & {\delta }_{12} & \displaystyle \frac{1}{2}{\rm{\Omega }}\\ \displaystyle \frac{1}{2}{\rm{\Omega }} & {\delta }_{12} & 0 & \displaystyle \frac{1}{2}{\rm{\Omega }}\\ 0 & \displaystyle \frac{1}{2}{\rm{\Omega }} & \displaystyle \frac{1}{2}{\rm{\Omega }} & 0\end{array}\right).\end{eqnarray} \tag{ 36 }$$

The formal solution of the Schrödinger equation is

$$\begin{eqnarray}&&| {\rm{\Psi }}(t)\rangle ={{\rm{e}}}^{-\tfrac{{\rm{i}}}{{\hslash }}{H}_{\mathrm{tot}}t}| {\rm{\Psi }}(0)\rangle .\end{eqnarray} \tag{ 37 }$$

The matrix exponential in (37) can be solved analytically. If initially both atoms are in the ground state, i.e. ${\rm{\Psi }}(0)=\left|{gg}\right\rangle$ the projection of Ψ(t) onto the doubly excited state is given by

$$\begin{eqnarray}&&\langle {ee}| {\rm{\Psi }}(t)\rangle =\displaystyle \frac{{{\rm{\Omega }}}^{2}}{4}\left(\displaystyle \frac{3{\delta }_{12}-A}{B}{{\rm{e}}}^{\tfrac{{\rm{i}}({\delta }_{12}-A)t}{2}}+\displaystyle \frac{3{\delta }_{12}+A}{C}{{\rm{e}}}^{\tfrac{{\rm{i}}({\delta }_{12}+A)t}{2}}\right)-\displaystyle \frac{1}{2},\end{eqnarray} \tag{ 38 }$$

with parameters

$$\begin{eqnarray}&&A=\sqrt{{\delta }_{12}^{2}+4{{\rm{\Omega }}}^{2}},\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&B=\displaystyle \frac{1}{2}{\left({\delta }_{12}-A\right)}^{3}-\left({\delta }_{12}^{2}+{{\rm{\Omega }}}^{2}\right)\left({\delta }_{12}-A\right)-{\delta }_{12}{{\rm{\Omega }}}^{2},\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&C=\displaystyle \frac{1}{2}{\left({\delta }_{12}+A\right)}^{3}-\left({\delta }_{12}^{2}+{{\rm{\Omega }}}^{2}\right)\left({\delta }_{12}+A\right)-{\delta }_{12}{{\rm{\Omega }}}^{2}.\end{eqnarray} \tag{ 41 }$$

The probability of a double-excitation ${P}_{\mathrm{ee}}(t)={\left|\langle {ee}| {\rm{\Psi }}(t)\rangle \right|}^{2}$ is plotted in figure 2 for varying interaction strength δ₁₂. An increasing interaction strength slows down the timescale of the excitation process. This is the crucial feature for the suppression of the photon emission by more than one atom and leads to a suppression of the simultaneous emission of more than one photon. A further interesting feature can be observed in figure 2, when P_ee(t) is compared with the product of the probabilities that respectively one of the atoms is excited, whereas the state of the other atom is arbritrary. For ${\delta }_{12}\ne 0$ these probabilities differ from each other. This discrepancy is interpreted as interaction-induced entanglement between the two atoms and will be analyzed in a separate work, [22]. In the following section the decay term is reintroduced, and the quantum jump approach is applied in order to derive corresponding g⁽²⁾-functions.

Zoom In Zoom Out Reset image size

In this section we analyze general situations which do not depend on a specific geometry, and in particular not on the presence of a nanotip. For that reason we regard the interaction strength δ₁₂ as free parameter. Moreover, we set both decay rates to the same value γ = γ₁ = γ₂, and – first – neglect the collective decay by setting γ₁₂ = 0. The latter assumption is justified in the sense that in real situations δ₁₂ can have arbitrary large values whereas the collective correction term is always limited to the range ${\gamma }_{12}\in \left[-\gamma ,\gamma \right]$. Time trajectories of the Hamiltonian (16) including the interaction terms are simulated using the quantum jump approach explained in section 3. In order to take the Poisson-distribution of atom numbers into account, the relative length of the corresponding trajectories is adjusted accordingly. We chose here an average atom number of μ = 1. With P₁(1) = 0.37 and P₁(2) = 0.18 the trajectory of the one-atom case is twice as long as that of the two atom case. The no-atom case with P₁(0) = 0.37 is not considered as it does only reduce the brightness, but leaves the g⁽²⁾-function unchanged. Cases with more than two atoms are also not considered, because their contribution is small. Each simulation results in a series of time-stamps when photons are emitted corresponding to the times a jump occurs in the simulation. Finally, the time trajectories for the one-atom and two-atom case are appended to each other, and from the resulting trajectory the g⁽²⁾-function is calculated. Corresponding g⁽²⁾-functions derived from such simulated time trajectories are exemplarily shown in figure 3. From time trajectories as shown in figure 3 the value of g⁽²⁾(0) is evaluated and plotted in figure 4. Without interactions g⁽²⁾(0) = 0.34 coincides in good approximation with the value of a Poissonian statistical mixture of N atoms as derived in equation (5). This justifies the validity of the approximation of neglecting emission from more than two atoms for the present case of μ = 1. For increasing interaction strength the value of g⁽²⁾(0) tends to zero, and the emission resembles more and more that of asingle atom. The brightness however does not tend to zero, but saturates at a value given by the fraction of the duration of the single-atom trajectory as compared to the total time. This corresponds to the fact that single photons are emitted with the single-emitter rate whenever a single atom is in the volume whereas at most one photon is emitted at times when two atoms are in the volume. In order to take the effect of collective decay into account, which was neglected in the main part of figure 4, we simulated time trajectories for fixed values of δ₁₂ and varying γ₁₂. The result is plotted in the inset of figure 4 for δ₁₂ = 5γ. We observe that the value of g⁽²⁾(0) is increased for ${\gamma }_{12}\to \pm \gamma$. This effect is more pronounced for small interaction strength δ₁₂. For large values of δ₁₂ the influence of γ₁₂ on the g⁽²⁾-function disappears. Concluding, photon antibunching can be reached in a statistical mixture of N interacting atoms without the drawback of a low brightness which is a known issue in standard non-deterministic single photon sources. The following section analyzes the specific situation of atoms that are excited by the light field emerging from a tapered nanotip as sketched in figure 1. Collective dispersive and radiative effects of the two atoms as well as the dispersive and radiative influence of the surface on the individual atoms are taken explicitely into account.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We consider two quantum emitters of transition frequency ω₀ = k₀c located at positions ${{\bf{r}}}_{1}={r}_{1}{\hat{{\bf{r}}}}_{1}$ and ${{\bf{r}}}_{2}={r}_{2}{\hat{{\bf{r}}}}_{2}$ in the vicinity of a spherical nanotip of radius R_tip (r₁, r₂ > R_tip), dielectric permittivity ε, and centered at the origin. The interatomic separation is ${r}_{12}=| {{\bf{r}}}_{1}-{{\bf{r}}}_{2}|$, as shown in figures 5(A) and (B). For large r₁₂, i.e. ${k}_{0}{r}_{12}\gg 1$, the collective correction to the decay rate and dipole–dipole interaction strength are negligible for any mutual orientation of the dipole moments, γ₁₂ = δ₁₂ = 0 . Conversely, in the limit of separations r₁₂ much smaller than the optical wavelength, i.e. ${k}_{0}{r}_{12}\ll 1$, the collective correction and the dipole–dipole interaction strength are given by the approximate expressions [23]

$$\begin{eqnarray}&&{\gamma }_{12}\approx \sqrt{{\gamma }_{1}{\gamma }_{2}}\left({\hat{{\bf{d}}}}_{1}\cdot {\hat{{\bf{d}}}}_{2}\right),\end{eqnarray} \tag{ 42 }$$

and

$$\begin{eqnarray}&&{\delta }_{12}\approx 3\sqrt{{\gamma }_{1}{\gamma }_{2}}\left[\left({\hat{{\bf{d}}}}_{1}\cdot {\hat{{\bf{d}}}}_{2}\right)-3\left({\hat{{\bf{d}}}}_{1}\cdot {\hat{{\bf{r}}}}_{12}\right)\left({\hat{{\bf{d}}}}_{2}\cdot {\hat{{\bf{r}}}}_{12}\right)\right]/4{\left({k}_{0}{r}_{12}\right)}^{3},\end{eqnarray} \tag{ 43 }$$

respectively, where ${{\bf{d}}}_{j}={d}_{{eg}}{\hat{{\bf{d}}}}_{j}$ and ${{\bf{r}}}_{12}={r}_{12}{\hat{{\bf{r}}}}_{12}$ with the dipole matrix element d_eg [23–26]. Our simulations are based on the general expression valid for arbitrary separations r₁₂ [23, 24]

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\gamma }_{12}}{\sqrt{{\gamma }_{1}{\gamma }_{2}}} & = & \displaystyle \frac{3}{2}\left\{\left[{\hat{{\bf{d}}}}_{1}\cdot {\hat{{\bf{d}}}}_{2}-\left({\hat{{\bf{d}}}}_{1}\cdot {\hat{{\bf{r}}}}_{12}\right)\left({\hat{{\bf{d}}}}_{2}\cdot {\hat{{\bf{r}}}}_{12}\right)\right]\displaystyle \frac{\sin \left({k}_{0}{r}_{12}\right)}{{k}_{0}{r}_{12}}\right.\\ & & \left.+\left[\left({\hat{{\bf{d}}}}_{1}\cdot {\hat{{\bf{d}}}}_{2}\right)-3\left({\hat{{\bf{d}}}}_{1}\cdot {\hat{{\bf{r}}}}_{12}\right)\left({\hat{{\bf{d}}}}_{2}\cdot {\hat{{\bf{r}}}}_{12}\right)\right]\times \left[\displaystyle \frac{\cos \left({k}_{0}{r}_{12}\right)}{{\left({k}_{0}{r}_{12}\right)}^{2}}-\displaystyle \frac{\sin \left({k}_{0}{r}_{12}\right)}{{\left({k}_{0}{r}_{12}\right)}^{3}}\right]\right\},\end{array}\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\delta }_{12}}{\sqrt{{\gamma }_{1}{\gamma }_{2}}} & = & \displaystyle \frac{3}{4}\left\{-\left[{\hat{{\bf{d}}}}_{1}\cdot {\hat{{\bf{d}}}}_{2}-\left({\hat{{\bf{d}}}}_{1}\cdot {\hat{{\bf{r}}}}_{12}\right)\left({\hat{{\bf{d}}}}_{2}\cdot {\hat{{\bf{r}}}}_{12}\right)\right]\displaystyle \frac{\cos \left({k}_{0}{r}_{12}\right)}{{k}_{0}{r}_{12}}\right.\\ & & \left.+\left[\left({\hat{{\bf{d}}}}_{1}\cdot {\hat{{\bf{d}}}}_{2}\right)-3\left({\hat{{\bf{d}}}}_{1}\cdot {\hat{{\bf{r}}}}_{12}\right)\left({\hat{{\bf{d}}}}_{2}\cdot {\hat{{\bf{r}}}}_{12}\right)\right]\times \left[\displaystyle \frac{\sin \left({k}_{0}{r}_{12}\right)}{{\left({k}_{0}{r}_{12}\right)}^{2}}+\displaystyle \frac{\cos \left({k}_{0}{r}_{12}\right)}{{\left({k}_{0}{r}_{12}\right)}^{3}}\right]\right\}.\end{array}\end{eqnarray} \tag{ 45 }$$

Equations (44) and (45) are well-known expressions for the collective correction to the decay rate and dipole–dipole interaction strength of two dipole emitters in vacuum [27]. However, here we are interested in the influence of a spherical nanotip on the two-atom dynamics, which is expected to be modified not only by the collective damping but also by the Purcell effect [28]. To this end, we use equations (44) and (45) and substitute the free-space decay rates by the decay rates γ₁ and γ₂ of the corresponding atoms not in free space but in the vicinity of the spherical nanoparticle. This means that the decay-rate enhancement induced by the dielectric nanofiber tip enters into equations (44) and (45) as a perturbation in the single-atom decay rates., i.e. ${\gamma }_{j}={\gamma }_{0}+{\gamma }_{j}^{\mathrm{ref}}$, where ${\gamma }_{j}^{\mathrm{ref}}$ is the reflected-field contribution associated with atom j = 1, 2. Since ${\gamma }_{12}\propto \sqrt{{\gamma }_{1}{\gamma }_{2}}$, the approximate collective parameters take into account both the emitted-field contribution from each atom and the overall reflected-field contribution from the nanofiber tip. We have verified that this approximation provides accurate results when compared with the exact solution as long as we consider moderate permittivities with Im(ε) = 0 and ε > 0, i.e. there is no ohmic losses on the dielectric nanotip. To simplify our discussion, we also consider that the spherical nanotip is much smaller than the emission wavelength, i.e. k₀R_tip≪1. Under this last assumption, the single-atom decay rate is given by [29]

$$\begin{eqnarray}&&{\gamma }_{j}={\gamma }_{j}^{\perp }{\left({\hat{{\bf{d}}}}_{j}\cdot {\hat{{\bf{r}}}}_{j}\right)}^{2}+{\gamma }_{j}^{| | }\left[1-{\left({\hat{{\bf{d}}}}_{j}\cdot {\hat{{\bf{r}}}}_{j}\right)}^{2}\right],\end{eqnarray} \tag{ 46 }$$

where

$$\begin{eqnarray}&&\displaystyle \frac{{\gamma }_{j}^{\perp }}{{\gamma }_{0}}\approx 1+\displaystyle \frac{3{k}_{0}^{3}}{2\pi }\mathrm{Im}\left\{\alpha ({\omega }_{0}){{\rm{e}}}^{2{{\rm{i}}{k}}_{0}{r}_{j}}\left[-\displaystyle \frac{1}{{\left({k}_{0}{r}_{j}\right)}^{4}}+\displaystyle \frac{2}{{\rm{i}}{\left({k}_{0}{r}_{j}\right)}^{5}}+\displaystyle \frac{1}{{\left({k}_{0}{r}_{j}\right)}^{6}}\right]\right\},\end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\gamma }_{j}^{| | }}{{\gamma }_{0}}\approx 1+\displaystyle \frac{3{k}_{0}^{3}}{8\pi }\mathrm{Im}\left\{\alpha ({\omega }_{0}){{\rm{e}}}^{2{{\rm{i}}{k}}_{0}{r}_{j}}\left[\displaystyle \frac{1}{{\left({k}_{0}{r}_{j}\right)}^{2}}-\displaystyle \frac{2}{{\rm{i}}{\left({k}_{0}{r}_{j}\right)}^{3}}-\displaystyle \frac{3}{{\left({k}_{0}{r}_{j}\right)}^{4}}+\displaystyle \frac{2}{{\rm{i}}{\left({k}_{0}{r}_{j}\right)}^{5}}+\displaystyle \frac{1}{{\left({k}_{0}{r}_{j}\right)}^{6}}\right]\right\},\end{eqnarray} \tag{ 48 }$$

with ${\gamma }_{0}={d}_{{eg}}^{2}{\omega }_{0}^{3}/3\pi {\varepsilon }_{0}{\hslash }{c}^{3}$ the Einstein A coefficient.

Zoom In Zoom Out Reset image size

We have verified that equations (47) and (48) are very good approximations for the exact expressions in [30] when the sphere is subwavelength. The polarizability α(ω) of a subwavelength nanosphere, accounting for radiation damping, is [29]

$$\begin{eqnarray}&&\alpha (\omega )=\displaystyle \frac{{\alpha }_{0}}{1-{\rm{i}}({k}^{3}/6\pi ){\alpha }_{0}},\end{eqnarray} \tag{ 49 }$$

where ${\alpha }_{0}\,=\,4\pi {R}_{\mathrm{tip}}^{3}(\varepsilon -1)/(\varepsilon +2)$ is the quasi-static polarizability of the nanosphere. The transmitted/scattered electric field by the spherical nanotip determines the dipole moment orientation. Within the fiber, we consider a z-polarized electromagnetic plane wave with wave vector ${{\bf{k}}}_{L}={k}_{L}\hat{{\bf{y}}}$ and amplitude E₀, where ω_L = k_Lc is the laser frequency. Since we work with light tuned close to the atomic resonance, we use ${k}_{L}\cong {k}_{0}$. Outside the spherical nanotip, the local electric field in usual spherical coordinates r = (r, θ, φ) is given by [31, p 411] as

$$\begin{eqnarray}&&{\bf{E}}({\bf{r}},{\omega }_{L})={E}_{0}\displaystyle \frac{\alpha ({\omega }_{L})}{4\pi {r}^{3}}{{\rm{e}}}^{{{\rm{i}}{k}}_{0}r}\left\{2\cos \theta \left(1-{{\rm{i}}{k}}_{0}r\right)\hat{{\bf{r}}}+\sin \theta \left[1-{{\rm{i}}{k}}_{0}r-{\left({k}_{0}r\right)}^{2}\right]\hat{{\boldsymbol{\theta }}}\right\},\end{eqnarray} \tag{ 50 }$$

where α(ω_L) is the effective polarizability of the spherical nanotip given in equation (49). Here, the near- and far-field contributions are proportional to 1/r³ and ${k}_{0}^{2}/r$, respectively, whereas the intermediate region are mainly governed by the term proportional to ik₀/r². Note that the electric field fixes the spherical coordinate system. In this case, the y axis is directed along the fiber and parallel to the horizontal axis of figure 1. This choice of coordinates implies E_φ = 0, which simplifies the expressions. Finally, we consider that the dipole moment of atom j is directed along the local electric field at r_j [28], i.e. ${{\bf{d}}}_{j}={d}_{{eg}}\{{\bf{E}}({{\bf{r}}}_{j},{\omega }_{L})\}/| | \{{\bf{E}}({{\bf{r}}}_{j},{\omega }_{L})\}| |$. This leads to the Rabi frequency

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{j}({{\bf{r}}}_{j},{\omega }_{L})=\displaystyle \frac{{d}_{{eg}}}{{\hslash }}| | \left\{{\bf{E}}({{\bf{r}}}_{j},{\omega }_{L})\right\}| | ,\end{eqnarray} \tag{ 51 }$$

where the electric field is given in equation (50). Note that we drop the time-harmonic dependence ${{\rm{e}}}^{-{\imath }{\omega }_{L}t}$ in equation (50), so that equation (51) is the usual definition of the Rabi frequency. Two specific geometries are considered in figures 5(A) and (B). In these geometries the interatomic distance is kept fixed at r₁₂ = 50 nm, and the position of the atoms relative to the tip is varied. The two geometries correspond to situations where the interatomic vector ${\hat{{\bf{r}}}}_{12}$ is parallel to the y- and z-axis, respectively. The results for the collective correction and the dipole–dipole interaction strength are plotted in figures 5(A/B-i/ii). The interaction strength is large, i.e. $\left|{\delta }_{12}\right|\gg {\gamma }_{0}$, in a large region around the tip, such that a reduction of the value of g⁽²⁾(0) can be expected according to figure 3. Note the fact that the collective parameters δ₁₂ and γ₁₂ are modulated as a function of the angle θ even when the atomic distance r from the tip is large. This is caused by the change in the relative orientation of the two dipoles following the direction of the electric field emitted from the tip.

In order to predict the result of a real experiment we average over many arbitrarily chosen atomic positions. For this purpose atoms are randomly placed around the tip in a half-sphere shell in 3D with inner radius given by the tip size, R_i = R_tip = 100 nm, and outer radius twice as large, R_o = 200 nm. Its volume is V = 2.9 × 10⁻¹⁴ cm³. Thus, a mean number of μ = 1 atom in the volume corresponds to an atomic density of ρ = 3.4 × 10¹³ cm⁻³ which can be reached in experiments with ultracold atoms. In order to take into account the poissonian distribution of atoms in the volume for μ = 1 where the probability of having one atom in the volume is twice as large as having two atoms in the volume, twice as many realizations are simulated for the one-atom case than for the two-atom case. We simulate 100 and 50 realizations, respectively. For each realization a quantum Monte Carlo simulation is carried out with the respective parameters γ₁, γ₂, γ₁₂, δ₁₂, Ω₁, and Ω₂ given by the positions of the atoms relative to the tip. The Rabi frequencies are an open parameter, because they are proportional to the laser power used in the experiment. In the simulation they are scaled such that a Rabi frequency of Ω = γ₀ is reached at the surface of the nanotip, i.e. at the position r = 100 nm, θ = 90°, and φ = 90° within the yz-plane. Moreover, in the one-atom case, the Rabi frequency acting on the second atom and the collective parameters are artificially set to zero, Ω₂ = 0, δ₁₂ = 0, and γ₁₂ = 0. Thus, the second atom is not excited and has got no influence on the first atom. A g⁽²⁾-function is calculated from the trajectory of each QMC-simulation. The individual g⁽²⁾-functions are weighted with their corresponding rate of photon emission and are averaged. Such averages are calculated for the one-atom cases and the two-atom cases, and for the statistical mixture of both. The results are plotted in figure 6. The averaged g⁽²⁾-function of the one-atom cases drops to zero for $\tau \to 0$, as expected for a single emitter. In the two-atom cases the value of the averaged g⁽²⁾(0) > 0.5. The dipole–dipole interaction of both atoms leads on average to an increase of g⁽²⁾(0) above the non-interacting case and thus to bunching, although values of ${g}^{(2)}(0)\ll 0.5$ (i.e. strong antibunching) can be reached for specific positions of the two atoms. However, the photon emission rate of the two-atom case is reduced, compared to the non-interacting situation. The averaged photon emission rates in figure 6 are given in the one-atom case by R_1A = 0. 30 γ₀, in the two-atom case by ${R}_{2{\rm{A}}}=0.42\,{\gamma }_{0}$, and in the statistical mixture by R_mix = 0. 34 γ₀. Without interactions the emission rate of the two-atom case would be ${R}_{2{\rm{A}}}=2{R}_{1{\rm{A}}}=0.6\,{\gamma }_{0}$. Thus, the value of g⁽²⁾(0) of the statistical mixture of atoms is reduced to a value of g⁽²⁾(0) = 0.25 smaller than without interactions. At the same time, the brightness B defined by (6) is with $B=\tfrac{{R}_{\mathrm{mix}}}{{R}_{1{\rm{A}}}}=1.13$ comparable to the one-atom case.

Zoom In Zoom Out Reset image size

In the simulations shown in figure 6 we have neglected the action of the Casimir–Polder (CP) potential on the atoms. Depending on the distance from the surface this potential changes the atomic transition frequencies and leads to an effective detuning. The strength of the CP potential close to a dieletric sphere can be calculated following [32]. In our situation, i.e. for Rubidium atoms at the surface of a dielectric nanosphere with refractive index n = 1.5 and sphere radius R = 100 nm, the CP-shift is comparable with the single atom decay rate γ_j from (46) at a distance of ∼40 nm. For shorter distances the atomic transition is tuned out of resonance with the exciting light field. Including the CP potential into our simulations we have checked that photon emission from these atoms is negligible. Furthermore, it has turned out that in specific two-atom cases when one of the atoms is detuned and the second is not (because its distance is larger), the g⁽²⁾-function resembles that of a single atom. Averaging over many arbitrary geometries the value of g⁽²⁾(0) is however only slightly smaller than that shown in figure 6 without the CP-shift, which is why we do not show these g⁽²⁾-functions.

The paper investigates the light statistics of the fluorescence from statistical mixtures of one and two atomic emitters within a nanoscale volume. The ${g}^{(2)}$-function of the emitted light is derived by a quantum Monte Carlo approach including dipole–dipole interactions between the atoms and collective decay. The main result of the paper is the fact that these interactions can improve the antibunching while the brightness of the source tends to a constant value. This result is interesting for the design of novel single photon sources that are not based on single emitters but on statistical ensembles. As application, we investigate the situation of cold atoms in a nanoscale volume at the apex of a nanofiber tip and predict a reduction of the g⁽²⁾-function due to atomic interactions to a value of g⁽²⁾(0) = 0.25. We expect that the quality of antibunching can be further improved by optimizing the radius of the nanotip, the atomic density and the light power used. Furthermore, manipulating the shape of the optical near-field, for instance by depositing tailored plasmonic nanostructures on the tip, may be a method to restrict the excitation of atoms to those positions where the antibunching is more pronounced. A second possibility not explored so far in order to manipulate the excitation zone would be to exploit the resonance shift of the atoms caused by the CP-potential and compensate it in a specific distance from the surface by the proper detuning. The fluorescence may be also collected by the nanotip which puts further restrictions on the atomic positions that contribute to the detected signal [14]. Using Purcell enhancement close to plasmonic nanotips a collection efficiency of more than 95% can be reached [13]. As a final remark we would like to add that the method we introduce can be used also for other types of interactions, for instance van-der Waals interactions between Rydberg atoms [33].

This work was performed in the context of the European COST Action MP1403 Nanoscale Quantum Optics. SS acknowledges support by the Fulbright-Cottrell Award. ES acknowledges support by the ColOpt—EU H2020 ITN 721465. TJA, RB, and PhWC acknowledge FAPESP (Grant Nos. 2015/21194-3, 2014/01491-0, and 2013/04162-5, respectively) for financial support. We acknowledge support by Deutsche Forschungsgemeinschaft and Open Access Publishing Fund of University of Tübingen. Finally, we would like to thank SY Buhmann and R Bennett for support on the Casimir–Polder potential at nanospheres.

Competing financial interests do not exist.