Advances in quantum technologies over the past two decades have made possible an exciting breadth of applications in fields such as communications [1], imaging [2], and computation [3]. Photon pair generation based on the spontaneous parametric downconversion (SPDC) [4] and four-wave mixing [5] processes has played an essential role in this revolution due to the ease with which the quantum entanglement characteristics of the emitted signal and idler photons may be tailored, and on account of their ability to propagate long distances either in free space or in optical fibers, with minimal interaction with the environment. However, a number of key challenges must be overcome in order for photon pair generation technology to achieve its true potential, including: (i) source miniaturization, enabling the eventual on-chip integration of source, optical manipulation, and detection [6,7]; (ii) increasing the conversion efficiency, thus permitting high-brightness photon pair emission with the lowest possible pump power; (iii) photon pair indistinguishability, including spectral factorizability, permitting photons from distinct sources to interfere [8]; and (iv) the reduction of the emission bandwidth to the megahertz, or submegahertz, level so as to be compatible with atomic electronic transitions [9]. The last requirement must be achieved in order to create single atom–single photon interfaces that will facilitate information transfer from photons, in the form of flying qubits, to atoms, thus constituting a quantum memory. This technology is essential for the further progress of quantum information processing based on photons.

An optical platform that has the potential to simultaneously meet all of these requirements is the ultrahigh $Q$ optical microresonator. As a result of the long photon lifetimes inside the cavities, extremely high circulating intensities are possible, resulting in very low optical thresholds for nonlinear behaviors [10–12]. In addition, the extremely narrowband resonances lead to nonlinear optical effects occurring at very well-defined frequencies. Because cavity-enhanced photon pair generation involves emission in well-defined cavity modes, by isolating a single cavity resonance for each of the signal and idler modes, the resulting two-photon state ends up being naturally spectrally factorizable, ensuring indistinguishability [13]. Lastly, microresonators are uniquely well positioned to permit source integration [14,15].

There are two possible routes for generating photon pairs based on spontaneous parametric processes that can occur in microresonators: the SPDC process based on second-order nonlinear materials [16–21], and the spontaneous four-wave mixing (SFWM) process based on third-order nonlinear materials [22–29]. Photon pair generation from a cavity-enhanced source, in which the nonlinear medium is contained by an optical cavity, based on either of these two processes with a narrowband pump, leads to a joint spectral amplitude in the form of a frequency comb expressed as a function of the frequency difference ${\omega }_s-{\omega }_i$, in terms of the signal ${\omega }_s$ and idler ${\omega }_i$ frequencies. While each resulting comb peak has a width that is inversely proportional to the cavity quality factor $Q$, the intrapeak separation [or free spectral range (FSR)] is inversely proportional to the cavity round-trip time [13,30]. Thus, an important benefit of a microresonator as compared to an extended cavity design is that the comb peaks end up being separated by a greater spectral distance, facilitating the possibility of addressing individual comb peaks. Note that in the classical nonlinear optics realm, four-wave mixing in optical resonators is likewise known to naturally give rise to the emission of frequency combs, with applications in the field of atomic clocks and generally in metrology [31,32].

Note that in the standard SPDC and SFWM processes (without the use of optical cavities), the spread of emission frequencies is limited mainly by phase-matching constraints and can be substantial. Such large bandwidths are useful in certain situations, e.g., they lead to narrow Hong–Ou–Mandel interference dips, which in turn permit a large resolution in quantum optical coherence tomography devices [33,34]. However, as already mentioned, the basic requirement for the development of single atom–single photon interfaces is the emission of narrowband photon pairs [35–37].

Cavity-enhanced SPDC sources have been demonstrated using integrated microresonators [17,19], nonlinear waveguide cavities [21], and free-space extended cavities [20,38–41]. In the case of SFWM, cavity-enhanced photon pair sources have likewise been based on microring [14,15,24,25,27–29,42] and microdisk [43] cavities. In this paper, we report the first demonstration of a photon pair source based on fused silica microspheres and present a theory for SFWM in these devices that leads to simulations that agree well with our measurements. In this work, we extend our previous cavity-enhanced SFWM theory [13], so as to include important characteristics relevant to our current experimental results such as: (i) all four waves at resonance in the cavity; (ii) pump varied in time to maintain resonance, leading to two-photon state in the form of a statistical mixture; and (iii) analysis carried out for micro- rather than extended cavities. Here we report ultranarrow photon pair generation, with emission bandwidths down to $\mathrm{\Delta }\nu =366\mathrm{kHz}$. While this small single-photon bandwidth is similar to those observed in certain extended-cavity SPDC sources (e.g., 666 kHz in Ref. [38] and 265 kHz in Ref. [9]), it represents an $\sim 43\times$ improvement with respect to the previously reported spectrally narrowest SFWM source (15.9 MHz) [44].

Here we are interested in studying photon pairs produced by SFWM in a fused silica microsphere, with radius $R$, as the nonlinear medium. The pump is assumed to be coupled evanescently from an elongated (tapered) fiber to a mode circulating on the sphere’s equator; see, for example, Ref. [45]. Photon pairs produced in the sphere can then couple out back to the tapered fiber, whence they may be directed as desired to an experiment of interest; see Fig. 1.

Each of the four waves participating in the SFWM process is assumed to propagate along the equator on the plane $\theta =\pi /2$, parallel to the unit vector ${\overrightarrow{e}}_{\varphi }$. Expressed in spherical coordinates, the Hamiltonian for this process may be written as $$\begin{gathered} \widehat{H}(t)=\frac{3}{4}{ϵ}_0{\chi }^{(3)}{\int }_0^{2\pi }\mathrm{d}\varphi {\int }_0^{\pi }\mathrm{d}\theta {\int }_0^R\mathrm{d}\rho {\rho }^2\sin \theta E_1^{(+)}(\overrightarrow{r},t)E_2^{(+)}(\overrightarrow{r},t) \\ \times {\widehat{E}}_s^{(-)}(\overrightarrow{r},t){\widehat{E}}_i^{(-)}(\overrightarrow{r},t)+\text{H.C.}, \end{gathered}$$(1)in terms of the third-order optical nonlinearity ${\chi }^{(3)}$, the permittivity of free-space ${ϵ}_0$, the positive frequency electric field operators for each of the two pumps $E_1^{(+)}(\overrightarrow{r},t)$ and $E_2^{(+)}(\overrightarrow{r},t)$, and the negative frequency electric field operators for the signal and the idler, ${\widehat{E}}_s^{(-)}(\overrightarrow{r},t)$ and ${\widehat{E}}_i^{(-)}(\overrightarrow{r},t)$.

We describe each mode supported by the sphere by an index vector $\overrightarrow{l}$ composed of the $l$, $m$, and $q$ individual indices related to azimuthal, polar, and radial coordinates, respectively. Each of the pump waves ($\nu =1$ and $\nu =2$), traveling in modes ${\overrightarrow{l}}_1$ and $\overrightarrow{l_2}$, respectively, is described classically as $$E_{\nu }^{(+)}(\overrightarrow{r},t)=\sum _{{\overrightarrow{l}}_{\nu }}{A}_{{\overrightarrow{l}}_{\nu }}{f}_{{\overrightarrow{l}}_{\nu }}(\rho ,\theta )\int \mathrm{d}\omega {\alpha }_{{\overrightarrow{l}}_{\nu }}(\omega )e^{-i[\omega t-k_{{\overrightarrow{l}}_{\nu }}(\omega )R\varphi ]},$$(2)where $A_{\overrightarrow{l_{\nu }}}$ is the amplitude, $f_{\overrightarrow{l_{\nu }}}(\rho ,\theta )$ is the transverse field distribution (for a fixed value of $\varphi$), and ${\alpha }_{{\overrightarrow{l}}_{\nu }}(\omega )$ is the spectral envelope for each of the two pump waves, $\nu =1$ and $\nu =2$. The electric field for the signal and idler photons is described quantum mechanically as follows: $${\widehat{E}}_{\nu }^{(+)}(\overrightarrow{r},t)=\sum _{\overrightarrow{l_{\nu }}}{f}_{\overrightarrow{l_{\nu }}}(\rho ,\theta )\int \mathrm{d}\omega \ell (\omega )e^{-i[\omega t-k_{\overrightarrow{l_{\nu }}}(\omega )R\varphi ]}{\widehat{b}}_{\overrightarrow{l_{\nu }}}(k_{\omega }),$$(3)where ${\widehat{b}}_{\overrightarrow{l_{\nu }}}(k_{\omega })$ is the annihilation operator for mode $\overrightarrow{l_{\nu }}$, for each of the signal ($\nu =s$) and idler ($\nu =i$); $f_{\overrightarrow{l_{\nu }}}(\rho ,\theta )$ represents the signal and idler transverse spatial distributions; and the function $\ell (\omega )$ is given as $\ell (\omega )={\omega }^{1/2}/[n(\omega )v_g(\omega )]$, with $n(\omega )$ and $v_g(\omega )$ the refractive index and the group velocity, respectively.

The resonant wavelengths in the sphere can be obtained solving numerically the following equation, derived from a Mie scattering approach, for $\lambda$ (for particular values of the indices azimuthal $l$, polar $m$, and radial $q$) [46]: $$\begin{gathered} \frac{1}{\lambda }=\frac{1}{2\pi R{n}_s(\lambda )}\{\nu +2^{-\frac{1}{3}}{\alpha }_q{\nu }^{\frac{1}{3}}-\frac{P(\lambda )}{{[n_s{(\lambda )}^2-1]}^{\frac{1}{2}}} \\ +\frac{3}{10}{2}^{-\frac{2}{3}}{\alpha }_q^2{\nu }^{-\frac{1}{3}}-\frac{2^{-\frac{1}{3}}P(\lambda )[n_s{(\lambda )}^2-\frac{2}{3}P{(\lambda )}^2]}{{[n_s^2(\lambda )-1]}^{\frac{3}{2}}}{\alpha }_q{\nu }^{-\frac{2}{3}}\}, \end{gathered}$$(4)where $R$ is the sphere radius, $P(\lambda )=n_s(\lambda )$ for a TE mode, and $P(\lambda )=1/n_s(\lambda )$ for a TM mode, with $n_s(\lambda )$ the fused silica index of refraction obtained from the Sellmeier equation [47]. In Eq. (4), $\nu =l+1/2$ and ${\alpha }_q$ represents the zeros of the Airy function $\mathrm{Ai}(-z)$; for simplicity, in this work we have considered only modes with $l=\mathrm{m}$ and $q=1$, so that $\overrightarrow{l}=(l,\mathrm{m},q)$ reduces to $(l,l,1)$. We can obtain the effective index for each whispering gallery mode $n_{\mathrm{eff}}(\lambda )$, with help of the resonance condition, as $n_{\mathrm{eff}}(\lambda )=l\lambda /L$, with $L=2\pi R$. The dispersion relation is subsequently obtained as $k(\omega )=n_{\mathrm{eff}}(\omega )\omega /c$.

The quantum state is then obtained, following a standard perturbative approach, as $$|\mathrm{\Psi }⟩\approx |0⟩+\frac{1}{i\hslash }{\int }_0^t\mathrm{d}{t}^{\prime }\widehat{H}(t^{\prime })|0⟩.$$(5)

Replacing the expressions for the electric field [Eqs. (2) and (3)] into the Hamiltonian Eq. (1), we obtain the quantum state $|\mathrm{\Psi }⟩=|0⟩+\eta |{\mathrm{\Psi }}_2⟩$, written in terms of a constant related to the source brightness $\eta$ and the two-photon component of the state $|{\mathrm{\Psi }}_2⟩$, the latter given by $$\begin{gathered} |{\mathrm{\Psi }}_2⟩=\frac{1}{i\hslash }{\int }_0^t\mathrm{d}{t}^{\prime }\widehat{H}(t^{\prime })|0⟩ \\ =\sum _{{\overrightarrow{l}}_1,{\overrightarrow{l}}_2,{\overrightarrow{l}}_s,{\overrightarrow{l}}_i}{A}_{{\overrightarrow{l}}_1}{A}_{{\overrightarrow{l}}_2}{\mathrm{\Theta }}_{{\overrightarrow{l}}_1{\overrightarrow{l}}_2{\overrightarrow{l}}_s{\overrightarrow{l}}_i}\int \mathrm{d}{\omega }_s\int \mathrm{d}{\omega }_i\{l({\omega }_s)l({\omega }_i) \\ \times \int \mathrm{d}{\omega }_1[{\alpha }_{{\overrightarrow{l}}_1}({\omega }_1){\alpha }_{{\overrightarrow{l}}_2}({\omega }_s+{\omega }_i-{\omega }_1)G_{{\overrightarrow{l}}_1{\overrightarrow{l}}_2{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_1,{\omega }_s,{\omega }_i)\times {\widehat{b}}_{{\overrightarrow{l}}_s}^{†}({\omega }_s){\widehat{b}}_{{\overrightarrow{l}}_i}^{†}({\omega }_i)|0⟩]\}. \end{gathered}$$(6)

Here we have assumed that the time interval between photon pair emission events is much greater than the characteristic time for each event, so that the limits of the temporal integral may be extended to $\pm \infty$. Note that in writing this expression for $|{\mathrm{\Psi }}_2⟩$, we have defined the field overlap ${\mathrm{\Theta }}_{{\overrightarrow{l}}_1{\overrightarrow{l}}_2{\overrightarrow{l}}_s{\overrightarrow{l}}_i}$, given by $${\mathrm{\Theta }}_{{\overrightarrow{l}}_1{\overrightarrow{l}}_2{\overrightarrow{l}}_s{\overrightarrow{l}}_i}=\int \mathrm{d}\rho \int \mathrm{d}\theta {\rho }^2\sin \theta f_1(\rho ,\theta )f_2(\rho ,\theta )f_s^{*}(\rho ,\theta )f_i^{*}(\rho ,\theta ),$$(7)and the function $G_{{\overrightarrow{l}}_1{\overrightarrow{l}}_2{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_1,{\omega }_s,{\omega }_i)$ expressed in terms of the sphere’s equatorial perimeter $L=2\pi R$, as $$G_{{\overrightarrow{l}}_1{\overrightarrow{l}}_2{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_1,{\omega }_s,{\omega }_i)=\mathrm{sinc}\left[\frac{L\mathrm{\Delta }{k}_{{\overrightarrow{l}}_1{\overrightarrow{l}}_2{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_1,{\omega }_s,{\omega }_i)}{2}\right]\exp \left[i\frac{L\mathrm{\Delta }{k}_{{\overrightarrow{l}}_1{\overrightarrow{l}}_2{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_1,{\omega }_s,{\omega }_i)}{2}\right],$$(8)defined in turn in terms of the phase mismatch function, $$\mathrm{\Delta }{k}_{{\overrightarrow{l}}_1{\overrightarrow{l}}_2{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_1,{\omega }_s,{\omega }_i)=k_{{\overrightarrow{l}}_1}({\omega }_1)+k_{{\overrightarrow{l}}_2}({\omega }_s+{\omega }_i-{\omega }_1)-k_{{\overrightarrow{l}}_s}({\omega }_s)-k_{{\overrightarrow{l}}_i}({\omega }_i).$$(9)

Let us now assume that the two pumps are degenerate (spatially and spectrally), as well as monochromatic at frequency ${\omega }_p$, i.e., $${\alpha }_{{\overrightarrow{l}}_p}(\omega )\equiv {\alpha }_{{\overrightarrow{l}}_1}(\omega )={\alpha }_{{\overrightarrow{l}}_2}(\omega )=\delta (\omega -{\omega }_p).$$(10)

Let us consider a short longitudinal section of the continuous-wave pump of length $\mathrm{\Delta }{z}_p$ (with $\mathrm{\Delta }{z}_p\ll L$). Performing the change of variables $\mathrm{\Omega }={\omega }_s-{\omega }_p$, we then obtain the state resulting from one cavity round trip of this longitudinal pump section of $\mathrm{\Delta }{z}_p$ length $$\begin{gathered} |{\mathrm{\Psi }}_2⟩=\sum _{{\overrightarrow{l}}_p,{\overrightarrow{l}}_s,{\overrightarrow{l}}_i}{A}_{{\overrightarrow{l}}_p}^2{\mathrm{\Theta }}_{{\overrightarrow{l}}_p{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}\int \mathrm{d}\mathrm{\Omega }\ell ({\omega }_p+\mathrm{\Omega })\ell ({\omega }_p-\mathrm{\Omega }) \\ \times g_{{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}(\mathrm{\Omega }){\widehat{b}}_{{\overrightarrow{l}}_s}^{†}({\omega }_p+\mathrm{\Omega }){\widehat{b}}_{{\overrightarrow{l}}_s}^{†}({\omega }_p-\mathrm{\Omega })|0⟩, \end{gathered}$$(11)where the function $g_{{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}(\mathrm{\Omega })$, which constitutes a reduced version of $G_{{\overrightarrow{l}}_1{\overrightarrow{l}}_2{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_1,{\omega }_s,{\omega }_i)$, can be expressed as $$g_{{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_p,\mathrm{\Omega })=\mathrm{sinc}\left[\frac{L\mathrm{\Delta }{\kappa }_{{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_p,\mathrm{\Omega })}{2}\right]\exp \left[i\frac{L\mathrm{\Delta }{\kappa }_{{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_p,\mathrm{\Omega })}{2}\right]$$(12)in terms of the reduced phase mismatch function, $$\mathrm{\Delta }{\kappa }_{{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_p,\mathrm{\Omega })=2k_{{\overrightarrow{l}}_p}({\omega }_p)-k_{{\overrightarrow{l}}_s}({\omega }_p+\mathrm{\Omega })-k_{{\overrightarrow{l}}_i}({\omega }_p-\mathrm{\Omega })-2\gamma P.$$(13)

Note that in the above expression, we have incorporated a nonlinear term in the phase mismatch associated with self- and cross-phase modulation $2\gamma P$ [47]. This term can be written in terms of the $Q$ parameter of the cavity, of the input power $P_{\mathrm{in}}$, of the nonlinear index of refraction $n_2$, and of the effective transverse mode area $A_{\mathrm{eff}}$ as $$2\gamma P=\frac{P_{\mathrm{in}}Q{n}_2}{\pi nR{A}_{\mathrm{eff}}}.$$(14)

Using the fact that the function $\ell (\omega )$ is a slow function of $\omega$, and assuming that each of the three waves (the degenerate pump, signal, and idler) travels in a single spatial mode, the quantum state can be simplified as $|\mathrm{\Psi }⟩=|0⟩+{\eta }^{\prime }|{\mathrm{\Psi }}_2⟩$, where the new constant ${\eta }^{\prime }$ incorporates the quantity $A_{{\overrightarrow{l}}_p}^2{\mathrm{\Theta }}_{{\overrightarrow{l}}_p{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}$, with $|{\mathrm{\Psi }}_2⟩$ given by $$|{\mathrm{\Psi }}_2⟩=\int \mathrm{d}\mathrm{\Omega }{g}_{{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}(\mathrm{\Omega }){\widehat{b}}_{{\overrightarrow{l}}_s}^{†}({\omega }_p+\mathrm{\Omega }){\widehat{b}}_{{\overrightarrow{l}}_s}^{†}({\omega }_p-\mathrm{\Omega })|\mathrm{vac}⟩.$$(15)

So far, the quantum state has been expressed in terms of the annihilation operators ${\widehat{b}}_{{\overrightarrow{l}}_s}$ and ${\widehat{b}}_{{\overrightarrow{l}}_i}$, which correspond to the modes resonating within the sphere. In any realistic experiment, we of course need to couple the photon pairs out of the resonating sphere so as to use them in an experimental setup of interest. This can be accomplished through evanescent coupling of the sphere modes to a given propagation mode of a tapered fiber placed in close proximity to the sphere. Let us denote as ${\widehat{a}}_{\nu }$ the extracavity mode in the tapered fiber (we assume that light couples into a single taper mode) for the signal ($\nu =s$) and idler ($\nu =i$), with $r_{\nu }$ representing the reflectivity of the sphere–taper interface (the probability amplitude corresponding to a photon remaining within the sphere), and ${t^{\prime }}_{\nu }$ representing the transmissivity (the probability amplitude corresponding to a photon coupling from the sphere to the taper); in a lossless cavity, energy conservation dictates ${|{t^{\prime }}_{\nu }|}^2+{|r_{\nu }|}^2=1$, with $|r_{\nu }|=1-\frac{l\pi }{Q}$. After $n+1$ iterations in the cavity the intracavity mode ${\widehat{b}}_{{\overrightarrow{l}}_s}$ and ${\widehat{b}}_{{\overrightarrow{l}}_i}$ can be expressed as follows: $${\hat{b}}_{\overrightarrow{l_v}}^{†}(\omega )\to r_v^{n+1}+e^{{\mathrm{ink}}_{{\overrightarrow{l}}_v}}(\omega )L{\hat{a}}_v^{†}(\omega )+\left[\frac{1-r_v^{n+1}{e}^{i(n+1)k_{\overrightarrow{l_v}(\omega )}L}}{1-r_v{e}^{i{k}_{\overrightarrow{l_v}}(\omega )}L}\right]t_v^{\prime }{\hat{a}}_v^{†}(\omega )\mathrm{.}$$(16)

Taking the limit $n\to \infty$, so that all SFWM light produced by a single cavity round trip of the longitudinal pump section of length $\mathrm{\Delta }{z}_p$ is allowed to escape the cavity, we may write the extracavity mode operator as follows: $$\underset{n\to \infty }{\lim }{\widehat{b}}_{{\overrightarrow{l}}_{\nu }}^{†}(\omega ){|0⟩}_{\nu }=A_{\nu }(\omega ){\widehat{a}}_{\nu }^{†}(\omega ){|0⟩}_{\nu },$$(17)written in terms of the Airy function $A_{\nu }(\omega )$, $$A_{\nu }(\omega )=\frac{{t^{\prime }}_{\nu }}{1-r_{\nu }{e}^{i{k}_{{\overrightarrow{l}}_{\nu }}(\omega )L}}.$$(18)

We may then write the two-photon state propagating in the tapered fiber modes $|{\mathrm{\Psi }}_2^{\prime }⟩$, described by annihilation operators ${\widehat{a}}_s$ and ${\widehat{a}}_i$, as $$\begin{gathered} |{\mathrm{\Psi }}_2^{\prime }⟩=\int \mathrm{d}\mathrm{\Omega }{g}_{{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}(\mathrm{\Omega })A_s({\omega }_p+\mathrm{\Omega })A_i({\omega }_p-\mathrm{\Omega }) \\ \times {\widehat{a}}_s^{†}({\omega }_p+\mathrm{\Omega }){\widehat{a}}_i^{†}({\omega }_p-\mathrm{\Omega })|0⟩. \end{gathered}$$(19)

Let us note that this is the quantum state produced by a single pass of the pump field through the cavity. In an experimental situation of interest, the pump field is resonant in the cavity, so that with a probability amplitude $t$ each pump photon couples from the taper to the sphere, and once within the cavity it remains with a probability amplitude $r_p$ at each pass through the taper–sphere interface. We can then write an expression for the two-photon state $|{\mathrm{\Psi }}_2^{\prime \prime }⟩$, resulting from $n$ iterations of the pump in the cavity, $$\begin{gathered} |{\mathrm{\Psi }}_2^{\prime \prime }⟩=t_p|{\mathrm{\Psi }}_2^{\prime }⟩+t_p{r}_p^2{e}^{i{k}_{{\overrightarrow{l}}_p}({\omega }_p)L}|{\mathrm{\Psi }}_2^{\prime }⟩+t_p{r}_p^4{e}^{i2k_{{\overrightarrow{l}}_p}({\omega }_p)L}|{\mathrm{\Psi }}_2^{\prime }⟩ \\ +\mathrm{...}+t_p{r}_p^{2n}{e}^{i2n{k}_{{\overrightarrow{l}}_p}({\omega }_p)L}|{\mathrm{\Psi }}_2^{\prime }⟩=\frac{t_p[1-r_p^{2(n+1)}{e}^{i{k}_{{\overrightarrow{l}}_p}({\omega }_p)L(n+1)}]}{1-r_p^2{e}^{i{k}_{{\overrightarrow{l}}_p}({\omega }_p)L}}|{\mathrm{\Psi }}_2^{\prime }⟩. \end{gathered}$$(20)

In this expression, upon each successive round trip of the pump longitudinal section $\mathrm{\Delta }{z}_p$ in the cavity, the electric-field amplitudes $A_{{\overrightarrow{l}}_1}$ and $A_{{\overrightarrow{l}}_2}$ are each reduced by $r_p$, so that a factor $r_p^2$ appears in addition to the phase term corresponding to one round trip for each pump (both from a single degenerate pump mode). In the limit $n\to \infty$ corresponding to a situation in which all the pump light from longitudinal section $\mathrm{\Delta }{z}_p$ initially coupled into the cavity has escaped, we may then write the resulting two-photon state $|{\mathrm{\Psi }}_2^{\prime \prime }⟩$ as follows: $$\underset{n\to \infty }{\lim }|{{\mathrm{\Psi }}^{\prime }}^{\prime }⟩=A_p({\omega }_p)|{\mathrm{\Psi }}^{\prime }⟩$$(21)in terms of the Airy function for the pump $A_p(\omega )$, expressed as $$A_p(\omega )=\frac{t_p}{1-r_p^2{e}^{i{k}_{{\overrightarrow{l}}_p}(\omega )L}}.$$(22)

We can then write the two-photon state, which incorporates the full effect of the cavity, as follows (note that the derivation shown here pertains to the case of degenerate pump waves; in the case of nondegenerate pumps, spatially and/or spectrally, two separate Airy functions will appear, one for each of the pumps): $$\begin{gathered} |{\mathrm{\Psi }}^{\prime \prime }({\omega }_p)⟩=|\mathrm{vac}⟩+{\eta }^{\prime }{A}_p({\omega }_p) \\ \times \int \mathrm{d}\mathrm{\Omega }f(\mathrm{\Omega };{\omega }_p){\widehat{a}}_s^{†}({\omega }_p+\mathrm{\Omega }){\widehat{a}}_i^{†}({\omega }_p-\mathrm{\Omega })|\mathrm{vac}⟩ \\ \equiv |\mathrm{vac}⟩+{\eta }^{\prime }{A}_p({\omega }_p)|\phi ({\omega }_p)⟩, \end{gathered}$$(23)in terms of the joint spectral amplitude function $f(\mathrm{\Omega };{\omega }_p)$, $$f(\mathrm{\Omega };{\omega }_p)=g_{{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}(\mathrm{\Omega };{\omega }_p)A_s({\omega }_p+\mathrm{\Omega })A_i({\omega }_p-\mathrm{\Omega }).$$(24)

Alternatively, it is useful for visualization purposes to write the two-photon state in an equivalent form involving the two-dimensional frequency generation space $\{{\omega }_s,{\omega }_i\}$, as follows: $$|{{\mathrm{\Psi }}^{\prime }}^{\prime }({\omega }_p)⟩=|\mathrm{vac}⟩+{\eta }^{\prime }{A}_p({\omega }_p)\int \mathrm{d}{\omega }_s\int \mathrm{d}{\omega }_i{f}_2({\omega }_s,{\omega }_i;{\omega }_p){\widehat{a}}_s^{†}({\omega }_s){\widehat{a}}_i^{†}({\omega }_i)|\mathrm{vac}⟩$$(25)in terms of the two-dimensional joint spectral amplitude, $$f_2({\omega }_s,{\omega }_i;{\omega }_p)=\delta ({\omega }_s+{\omega }_i-2{\omega }_p)G_{{\overrightarrow{l}}_p{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_p,{\omega }_s,{\omega }_i)A_s({\omega }_s)A_i({\omega }_i).$$(26)

It is convenient to re-express this joint amplitude in terms of “rotated” variables, $\mathrm{\Omega }=({\omega }_s-{\omega }_i)/2$ (already introduced) and ${\omega }_{+}=({\omega }_s+{\omega }_i)/2$ as $$f_2^{\prime }({\omega }_{+},\mathrm{\Omega };{\omega }_p)=\delta ({\omega }_{+}-{\omega }_p)\times G_{{\overrightarrow{l}}_p{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_p,{\omega }_{+}+\mathrm{\Omega },{\omega }_{+}-\mathrm{\Omega })A_s({\omega }_{+}+\mathrm{\Omega })A_i({\omega }_{+}-\mathrm{\Omega }).$$(27)

In our experiments, while the pump wave is in the form of a continuous wave with a narrowband linewidth ${\delta }_p$, the pump frequency is varied in time according to a triangular wave with amplitude ${\mathrm{\Delta }}_p$ and frequency $f_p$. As ${\omega }_p$ is swept within the spectral window of width ${\mathrm{\Delta }}_p$, a bandwidth ${\mathrm{\Delta }}_r$ corresponding to the spectral width of the cavity resonance function $|A_p({\omega }_p{)|}^2$ can couple into the cavity. There are therefore three bandwidths of interest that govern the pump and its interaction with the cavity: (i) pump linewidth ${\delta }_p$, (ii) pump resonance bandwidth ${\mathrm{\Delta }}_p$, and (iii) pump sweeping range ${\mathrm{\Delta }}_r$. In our experiments (see below) the relationships ${\delta }_p\ll {\mathrm{\Delta }}_p\ll {\mathrm{\Delta }}_r$ are fulfilled with approximate values ${\delta }_p\approx 200\mathrm{kHz}$, ${\mathrm{\Delta }}_p\approx 20\mathrm{MHz}$, and ${\mathrm{\Delta }}_r\approx 25\mathrm{GHz}$.

Let us now express the two-photon state produced by the sphere as a statistical mixture of all the pure states produced by each individual pump spectral component ${\omega }_p$ that can couple into the sphere, as follows: $$\widehat{\rho }=\int \mathrm{d}{\omega }_p{|A_p({\omega }_p)|}^2|\phi ({\omega }_p)⟩⟨\phi ({\omega }_p)⟩|.$$(28)

We can now write down an expression for the spectral intensity (SI) for the idler photon $R_i(\mathrm{\Omega })$, in terms of the photon number operators $\widehat{n}({\omega }_{\nu })=a^{†}({\omega }_{\nu })a({\omega }_{\nu })$ (with $\nu =s,i$), as follows: $$\begin{gathered} R_i(\mathrm{\Omega })=\int \mathrm{d}{\omega }_s⟨\widehat{n}({\omega }_s)\widehat{n}({\omega }_p-\mathrm{\Omega })⟩ \\ =\mathrm{Tr}[{\widehat{a}}^{†}({\omega }_s){\widehat{a}}^{†}({\omega }_p-\mathrm{\Omega })\widehat{a}({\omega }_p-\mathrm{\Omega })\widehat{a}({\omega }_s)\widehat{\rho }] \end{gathered}$$(29)$$=\int \mathrm{d}{\omega }_p{|A({\omega }_p)f(\mathrm{\Omega };{\omega }_p)|}^2.$$(30)

An explicit version of the above equation is as follows, where in the second line we have used the approximation that the state is defined by the three Airy functions with a negligible effect of the phase-matching function: $$\begin{gathered} R_i(\mathrm{\Omega })=\int \mathrm{d}{\omega }_p{|g_{{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}(\mathrm{\Omega })|}^2{\mathcal{A}}_p({\omega }_p){\mathcal{A}}_s({\omega }_p+\mathrm{\Omega }){\mathcal{A}}_i({\omega }_p-\mathrm{\Omega }) \\ \approx \int \mathrm{d}{\omega }_p{\mathcal{A}}_p({\omega }_p){\mathcal{A}}_s({\omega }_p+\mathrm{\Omega }){\mathcal{A}}_i({\omega }_p-\mathrm{\Omega }) \end{gathered}$$(31)in terms of the intensity Airy functions ${\mathcal{A}}_{\nu }(\omega )\equiv |A_{\nu }(\omega {)|}^2$ (with $\nu =p,s,i)$. It is of interest to write equivalent expressions in time domain, for which we define time-domain annihilation operators as follows: $${\widehat{a}}_t(t)=\frac{1}{2\pi }\int \mathrm{d}\omega \widehat{a}(\omega )e^{-i\omega t}.$$(32)

We can then write down an expression for the resulting two-photon times of emission distribution (TEDs) $\tilde{R}(t_s,t_i)$ as $$\begin{gathered} \tilde{R}(t_s,t_i)=⟨{\widehat{n}}_t(t_s)n_t(t_i)⟩=\mathrm{Tr}[{\widehat{a}}_t^{†}(t_s){\widehat{a}}_t^{†}(t_i){\widehat{a}}_t(t_i){\widehat{a}}_t(t_s)\widehat{\rho }] \\ =\int \mathrm{d}{\omega }_p{\mathcal{A}}_p({\omega }_p)\int \mathrm{d}\mathrm{\Omega }{|f(\mathrm{\Omega };{\omega }_p)|}^2{e}^{-i(t_s-t_i)\mathrm{\Omega }} \\ =\int \mathrm{d}\mathrm{\Omega }\int \mathrm{d}{\omega }_p{\mathcal{A}}_p({\omega }_p){|f(\mathrm{\Omega };{\omega }_p)|}^2{e}^{-i(t_s-t_i)\mathrm{\Omega }} \\ =\int \mathrm{d}\mathrm{\Omega }{R}_i(\mathrm{\Omega })e^{-i(t_s-t_i)\mathrm{\Omega }}. \end{gathered}$$(33)

Note that in the third line, we have interchanged the order of integration, leading to a Fourier transform relationship between the SI and the TED (fourth line).

With the help of the above theory, we can now describe the spectral and temporal properties of the two-photon state produced by SFWM in a specific source design. In this section we present simulations of the photon-pair properties of interest, assuming experimental parameters consistent with our experiments described in Section 4. Let us consider a fused silica sphere of radius $R=135\mathrm{\mu m}$. For this radius, the degenerate SFWM frequency ${\omega }_s={\omega }_i={\omega }_p=2\pi c/1550.92\mathrm{nm}$ corresponds most closely with azimuthal index values $l_s=l_i=774$ [i.e., this value of $l$ comes closest to fulfilling the resonance condition $k({\omega }_l)L=2\pi l$].

In Figs. 2(a) and 2(b) we plot for two values of the $Q$ parameter (namely $Q{=10}^6$ and $Q{=10}^8$), for a radius of 135 μm, the phasematching strength $|g_{{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_p,\mathrm{\Omega }{)|}^2$, as a function of the pump frequency ${\omega }_p$ in the horizontal axis and the signal-photon generation frequency (detuned from the pump) $\mathrm{\Omega }$ in the vertical axis. Note that fixing the input pump power to $P_{\mathrm{in}}=7\mathrm{mW}$, the $Q$ value defines $2\gamma P$ through Eq. (14). In each of the plots we also show the phasematching contours $L\mathrm{\Delta }k=0$ and $L\mathrm{\Delta }k=0.49$, the latter value chosen because it yields $\mathrm{sinc}(L\mathrm{\Delta }k/2)\approx 0.99$. In addition, a rectangle appearing close to the center of each plot represents the pump/generation spectral region of interest in our experiments. From these plots we can draw the following two conclusions: (i) the phasematching term $g_{{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_p,\mathrm{\Omega })$ is essentially constant within the spectral area of interest; (ii) changes in the nonlinear term $2\gamma P$ (within the experimental range of interest) do not affect the resulting photon-pair properties. What this means is that we are able to approximate $g_{{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_p,\mathrm{\Omega })\approx 1$, with the implication that the two-photon state will be entirely determined by the cavity properties through the Airy functions ${\mathcal{A}}_p(\omega )$, ${\mathcal{A}}_s(\omega )$, and ${\mathcal{A}}_i(\omega )$, as in the second line of Eq. (31).

In this context, so as to illustrate the effect of using extended rather than micro-cavities, in Figs. 2(c) and 2(d) we show, for cavity radii of 1.35 cm and 1.35 m and a $Q$ value of $Q={10}^8$ the phasematching strength $|g_{{\overrightarrow{l}}_p{\overrightarrow{l}}_s{\overrightarrow{l}}_i}({\omega }_p,\mathrm{\Omega }{)|}^2$, as a function of ${\omega }_p$ and $\mathrm{\Omega }$. As in panels Figs. 2(a) and 2(b) we show with a rectangle the spectral area of interest in our experiments. Note that in contrast with microcavities which are the focus of this work, for extended cavities phasematching does become a relevant factor in defining the two-photon state.

Considering the discussion above, it is the cavity resonances for the signal and idler modes, as well as for the pump, that determine the two-photon state. Note that along the signal frequency ${\omega }_s$ axis, the two-photon amplitude can be non-zero only within each of the resonances (each corresponding to a particular value of $l_s$), and likewise for the idler frequency ${\omega }_i$ axis. Therefore, the two-photon state can be non-zero in $\{{\omega }_s,{\omega }_i\}$ space only within each particular mode belonging to a matrix of modes, defined by all combinations of $l_s$ and $l_i$, centered around $l_s=l_i=774$.

We point out that the spectral separation between two neighboring spectral modes, also known as FSR, has a slight dependence on frequency, as plotted in Fig. 3(a) from Eq. (4). The effect of such a spectral drift in the FSR on the two-photon state structure is illustrated in Figs. 3(b)–3(e). In Fig. 3(b) we present an illustration in $\{{\omega }_s,{\omega }_i\}$ space of the matrix of generation modes under the assumption of a constant FSR. We also indicate in this plot the rotated axes ${\omega }_{+}$ and $\mathrm{\Omega }$, where a particular monochromatic pump corresponds to fixing the value of ${\omega }_{+}$ as ${\omega }_{+}={\omega }_p$. An appropriate choice of ${\omega }_p$ so as to ensure overlap with the vertices of the resulting squares leads to a state in the form of a frequency comb along the main diagonal of the generation-mode matrix, as indicated in the figure. In Fig. 3(c) we show the same information in a rotated frequency space $\{{\omega }_{+},\mathrm{\Omega }\}$, so that the generation modes appear along a horizontal line. Figures 3(d) and 3(e) are analogous to Figs. 3(b) and 3(c), except that we have included the effect of the spectral drift of the FSR. It can be appreciated that the two-photon state structure now increasingly departs from the diagonal for large $|\mathrm{\Omega }|$ (with a concave locus of all possible generation modes).

In our two-photon state modeling, a prominent role is played by the resonance function for the pump ${\mathcal{A}}_p(\omega )$. This function may be characterized experimentally; see Section 4.A. While in Fig. 3 we present sketches designed to explain the geometry of the resulting SFWM two-photon state, in Fig. 4 we present an actual plot of the joint spectral intensity (SI) $|f_2({\omega }_s,{\omega }_i;{\omega }_p{)|}^2$ [see Eq. (25)] for the specific case of a sphere of $R=135\mathrm{\mu m}$ radius, assuming $Q$ values for the signal and idler modes of $Q_s=Q_i=1\times {10}^8$, and a full width at half-maximum width for function ${\mathcal{A}}_p(\omega )$ of 20.4 MHz. Note that we have plotted the joint SI at each location on the main diagonal of the generation-mode matrix, within a square of 300 MHz side, and a center-to-center square separation of 1.52 THz. In this plot, we have included along each of the ${\omega }_s$, ${\omega }_i$, and $\mathrm{\Omega }$ (diagonal) axes the resulting marginal intensity distributions. In the inset, we show the central square corresponding to $l_s=l_i=774$, where we also include plots of ${\mathcal{A}}_s(\omega )$, ${\mathcal{A}}_i(\omega )$, and ${\mathcal{A}}_p(\omega )$, along each of the ${\omega }_s$, ${\omega }_i$, and ${\omega }_{+}$ axes.

In order to further clarify the two-photon state structure, we plot in Fig. 5 within similar square regions as in Fig. 4 (this time defined in the rotated variables $\{\mathrm{\Omega },{\omega }_{+}\}$), the function ${\mathcal{A}}_s({\omega }_{+}+\mathrm{\Omega })$ in Fig. 5(a), the function ${\mathcal{A}}_i({\omega }_{+}-\mathrm{\Omega })$ in Fig. 5(b), the product ${\mathcal{A}}_s({\omega }_{+}+\mathrm{\Omega }){\mathcal{A}}_i({\omega }_{+}-\mathrm{\Omega })$ in Fig. 5(c), and the product ${\mathcal{A}}_s({\omega }_{+}+\mathrm{\Omega }){\mathcal{A}}_i({\omega }_{+}-\mathrm{\Omega }){\mathcal{A}}_p({\omega }_{+})$ in Fig. 5(d). Here, we may appreciate the effect already discussed through the sketches in Fig. 3, in which the generation modes increasingly depart from the axis ${\omega }_{+}={\omega }_p$ for larger values of $|\mathrm{\Omega }|$. In Fig. 5(c), we also indicate with two dotted lines the width of the pump resonance function ${\mathcal{A}}_p({\omega }_{+})$. The result is that the generation modes lie increasingly outside of the pump resonance and are therefore suppressed for large $|\mathrm{\Omega }|$, as is clear in Fig. 5(d). In Fig. 5(e), we show the resulting photon pair SI in the form of a frequency comb, as is obtained from Eq. (31).

In Fig. 6, we show the idler–photon SI, similar to that shown in Fig. 5(e), for three different signal/idler $Q$ values ($1\times {10}^6$, $1\times {10}^7$, and $1\times {10}^8$), where we also indicate the corresponding cavity reflectivity coefficients. Note that for a comparatively smaller $Q$ parameter, resulting in spectrally broader generation modes, the envelope decays more slowly because spectral overlap is retained with the pump resonance function ${\mathcal{A}}_p(\omega )$ over a larger span of $\mathrm{\Omega }$ values.

It is of interest to provide a temporal description of the two-photon state in addition to the spectral description already provided. In Fig. 7 we show, as calculated numerically from Eq. (33) and from the traces in Fig. 6, the TED for the same three $Q$ values [Figs. 7(a)–7(c)]. Note that each of these three functions exhibits an envelope with closely spaced oscillations (with a period corresponding to the cavity round-trip time), which cannot be resolved in those plots. In Fig. 7(d), we show a closeup of the curve shown in Fig. 7(a), within the region of the maximum showing these temporal oscillations.

As can be appreciated from Fig. 6, the resulting single-photon SI is in the form of a frequency comb, with the relative heights of the peaks modulated by an envelope function. While the FSR exhibits a slow frequency dependence [as indeed has been shown in Fig. 6(a)], within a restricted spectral window we may model the SI function as a fixed FSR comb function as follows: $$R_i(\mathrm{\Omega })=H(\mathrm{\Omega })·h(\mathrm{\Omega })*{\mathrm{comb}}_{\delta \mathrm{\Omega }}(\mathrm{\Omega }),$$(34)where $H(\mathrm{\Omega })$ represents the envelope function, $h(\mathrm{\Omega })$ represents an individual peak in the comb, and the symbol $*$ denotes a convolution. In addition, this equation is written in terms of the Dirac–delta comb function ${\mathrm{comb}}_{\mathrm{\Delta }}(x)$, defined as follows: $${\mathrm{comb}}_{\mathrm{\Delta }}(x)=\sum _{j=-\infty }^{\infty }\delta (x-j\mathrm{\Delta }).$$(35)

Let functions $\tilde{H}(T)$ and $\tilde{h}(T)$ represent the Fourier transforms of $H(\omega )$ and $h(\omega )$, respectively. Under the assumption that the function $\tilde{H}(T)$ is much narrower than function $\tilde{h}(T)$, we can then write the joint temporal intensity (or temporal emission distribution, TED, function) $\tilde{R}(T)$, with $T=t_s-t_i$, as follows: $$\tilde{R}(T)=\tilde{h}(T)·\tilde{H}(T)*{\mathrm{comb}}_{1/\delta \mathrm{\Omega }}(T).$$(36)

Thus, the joint temporal amplitude function is composed of a comb function in the temporal variable $T=t_s-t_i$ with individual comb peaks defined by the function $\tilde{H}(T)$, with a peak-to-peak separation of $1/\delta \mathrm{\Omega }$, and an envelope function $\tilde{h}(T)$ that describes the roll-off in amplitude for large $|T|$. Interestingly, the roles of the function $H(\mathrm{\Omega })$ and $h(\mathrm{\Omega })$ in the frequency domain and $\tilde{H}(T)$ and $\tilde{h}(T)$ in the temporal domain are reversed. Specifically, the spectral envelope defines the functional dependence of each individual temporal comb peak, and the functional dependence of each individual frequency comb peak defines the envelope of the resulting comb in the temporal domain.