Direct generation of entangled photon pairs in nonlinear optical waveguides

2.1 Down-conversion efficiency in cylindrical waveguides

To quantify the SPDC efficiency, we compute the probability of the inverse process: SF generation produced by two counter-propagating guided photons of frequencies ω _i and ω _{i ^′} which combine to generate a photon frequency ω _{ii ^′} = ω _i + ω _{i ^′} that is normally emitted from the waveguide. In virtue of reciprocity, the per-photon probabilities for the two processes (SPDC and SF generation) are identical. In practice, we calculate the efficiency by considering two photons within counter-propagating guided modes i and i′, prepared as long pulses of length L and space/time-dependent electric fields E _i (r, t) and E _{i ^′} (r, t) (Figure 1(a)) that comprise frequency components that are tightly packed around ω _i and ω _{i ^′} . Through the SF second-order susceptibility tensor χ⁽²⁾, a polarization density P _{ii ^′} (R) is produced within a narrow frequency range around ω _{ii ^′} . More precisely,

where the indices {a, b, c} run over Cartesian components, E _i (R) gives the profile of mode i in the transverse plane R = (x, y), and we normalize the susceptibility to the quantity

For simplicity, we consider the wave vectors of the two modes to satisfy the condition q _i + q _{i ^′} = 0, so that the SF photons are emitted with zero wave vector component parallel to the waveguide (i.e., along normal directions). The SF polarization density generates a field that we compute at long distances from the waveguide using the electromagnetic Green tensor of the system G(r,r′,ω) [22], from which we calculate the far-field Poynting vector, whose radial component is in turn integrated over time and directions of emission to produce the emitted energy. We then divide this energy by ℏω _{ii ^′} to obtain the number of emitted photons N _{ii ^′} . Likewise, we calculate the Poynting vector associated with each of the pulses and integrate the component parallel to the waveguide over time and transverse spatial directions to yield the number of photons incident in each pulse, N _i and N _{i ^′} . Finally, the ratio of the emitted number of photons to the number of photons in each pulse is interpreted as the probability η _{ii ^′} = N _{ii ^′} /N _i N _{i ^′} that two colliding quanta combine into one emitted SF quantum (again, identical with the probability that an externally incident photon produces a pair of counter-propagating photons within modes i and i′). In the long L limit, the incident pulses become monochromatic and η _{ii ^′} turns out to be independent of L. For convenience, we decompose the up-conversion efficienty into contributions associated with emission along different azimuthal angles φ (see coordinate system in Figure 1(a)) as

After a lengthy calculation (see a detailed self-contained derivation in Methods), we obtain the following result for the angle-resolved efficiency:

where w _i = ω _i a/c, β _i = v _i /c, v _i = ∂ω _i /∂q _i is the group velocity in mode i, H is the magnetic field, and g(φ − φ′, R′, ω) is the amplitude of the electromagnetic Green tensor in the far-field limit defined through G(r,r′,ω)→(eiϵhωR/c/R)g(φ−φ′,R′,ω) for normal emission (see Eq. (34) below for an explicit expression). Here, I _i and I _{ii ^′} are proportional to the number of incident photons in waveguide mode i and emitted outside the waveguide, respectively. Incidentally, these coefficients are normalized in such a way that they are independent of the waveguide radius a, so the efficiency η _{ii ^′} only depends on a through an overall factor 1/a⁴ for a fixed value of ω _i a/c. In brief, η _{ii ^′} represents the ratio of SF photons produced per two incident photons (one in each waveguide mode), that is, the SF matrix element for ω _i + ω _{i ^′} → ω _{ii ^′} , which must be equal to the SPDC matrix element corresponding to ω _{ii ^′} → ω _i + ω _{i ^′} . The latter affects each incident photon separately, so it can be interpreted as the fraction of incident photons that undergo SPDC, and we thus equate it to the fraction of down-converted power.

For the cylindrical waveguides under consideration, we can multiplex the mode labels as i = {q _i , m _i , l _i , σ _i }, where q _i is the wave vector, m _i is the azimuthal angular momentum number, l _i refers to different radial resonances, and σ _i runs over polarization states (i.e., TE0li and TM0li for m _i = 0, and hybrid modes HEmili and EHmili for m _i ≠ 0, see Section 4.1.1). Given the symmetry of the waveguide, the radial and azimuthal components of the transverse field associated with each mode only depend on radial distance R, apart from an overall phase factor eimiφ. For simplicity, we consider a second-order response tensor χ⁽²⁾ that also preserves the cylindrical symmetry, so that the angular integral in Eq. (5c) leads to angular momentum conservation (m _{ii ^′} = m _i + m _{i ^′} for the emitted photons). In particular, we assume a nonlinear tensor dominated by the χzzz(2) component (e.g., a LiNbO₃ waveguide with the z axis aligned along the waveguide), which implies that the TE modes and the TE component of the hybrid modes do not couple to the incident field through χ⁽²⁾. In fact, as we show in the Supplementary Information, the mode profiles and dispersion properties of cylindrical waveguides are remarkably similar to their counterparts with rectangular geometries, such as those more commonly employed in the construction of waveguides from highly nonlinear crystals in which a single component of the χ⁽²⁾ tensor dominates [37, 38].

2.2 Availability and efficiency of different down-conversion channels

We are now equipped to discuss the generation of entangled photon pairs through SPDC in our waveguide. Assuming the above conditions, the lowest-frequency modes that possess a nonzero z component of the electric field, and can consequently couple to normally impinging external light, are HE₁₁ and TM₀₁ (see Figure 1(b)). We identify two relevant frequencies in this region (see Figure 1(c)): the threshold of the TM₀₁ mode at ω₁ (satisfying ω1a/c=2.4048/ϵ1−ϵh, see Section 4.1), and the frequency ω₂ of mode HE₁₁ with the same wave vector. Upon inspection, we find that for an incident light frequency ω < ω₁ + ω₂, the only SPDC channel available corresponds to the generation of two HE₁₁ modes of frequency ω/2 and opposite wave vectors. This situation already allows us to produce entangled photon pairs of the form given in Eq. (1), where i and i′ now refer to the azimuthal numbers m _i , m _{i ^′} ∈ {−1, 1} for each of the emitted photons. In particular, if the waveguide is symmetrically illuminated along different azimuthal directions (see Section 2.3 below), it is possible to select only the m = 0 component from the external light, so that conservation of azimuthal angular momentum leads to the condition m _i + m _{i ^′} = 0, and therefore, the emitted photon pair forms an entangled state |L₋₁R₁⟩ + |L₁R₋₁⟩, where the subindices indicate the values of m _i and m _{i ^′} for the L and R emission directions, all of them sharing the same frequency ω/2 and polarizations HE_±1,1, so we refer to this channel as HE₁₁ + HE₁₁.

Another interesting range of incidence frequencies is ω₁ + ω₂ < ω < 2ω₁ (blue area in Figure 1(c)), where the HE₁₁ + HE₁₁ channel is now supplemented by two additional possibilities in which the two generated photons have different frequencies (with the sum satisfying ω = ω _i + ω _{i ^′} ) and lie in different bands (HE₁₁ or TM₀₁). This is indicated by the two pairs of color-matched blue and green dots in Figure 1(b) and (c), where the condition of opposite wave vectors is obviously satisfied. Again, it is possible to select a specific SPDC channel by illuminating with a fixed m number (see below), and in particular, by setting m = m _i + m _{i ^′} = 1, the HE₁₁ + HE₁₁ channel is eliminated (because the overall azimuthal number obtained by combining two HE_±11 modes is 0 or ±2), so that we obtain again a maximally entangled state of the form given in Eq. (1) with i and i′ now referring to TM₀₁ and HE₁₁ (i.e., |L_TMR_HE⟩ + |L_HER_TM⟩, with azimuthal numbers m _i and m _{i ^′} taking the values 0 and 1 in the TM and HE components, respectively).

The formalism presented in Section 2.1 allows us to calculate the SPDC efficiency for the production of specific photon pair-states, using external illumination prepared with an azimuthal number m = m _i + m _{i ^′} and a polarization state determined by the time reversal of the SF generation state considered in the derivation of these results. Under the assumed conditions of incidence along transverse directions, and considering a zzz-dominant component in the second-order susceptibility tensor, the profile of the applied light amplitude as a function of azimuthal angle φ is therefore taken to be e^imφ, with the field oriented parallel to the waveguide direction. These conditions can be met by combining several incident light beams, as we discuss below. The efficiencies calculated for different SPDC channels in this scheme are shown in Figure 2, normalized to (ℏc/a4)|χ̄(2)|2 in order to present universal, dimensionless results as a function of the scaled incident light frequency ωa/c. For an ϵ₁ = 5 waveguide in air (Figure 2(a)), we find efficiencies that generally grow with the order of the waveguide modes, exhibiting resonances as a function of the incident frequency ω = ω _i + ω _{i ^′} . These resonances are inherited from the two-dimensional transmission coefficients (see Section 4.1) and can be understood as coupling of the incident light to leaky cavity modes at the incident light frequency. We indicate in white the area in which there is only one decay channel (HE₁₁ + HE₁₁, see above), whereas the area with three decay channels (two additional ones corresponding to TM₀₁ + HE₁₁ and HE₁₁ + TM₀₁) is highlighted in blue. The green region at higher frequencies contains an increasing number of channels, which could be also exploited to generate more complex entangled mixtures, involving multiple output states in each direction (L and R) and higher-order modes.

Figure 2:

Down-conversion efficiency for different output channels. (a) Normalized SPDC efficiency for a waveguide with ϵ₁ = 5 and ϵ_h = 1 as a function of incident light frequency ω = ω _i + ω _{i ^′} . Different output channels i + i′ are indicated by labels, while the regions highlighted in white, blue, and green indicate 1, 3, and >3 available channels (each of them degenerate in the sign of both the wave vectors and the azimuthal numbers). (b) Efficiencies corresponding to the three lowest channels (one in HE₁₁ + HE₁₁ and two in HE₁₁ + TM₀₁) for different values of the waveguide permittivity ϵ₁ (see color-coded labels). Colored areas highlight the respective regions in which three channels exist, while the frequency threshold for each of the down-conversion channels is indicated by colored circles.

The spectral evolution of the efficiencies is roughly maintained when varying the waveguide permittivity ϵ₁ (Figure 2(b)), but we observe a general increase in η _{ii ^′} with increasing ϵ₁ in the region of interest, as well as a spectral shift of the region with three output channels (highlighted in shading colors and evolving toward lower frequencies as we increase the permittivity, in agreement with the single-mode-fiber cutoff condition). Interestingly, we find a crossover in the efficiency of HE₁₁ + HE₁₁ relative to that of HE₁₁ + TM₀₁: the former dominates over the latter within the three-channel region at high ϵ₁, whereas the opposite behavior is found at lower permittivities.

Quantitatively, our results indicate that the current scheme is feasible for producing a reasonable rate of entangled photon pairs, taking into account that they are already prepared within waveguide modes [31, 33]. In particular, for values of |χ̄(2)|∼10−10 m V⁻¹ found in good nonlinear materials such as LiNbO₃ [17, 39] and a waveguide radius ∼100 nm, the scaling factor in Figure 2 is (ℏc/a4)|χ̄(2)|2∼10−10 upon conversion of the arguments to the Gaussian-GCS unit system adopted here, which yields a power fraction of 10⁻¹¹ for SPDC when it is multiplied by a scaled efficiency of ∼0.1 (Figure 2). Considering photon energies ∼1 eV and an incident light power ∼1 mW, the estimated power fraction amounts to a generation rate of ∼10⁵ entangled photon pairs per second. As an additional possibility, the efficiency could be increased by incorporating resonant elements to amplify the external light in the region surrounding the waveguide, such as planar Fabry–Perot resonators, which is a natural option for waveguides fabricated on a substrate [34, 37]. Alternatively, the SPDC efficiency can be improved by physically or chemically doping the waveguide with active centers to enhance its intrinsic nonlinear response, or by incorporating different materials in lithographically patterned waveguide geometries.

2.3 Selection of down-conversion channels through illumination interference

A p-polarized electromagnetic plane wave of amplitude E₀ incident through the host medium with a wave vector kh⊥ẑ normal to the waveguide (Figure 3(a)) contributes with a broad range of azimuthal numbers m = m _i + m _{i ^′} according to the decomposition

in terms of cylindrical waves Eh,0mpJ (see Supplementary Information). This situation leads to entangled states that combine more than two polarizations for each of the two waveguiding directions (L and R). For example, in the single HE₁₁ + HE₁₁ channel region (at incident light frequency ω < ω₁ + ω₂), we can have all combinations of m _i = ±1 and m _{i ^′} = ±1, thus reducing the degree of entanglement. One way to address this issue is to combine illumination from different azimuthal directions (e.g., in an interferometric setup involving beamsplitters and mirrors, as illustrated in Figure 3(b)). In particular, when illuminating with two in-phase counter-propagating waves (Figure 3(c)), only even values of m survive, whereas only odd m’s are selected if the waves have a π relative phase difference (Figure 3(d)). In general, we can consider an arbitrary number of plane waves (Figure 3(e)), so that the total field acting on the waveguide has the same form as in Eq. (6), but with E0e−imφkh substituted by ∑jEje−imφj, where the sum runs over plane waves j of amplitude E _j directed along azimuthal directions φ _j .

Figure 3:

Mode selection through light interference. A single incident light plane wave (a) contains all possible values of the azimuthal number m in the external field, and can thus excite all available SPDC channels for the chosen input frequency. We can combine illumination from different directions through beam splitters and mirrors, as indicated in (b) for two-plane-wave irradiation, leading to a selection of the m values from the external light. Examples of selection by irradiation with two in-phase and out-of-phase counter-propagating plane waves are shown in (c) and (d). More stringent selection of m is possible by combining multiple plane waves of amplitudes E _j along different azimuthal directions φ _j , with j = 0, …, n.

In the one-channel regime (HE₁₁ + HE₁₁ output), we can generate the maximally entangled state |L₋₁R₁⟩ + |L₁R₋₁⟩ by selecting an incident m = 0 component and eliminating m = ±2 contributions, as other values of m do not couple to the output modes that are available in that region. This selection requires a minimum of three external plane waves (e.g., with equal amplitudes and azimuthal angles of 0 and ±π/3). Likewise, a more stringent selection of m contributions is possible by resorting to more incident plane waves, therefore opening a vast range of possible entangled photon pairs prepared in higher-order modes. In particular, a simple discrete-transform analysis leads to the conclusion that we can select a specific m = m₀ component while canceling out those of all other m’s in the |m| ≤ N range by combining 2N + 1 incident waves along azimuthal directions φ _j = 2πj/(2N + 1) with amplitudes Ej=E0eim0φj for 0 ≤ j ≤ 2N. For example, in the three-channel region of Figure 2, applying this procedure with N = 2, five plane waves can be used to select m = 1 (or m = −1) while discarding the undesired m = 0, ±2 and m = −1 (or m = 1) contributions to obtain a maximally entangled state |L_TMR_HE⟩ + |L_HER_TM⟩ that combines the TM₀₁ and HE₁₁ (or HE₋₁₁) modes.

4.1 Electromagnetic waves in a cylindrical waveguide

To describe electromagnetic waves in a cylindrical geometry, we first decompose the electric field into cylindrical waves following the prescription of Ref. [40]. More specifically, adopting a cylindrical coordinate system r = (R, φ, z), we consider a homogeneous, isotropic dielectric medium (labeled j) free of external charges and currents that is characterized by a permittivity ϵ _j (setting the magnetic permeability to μ = 1) and express the electric field in cylindrical waves indexed by their azimuthal number m, wave vector q along ẑ, and polarization σ ∈ {s, p} according to

where we define kj=ϵjω/c and Qj=kj2−q2+i0+ (with the square root yielding a positive real part), while the primes on the Bessel functions denote differentiation with respect to the argument. From the orthogonality of the Bessel functions ∫0∞xdxJm(x)Jm(ax)=δ(a−1), it is easy to show that these fields satisfy the orthonormality relation ∫d2REj,qmσJ⋅Ej,q′m′σ′J*=2πδmm′δσσ′δ(q−q′)/q, while the field of modes with different polarizations are related as

We now discuss a cylindrical wave emanating from the interior of a cylindrical waveguide of radius a that is infinitely extended in the z direction, comprised of a dielectric material of permittivity ϵ₁ (medium j = 1), and embedded in a host medium j = h (permittivity ϵ_h). Using the notation introduced above, the electric field is expressed as

where

are outgoing waves similar to the propagating waves in Eq. (7), but with the Bessel functions J _m substituted by Hankel functions Hm(1), while the reflection and transmission coefficients r_m,σσ^′ and t_m,σσ^′ are given by (see Supplementary Information)

with the matrix M defined as [40]

with k = ω/c. The above result is equivalent to other textbook forms of the dispersion relation for cylindrical waveguide modes [41–44], typically expressed in terms of modified Bessel functions in lieu of Hankel functions. The range of wavelengths λ for which only a single mode exists is determined by the condition (a/λ)ϵ1−ϵh<(α0/2π)=0.3827, where α₀ is the first zero of J₀, [45] thus setting a wavelength threshold for the multimode waveguide here considered.

4.2 Field produced by an inner line dipole in the region outside the waveguide

We consider a line dipole placed at a transverse position R₀ = (x₀, y₀) within the waveguide and represent it by a dipole density peiqz0 (dipole per unit length) extending along the line defined by varying z₀ in (R₀, z₀). It is useful to begin by calculating the electric field produced in a homogeneous medium with the same permittivity ϵ₁ as the waveguide material, expressed as the integral over z₀ of the field due to a point dipole [30]:

where ∇ is understood to act on r, whereas the integral can be evaluated using the identity ∫dz0exp(ik1|r−r0|+iqz0)/|r−r0|=iπeiqzH0(1)Q1|R−R0| with Q₁ defined as in Eq. (7). The line dipole should generate a set of outgoing cylindrical waves (therefore the Hankel functions) centered at R₀, so in order to capitalize on the axial symmetry of the waveguide, we need to express the field in terms of waves centered at the origin R = 0. To this end, we invoke Graf’s theorem (see Eq. 9.1.79 in Ref. [47]), H0(1)Q1|R−R0|=∑mHm(1)(Q1R)Jm(Q1R0)eim(φ−φ0), which holds for |R| > |R₀| and can thus be used to describe the dipole field in the waveguide surface region R = a > R₀, through which the outgoing waves are partially transmitted outside the waveguide. This translation formula allows us to recast Eq. (17) into

Now, projecting the dipole as p=∑±p±(x̂±ŷ)/2+pzẑ, where

where the fields E1,qmσH(R) are defined in Eq. (7).

4.3 Sum-frequency generation by counter-propagating waveguided pulses

We now introduce counter-propagating pulse fields E _i (r, t) and E _{i ^′} (r, t) of the form given in Eq. (15), oscillating at frequencies ω _i and ω _{i ^′} , respectively. Through the second-order nonlinearity of the waveguide material χabc(2), where the subscripts a, b, and c denote Cartesian components, a polarization density P _{ii ^′} is produced at frequency ω _{ii ^′} = ω _i + ω _{i ^′} . More precisely,

where |χ̄(2)|, defined in Eq. (3), is introduced to quantify the strength of the SF susceptibility, while the normalized polarization density can be separated as

by defining P̃ii′(R) as in Eq. (2), as well as the (z, t)-dependent factor

Eventually, we set q _i = −q _{i ^′} (i.e., waveguided photons with opposite wave vectors, leading to normal emission of SF photons), so v _i and v _{i ^′} also have opposite signs inherited from q _i and q _{i ^′} . In addition, the nonlinear susceptibility χabc(2)(r) is taken to be constant over the range of frequencies under consideration, as well as uniform inside the waveguide and zero outside it. We note that the present analysis could be trivially extended to consider a nonlinear cladding instead.

The SF field E _{i i ^′} (r, t) is thus produced by the nonlinear polarization density according to

where we have introduced the three-dimensional electromagnetic Green tensor G(r,r′,ω), defined in such a way that E(r,r′)=G(r,r′,ω)⋅p gives the field produced at r by a point dipole of strength p placed at r′. To quantify SF generation in the far-field limit (i.e., at k_hr ≫ 1), we exploit the translational invariance of the polarization source to write

which allows us to recast Eq. (26) as

in terms of the frequency-space electric far-field amplitude

Finally, as shown in Supplementary Information, we can relate the far-field three-dimensional amplitude obtained from the Green tensor in Eqs. (27) and (29) to the two-dimensional one presented in Section (4.2.1) as

with j taking the values +, −, and z. Again, explicit expressions for the components of S(R̂,R′,q,ω) are offered in Section 4.2.1.

4.4 Up- and down-conversion efficiency

The efficiency of the SPDC process in which a photon impinging normally to the waveguide direction produces a pair of guided photons moving in opposite directions away from one another is argued to be identical with that of an up-conversion process involving SF generation of two guided photons moving toward one another, provided that reciprocity applies. We then calculate the SF efficiency η _{ii ^′} = N _{ii ^′} /N _i N _{i ^′} as the ratio of the number of emitted SF photons N_ii′ to the number of incident photons N _i and N_i′ in both guided pulses.

The number of photons carried by a guided mode pulse with field profile E _i (r, t) is expressed as

where we evaluate the flux carried by the Poynting vector S _i (r, t) = (c/4π)E _i (r, t) ×H _i (r, t) in the R plane. Inserting the fields given by Eqs. (15) and (16) and integrating over time, we obtain

Likewise, the number of photons produced in the far field via SF generation from the radial component of the energy emanating from the waveguide can be in turn obtained from the far-field Poynting vector Sii′∞ (evaluated from the field in Eq. (28)) as

In the derivation of the above expression, we have used the fact that fii′(r̂,ω)⋅r̂=0 (i.e., the far field is transverse). For long pulses (see Eq. (15)), only the first term in Eq. (24) (peaked around frequencies ω ∼ ω_ii′) contributes to fii′(r̂,ω) over the ω > 0 integral in Eq. (31), and therefore, using Eqs. (23), (24), and (29), we can write