Pure down-conversion photons through sub-coherence length domain engineering

TL;DR

This paper introduces novel crystal nonlinearity engineering techniques using sub-coherence-length domains to produce highly pure heralded single photons, significantly improving photon purity for short nonlinear crystals with femtosecond lasers.

Contribution

It presents new methods combining deterministic and probabilistic domain engineering at sub-coherence lengths to enhance photon purity in parametric down-conversion sources.

Findings

Sub-coherence-length domain engineering outperforms previous methods.

Techniques achieve high spectral purity with short nonlinear crystals.

Methods are particularly effective with femtosecond laser sources.

Abstract

Photonic quantum technology relies on efficient sources of coherent single photons, the ideal carriers of quantum information. Heralded single photons from parametric down-conversion can approximate on-demand single photons to a desired degree, with high spectral purities achieved through group-velocity matching and tailored crystal nonlinearities. Here we propose crystal nonlinearity engineering techniques with sub-coherence-length domains. We first introduce a combination of two existing methods: a deterministic approach with coherence-length domains and probabilistic domain-width annealing. We then show how the same deterministic domain-flip approach can be implemented with sub-coherence length domains. Both of these complementary techniques create highly pure photons, outperforming previous methods, in particular for short nonlinear crystals matched to femtosecond lasers.

Figures

Equations

$${|\psi ⟩}_{s,i}=\mathcal{T}\exp \left(2\pi BO\iint d{\omega }_{s}d{\omega }_i\varphi (\Delta k({\omega }_s,{\omega }_i))\alpha ({\omega }_s+{\omega }_i){\hat{a}}_s^{†}({\omega }_s){\hat{a}}_i^{†}({\omega }_i)-h.c.\right){|0⟩}_{s,i},$$

$$\varphi (\Delta k({\omega }_s,{\omega }_i))={\int }_{-\infty }^{+\infty }g(z)e^{i\Delta k({\omega }_s,{\omega }_i)z}dz,$$

$${|{\psi }_{\text{pair}}⟩}_{s,i}=2\pi BO\iint d{\omega }_{s}d{\omega }_{i}f({\omega }_s,{\omega }_i){\hat{a}}_s^{†}({\omega }_s){\hat{a}}_i^{†}({\omega }_i){|0⟩}_{s,i},$$

$$A(z,\Delta k)=-i{\int }_0^{z}g(z^{\prime })e^{i\Delta k{z}^{\prime }}d{z}^{\prime },$$

$$A_{\text{target}}(z,\Delta k)=-i{\int }_0^z{g}_{\text{target}}(z^{\prime })e^{i\Delta k{z}^{\prime }}d{z}^{\prime }.$$

$$e(x+w)=A_{\text{target}}(z+w,\pi /{\ell }_c)-A(z,\pi /{\ell }_c),$$

$$A_{\text{target}}\left(z,\pi /{\ell }_c\right)=-ic\left(\text{erf}\left(\frac{L-2z}{2\sqrt{2}\sigma }\right)-\text{erf}\left(\frac{L}{2\sqrt{2}\sigma }\right)\right),$$

$$E={\left(\sum _{\Delta k}{\left[|{\varphi }_{\text{target}}\left(\Delta k\right)|-|\varphi \left(\Delta k\right)|\right]}^2\right)}^{1/2}\text{.}$$

$$A_m\left(\{s_n{\}}_{n=1}^m,\pi /{\ell }_c\right)=$$

$$=$$

$$e_m(\{s_n{\}}_{n=1}^m)=|A_{\mathrm{target}}(mw)-A_m(\{s_n{\}}_{n=1}^m)|.$$

Full text

I Introduction

The generation of pure single photons on demand is a key requirement for photonic quantum technologies. Single-photon sources can be approximated using either nonlinear parametric processes or single quantum emitters such as quantum dots Lounis and Orrit (2005); Shields (2007); Michler (2009). The most widely used technique for generating single photons is parametric down-conversion (PDC) in nonlinear optical crystals Burnham and Weinberg (1970); Kwiat et al. (1995); Tanzilli et al. (2001); Mosley et al. (2008). In a PDC process, a pump laser is sent through a crystal with large $\chi^{(2)}$ nonlinearity, where $\chi^{(2)}$ is the second order nonlinear electric susceptibility. A pump photon has a small probability of being annihilated while creating pairs of “down-converted” single photons under conservation of energy and momentum. Due to the spontaneous nature of photon-pair creation, the complete generated state is described as a superposition of a large vacuum component, the desired photon pairs, as well as typically undesired higher-order photon-pair events Jennewein et al. (2011); Christ and Silberhorn (2012).

PDC can be used as a heralded source of single photons, in which one photon of a pair is sacrificed in detection to herald the presence of a photon in the other mode. The heralded photon generally emerges in a spectrally mixed state due to strong (anti)correlation in the joint spectrum of PDC photon pairs resulting from the conservation of energy and momentum. To create pure photons, we therefore need to eliminate correlations without introducing mixture in other degrees of freedom.

The easiest way to reduce spectral correlations is to employ narrowband spectral filters. This approach, however, severely compromises the heralding efficiency (i.e. the probability of detecting a photon knowing that the other photon of the pair has been detected) of the source even with ideal filters Meyer-Scott et al. (2017) and decreases the absolute flux of the heralded photons. At high pump powers, narrowband filtering can also introduce mixing in the photon-number degree of freedomBrańczyk et al. (2010).

A lossless method to remove frequency correlations is to use group-velocity matching (GVM) Keller and Rubin (1997); Grice et al. (2001); Mosley et al. (2008); U’Ren et al. (2006). Starting with a set of desired pump and down-conversion wavelengths, one can find phasematching conditions in certain crystals that allow the inverse group velocity of the pump laser in the nonlinear crystal to either match one of the PDC photon’s inverse group velocities or to match the average of the two PDC photons’ inverse group velocities. This erases the timing information that otherwise arises between the pump and PDC photons which in consequence leads to more separable joint spectra and high purity photons. GVM can be achieved both in bulk crystals and in periodically poled crystals, where the crystal’s axes are flipped every coherence length (i.e. the distance over which the phase between the pump and the PDC photons changes by $\pi$) to fulfil the quasi-phase matching condition Fejer et al. (1992), resulting in a periodical configuration of domains with alternating orientation.

What group-velocity matching cannot address is that PDC in standard nonlinear crystals with a rectangular nonlinearity profile produces sinc-shaped joint spectral amplitudes (JSA) which further restricts photon purity Brańczyk et al. (2010). The first step towards overcoming this issue was proposed by Brańczyk et al., who showed that the nonlinearity profile of the PDC crystal could be suitably shaped via domain engineering Brańczyk et al. (2011a) in order to achieve a Gaussian phase matching function (PMF). Considering a poled crystal, it is indeed possible to individually choose each domain orientation such that the effective PMF fits an almost arbitrary target function. A proof-of-principle experiment Brańczyk et al. (2011a) showed a high overlap between the design and the experimentally determined PMF, proving the feasibility of tailored nonlinearities for this purpose.

This technique was subsequently refined Dixon et al. (2013); Dosseva et al. (2016); Tambasco et al. (2016), leading to even better approximations to ideal nonlinearity profiles. Simulations showed that these methods would increase photon purities without compromising the heralding efficiency of the source. Chen et al. Chen et al. (2017) showed that the spectral response of a tailored PDC crystal had a low amount of entanglement in the frequency degree of freedom. However, a direct measurement of the indistinguishability (such as Hong-Ou-Mandel (HOM) interference) of photons produced by two independent PDC processes in engineered crystals has not been done yet. Note that HOM interference between photons generated from a single pulse in one crystal is not a good indicator for single photon purity because one can achieve perfect interference also with a standard periodically poled crystal Hong et al. (1987).

All domain engineering methods developed so far are intrinsically related to a coherence-length structure typical of periodically poled crystals. In Dixon et al. (2013) the domain-periodicity is restricted to twice the coherence length of the crystal; in the other methods Brańczyk et al. (2011a); Dosseva et al. (2016); Tambasco et al. (2016), the domain-width is kept fixed to the coherence length. While this is an intuitive choice and allows a simple approach to PMF shaping, the fixed structure substantially limits the performance and versatility of domain engineering techniques. In this work we move beyond the coherence-length domain boundaries to allow for more fine-grained shaping of the nonlinearity. To this aim we explore two complementary approaches.

In the first approach, we modify a deterministic method for domain engineering Tambasco et al. (2016) and use it for tailoring a crystal with fixed domain widths: we then shift the boundaries of the domains by means of a previously existing annealing method developed for classical applications in higher-harmonic generation Reid (2003); Kornaszewski et al. (2008). In the second approach, we generalise the algorithm in Tambasco et al. (2016) to arbitrarily small, but constant, domain widths (not necessarily equal to the coherence length) and to a complex field amplitude target function.

Both methods lead to a better approximation of the desired phase matching function and to an enhanced heralded single-photon purity, especially in the short-crystal regime.

II Joint spectral amplitude

The output of a second-order nonlinear crystal, pumped by a bright coherent field with envelope function $\alpha(\omega_{s}+\omega_{i})$, can be written as Harder (2016):

[TABLE]

where “$p$” labels the pump field and “$s$” and “$i$” label the two down-converted “signal” and “idler” photons respectively, and $\mathcal{T}$ is the time-ordering operator Brańczyk et al. (2011b); Christ et al. (2013); Quesada and Sipe (2014, 2015). We call $\phi(\Delta k(\omega_{s},\omega_{i}))$ the phase matching function and define it below. The constant $O$ is related to the transverse modes of the photons Harder (2016), and ${B=\chi^{(2)}_{0}(\hbar\omega_{0p}\omega_{0s}\omega_{0i}/\epsilon_{0}\pi^{3}c^{3}n_{p}n_{s}n_{i})^{1/2}}/2$, where $\omega_{0,j}$ (with $j=p,i,s$) are the central frequencies of the fields, $n_{j}$ are the corresponding refractive indices and $\chi^{(2)}_{0}$ is the crystal’s second order nonlinearity. There exist a number of related, but slightly different, definitions of the phase matching function in the literature. We define it as:

[TABLE]

where $\Delta k(\omega_{s},\omega_{i})=k_{p}(\omega_{s}+\omega_{i})-k_{s}(\omega_{s})-k_{i}(\omega_{i})$ is the phase mismatch, which depends on the material dispersion, and $g(z)=\chi^{(2)}(z)/\chi^{(2)}_{0}$ is the normalised nonlinearity along the crystal which will be effectively “tailored” via domain engineering. We note that $g(z)$ is a real function. When defined as above, the PMF for two consecutive domains is the sum of the PMFs for the individual domains111Take for example a block of two domains. If the normalized nonlinearity function for the $i$th domain ($i=1,2$) is $g_{i}(z)$, then the normalized nonlinearity function for both domains is $g_{1,2}(z)=g_{1}(z)+g_{2}(z)$. According to our definition, the PMF for both domains is then $\phi_{12}(\Delta k)=\int_{-\infty}^{+\infty}g_{12}(z)\ e^{i\Delta kz}dz=\int_{-\infty}^{+\infty}g_{1}(z)\ e^{i\Delta kz}dz+\int_{-\infty}^{+\infty}g_{2}(z)\ e^{i\Delta kz}dz=\phi_{1}(\Delta k)+\phi_{2}(\Delta k)$, where $\phi_{i}(\Delta k)$ is the PMF for the $i$th domain. Take on the other hand another popular definition for the PMF, $\phi^{\text{other}}(\Delta k)=L^{-1}\int_{-\infty}^{+\infty}g(z)\ e^{i\Delta kz}dz$. The function $\phi^{\text{other}}(\Delta k)$ is dimensionless, but in this case $\phi^{\text{other}}_{12}(\Delta k)=(2l_{c})^{-1}\int_{-\infty}^{+\infty}g_{12}(z)\ e^{i\Delta kz}dz\neq\phi^{\text{other}}_{1}(\Delta k)+\phi^{\text{other}}_{2}(\Delta k)=l_{c}^{-1}\int_{-\infty}^{+\infty}g_{1}(z)\ e^{i\Delta kz}dz+l_{c}^{-1}\int_{-\infty}^{+\infty}g_{2}(z)\ e^{i\Delta kz}$. (this definition has units of length, while in other popular definitions, the PMF is dimensionless).

The joint spectrum of the down-converted light is given by the product of the PMF and the pump envelope function. We therefore define the joint spectral amplitude as ${f(\omega_{s},\omega_{i})\equiv\phi(\Delta k(\omega_{s},\omega_{i}))\alpha(\omega_{s}+\omega_{i})}$. For the purposes of heralding single photons, we can restrict our analysis to the two-photon state

[TABLE]

given by the first-order term in the Dyson series expansion of Eq. (1).

The spectral purity of the heralded single-photon state is defined as $P=\sum_{k}b_{k}^{4}$ where $b_{k}$ are the Schmidt coefficients defined by the Schmidt decomposition $f(\omega_{s},\omega_{i})=\sum_{k}b_{k}q_{k}(\omega_{s})r_{k}(\omega_{i})$. The purity is therefore limited by correlations in the joint spectral amplitude. We thus seek a separable joint spectral amplitude, that is, a function of the form $f(\omega_{s},\omega_{i})=q(\omega_{s})r(\omega_{i})$. For a symmetric joint spectral amplitude, the standard approach is to first choose group velocities that make $\phi(\Delta k(\omega_{s},\omega_{i}))$ intersect $\alpha(\omega_{s}+\omega_{i})$ at right angles (when plotted as a function of $\omega_{i}$ and $\omega_{s}$), then match the widths of the pump envelope and phase matching functions Laudenbach et al. (2016). This removes most of the correlations, but those that are due to the shape of the PMF remain. Spectral filtering can be used to remove the final correlations, but very high spectral purities can only be achieved at the expense of the heralding efficiency. To achieve high purities and high heralding efficiencies, the phase matching function must be suitably engineered.

III Engineering the phase matching function

In a periodically-poled crystal, $g(z)$ alternates between $+1$ and $-1$ every coherence length: for a large number of periods $g(z)$ can be approximated as constant and the PMF is consequently proportional to the sinc function.

Aperiodic poling is also possible. Although such techniques still constrain $g(z)$ to values of $+1$ and $-1$, the non-trivial structure makes it possible to shape the effective nonlinearity of the crystal and consequently to customise the PMF. In particular, non-trivial poling can make the PMF approximate a Gaussian function, which is necessary to make $f(\omega_{s},\omega_{i})$ separable.

Domain engineering methods have long been studied in nonlinear optics, for example to compress or shape pulses in second harmonic generation Arbore et al. (1997); Imeshev et al. (2000); Shur et al. (2001); Zhang and Gu (2001); Reid (2003); Kornaszewski et al. (2008). These techniques have only recently been adapted for single-photon generation. In this context, existing methods for non-trivial poling fall into two categories. Those that vary the domain widths of a predefined poling configuration Dixon et al. (2013), and those that keep the domain widths equal, but vary their relative orientations Brańczyk et al. (2011a); Dosseva et al. (2016); Tambasco et al. (2016). All methods in the latter category have thus far used domains equal in width to the crystal’s coherence length, where the coherence length $\ell_{c}$ is defined in terms of the phase mismatch $\Delta k(\omega_{s,0},\omega_{i,0})$ at the central PDC frequencies as follows: $\ell_{c}=\pi/\Delta k(\omega_{s,0},\omega_{i,0})$.

In the following, we describe two methods for engineering $g(z)$. Both use the method recently introduced by Tambasco Tambasco et al. (2016) as a starting point, but deviate from the method by allowing domain widths smaller than the crystal’s coherence length. This move toward sub-coherence length structures allows much greater accuracy in tailoring the phase matching function.

III.1 Deterministic domain engineering with fixed domain widths

We begin by summarising the method introduced by Tambasco et al. Tambasco et al. (2016). In this method the width of each domain in the grating is fixed, and the nonlinearity profile is shaped by choosing the orientation of successive “poling blocks” along the crystal, starting with the crystal’s input face (the term “poling block” refers to two consecutive domains of fixed length $w=\ell_{c}$, where $\ell_{c}$ is the coherence length and $w$ is the domain width). The decision to flip (or not flip) a given block is determined by which option gives a closer approximation to the target signal and idler joint field at that point in the crystal. The normalised field inside the crystal at position $z$ is defined as

[TABLE]

where we have omitted the amplitude prefactor Boyd (2003).

The goal is to design a crystal with a $g(z)$ that gives the desired PMF according to Eq. (2). We therefore define a $g_{\text{target}}(z)$, which in turn defines a target field amplitude

[TABLE]

It is sufficient to match the the amplitude at a single value of $\Delta k$. At $\Delta k=\pi/\ell_{c}$, the propagation of the pump for a distance equal to one domain width changes the real part of the generated field amplitude by a factor $\Delta A=\pm 2w/\pi$, where the sign depends on the overall previous domain configuration in the crystal. The imaginary part of the field is always equal to zero at domain boundaries. The poling pattern can be chosen so that the generated field amplitude $A(z,\pi/\ell_{c})$ along the crystal fits as well as possible the target field amplitude $A_{\text{target}}(z,\pi/\ell_{c})$. For simplicity of notation we will refer to a domain pointing up by “up” and to a domain down by “down”.

The poling algorithm proposed by Tambasco et al. involves three different building blocks which are chosen to track the target field amplitude: an up-up block leaves the average field amplitude unchanged, an up-down block increases the average field amplitude by a factor $2\Delta A$, while down-up decrease the average field amplitude by $2\Delta A$. These three possible configurations are graphically represented in Fig. 1 (a).

Tambasco et al. framed their algorithm in terms of blocks to ensure that the inverted regions were of equal width: this is not a necessary requirement—one might just as well choose each individual domain’s orientation. We therefore consider the “domain by domain” error function

[TABLE]

which quantifies the difference between the generated field amplitude at a certain position $z$ and the target field amplitude at $z+w$ (i.e. after one domain). Whether to use a domain up or a domain down is determined by the following conditions:

• •

if $e(z+w)\geq 0$ and $A(z,\pi/\ell_{c})\geq A(z-w,\pi/\ell_{c})$ (i.e. in the previous domain the field amplitude was increasing), flip the domain orientation with respect to the previous domain: in this way the amplitude field will continue to increase;

• •

if $e(z+w)\geq 0$ and $A(z,\pi/\ell_{c})\leq A(z-w,\pi/\ell_{c})$ (i.e. in the previous domain the field amplitude was decreasing), keep the same orientation of the previous domain: in this way the amplitude field will start increasing;

• •

if $e(z+w)<0$ and $A(z,\pi/\ell_{c})\geq A(z-w,\pi/\ell_{c})$, keep the same orientation of the previous domain: in this way the amplitude field will start decreasing;

• •

if $e(z+w)<0$ and $A(z,\pi/\ell_{c})\leq A(z-w,\pi/\ell_{c})$, flip the domain orientation with respect to the previous domain: in this way the amplitude field will continue to decrease.

With this technique, it is possible to generate a wide range of PMF shapes, as long as the corresponding field amplitude does not vary too quickly Tambasco et al. (2016). The computational runtime of this algorithm is minimal (generally of the order of a few seconds).

As discussed above, the case of a Gaussian modulation proportional to $\allowbreak g(z)=\exp(-(z-L/2)^{2}/2\sigma^{2})\cos((\pi/\ell_{c})z)$, where $L$ is the length of the crystal, is the most interesting for engineering a pure heralded photon source. Ignoring the quickly oscillating terms, the resulting target field amplitude is

[TABLE]

where $L$ is the length of the crystal, $\sigma$ is the width of the Gaussian nonlinearity profile and $c$ is a scale factor. As suggested in Tambasco et al. (2016), a value of $\sigma\approx L/4$ allows to reach high single-photon purity without significant degradation in the source brightness. Choosing the prefactor $c=\sqrt{2/\pi}\sigma$ guarantees the optimal PMF: a higher value would compromise the purity because the target amplitude’s gradient is too high to be efficiently traced by the actual field while a lower value would reduce the effective nonlinearity of the crystal.

Notice that the target field amplitude in Eq. (7) is purely imaginary, while the method described above allows only the real part of the field amplitude to be customised. In principle, it is therefore impossible to use this technique to approximate Eq. (7), but in practice, we can approximate Eq. (7) up to a phase factor of $i$. This results in a $\pi/2$ phase shift in the exponent of Eq. (1), but when considering just the two-photon term, it can be regarded as an irrelevant global phase.

Fig. 1 shows a comparison between the algorithm in Tambasco et al. (2016) and our modified algorithm described above. The modified algorithm yields a closer approximation to $A_{\text{target}}$. Consequently, the PMF will be closer to the target and, in the case of the Gaussian shaping, it will lead to higher purities.

III.2 Domain width annealing

The previous methods consider the case of a poled crystal with constant domain width, which allows to easily approach the problem. However, there is no evidence to suggest that a fixed-domain structure leads to an optimal result, and it is reasonable to ask if it is possible to improve the PMF shape by slightly varying the width of each domain. To this aim we use an adapted version of the simulated annealing algorithm introduced by Reid et al. in Reid (2003); Kornaszewski et al. (2008). Annealing algorithms are commonly exploited for finding a global minimum of a given function dependent on multiple parameters Kirkpatrick et al. (1983); Kirkpatrick (1984) by slightly perturbing the system from a suitable starting point, calculating the relative cost function (commonly called energy in analogy with the internal energy in a physical annealing process) and accepting the change with some probability—the higher the energy the lower is the probability of accepting the new configuration. These kinds of algorithms are probabilistic and may require several runs to get an optimal result.

Given a target PMF $\phi_{\text{target}}(\Delta k)$, we first find the initial configuration of domain orientations by means of the method described above. Adjacent domains with the same orientation are grouped together into bigger blocks. Secondly, we define a starting temperature $T$ and a temperature step $\Delta T$ for the annealing algorithm: we found that a value of $T\epsilon\{0.1;10\}$ and $\Delta T=T/100000$ lead to slightly better results. The width of each block is then perturbed by up to 1% and the relative PMF is computed. The perturbation value 1% is empirically determined: too small values would lead to slow convergence of the algorithm while if the perturbation is too big the algorithm is unstable and doesn’t converge to a correct result. It is then possible to find the system energy defined as

[TABLE]

If this energy is smaller than the minimum energy recorded so far the new domain widths are recorded as the best configuration, and they are accepted with a probability of $\exp{\left(-E/T\right)}$ (even if they are not the best configuration). Finally, the temperature is decreased by $\Delta T$ and the algorithm is repeated until the energy becomes smaller than a chosen threshold or the temperature reaches [math]. A more detailed description of the algorithm is provided in Appendix A. This algorithm is computationally demanding and may require a few hours for converging to a solution for a crystal having a large number of domains.

One may think to directly use the annealing procedure without using a pre-processed initial configuration as a seed: however, our annealing algorithm doesn’t allow the flipping of domain orientations, and thus a near-optimal result cannot be obtained without a suitable starting configuration. Moreover, it’s reasonable to think that the method discussed above provides a good initial domain-orientation because we know that the generated amplitude has to fit the ideal case of Eq. (7).

III.3 Domain engineering with arbitrary small sub-coherence length domains

Until now we have considered domain-widths equal to the coherence length with small variations of this configuration. In this section we will see that pushing the algorithm beyond this constraint by allowing a finer discretisation of the domain-structure leads to an even better approximation of the target function.

According to Eq. (5), we consider the field amplitude at the end of $m$-th domain:

[TABLE]

where $w$ is the domain width and $s_{n}=\pm 1$ is the orientation of the $n$-th domain. Note that $A_{m}\left(\{s_{n}\}_{n=1}^{m},\pi/\ell_{c}\right)$ depends on the orientations of all domains that come before it (but not those that come after).

For domain-widths equal to the coherence length, the imaginary part of the field is always zero at the domain boundaries, and therefore it always oscillates about zero (see Fig. 1). If the domain-widths differ from the coherence length, however, the phase might get flipped at a place where the imaginary part is non-zero, providing control over both the real and imaginary parts of the amplitude. With this modification, it is now possible to approximate complex field amplitudes. In particular, the field in Eq. (7) can be approximated up to the correct phase. But to maintain consistency, we will continue to work with the real target function used in the previous section.

To account for the complex nature of the field amplitude, we define a cost function that we want to minimize at each domain:

[TABLE]

Since the target function has a zero imaginary part, the algorithm with the new cost function will force the imaginary component close to zero.

The algorithm can be summarised as follows. First define a domain width $w$, the coherence length $\ell_{c}$, and the number of domains $N$. Then, starting from a crystal having just one domain, compute the cost function $e_{\textsc{up}}$ for the case where a domain up is added and the cost function $e_{\textsc{down}}$ when a domain down is added. Next, compare the two cost functions: if $e_{\textsc{up}}<e_{\textsc{down}}$ keep the configuration where the up domain was added, otherwise keep the configuration where the down domain was added. Repeat for each subsequent domain. A more detailed description of the algorithm is provided in Appendix A. The running time of this algorithm can vary from a several seconds to several minutes depending on the length of the crystal and on the domain-widths.

III.4 Comparison of the algorithms

In this section we benchmark our domain engineering techniques in a realistic scenario. We consider a short potassium titanyl phosphate (KTP) crystal because it can be matched with femtosecond lasers Laudenbach et al. (2016), and its optimal GVM condition (pump wavelength of 791 nm and signal/idler wavelengths of 1582 nm) produces photons compatible with telecom technologies.

We tailor a short KTP crystal ($\sim 2$ mm) with a Gaussian spectral response using three different domain engineering algorithms: the method by Tambasco et al. Tambasco et al. (2016) and our two new methods. The comparison results are presented in Fig. 2. In this regime, the single photon purity increases from a value of $85.4\%$ for a standard periodically poled crystal to $97.3\%$ for the Tambasco algorithm, while a $99.0\%$ photon purity is observed for the width annealing and a value of $99.4\%$ is achieved for arbitrary small sub-coherence domains. A practical limitation for sub-coherence length structures is that poling domains cannot be smaller than a couple of microns. The precision to which such domains can be grown is however very high. These practical aspects and potential tradeoffs coming from uncertainties in domain wall placement and diffusion will be addressed in future work. Here we limited our simulations to domain widths of one tenth of the coherence length.

The degree of separability of the JSA has been estimated through the Schmidt method by using the singular value decomposition Law et al. (2000); Laudenbach et al. (2016). It’s important to notice that this method strictly depends on the spectral range under examination: too small ranges lead to inaccurate results, too wide ranges are computationally intractable. It also depends on the discretisation of the JSA: too course a discretisation make the purity artificially high. An empirically-found reasonable choice for these parameters is a range of about eight times the PMF spectral width, and a JSA discretised as a $100$ by $100$ matrix.

In Fig. 3 we compare the single-photon purity for different methods—our two new methods, and those by Dixon et al. Dixon et al. (2013) and Tambasco et al. Tambasco et al. (2016)—as a function of crystal length (for $\sigma=L/4$). Note that we do not compare with the Dosseva et al. method Dosseva et al. (2016) as this method uses simulated annealing to solve a problem to which the Tambasco et al. method Tambasco et al. (2016) finds the optimal solution. In other words, the purity given by Dosseva et al. (2016) will be bounded from above by the purity given by Tambasco et al. (2016) for any set of parameters.

For all algorithms, the purity is highest for long crystals. The method proposed by Dixon et al. Dixon et al. (2013) has an almost constant value of $97.9\%$ over the range we considered. The purity for all other methods saturates at the same value of $99.5\%$. Note that what limits the purity here is the choice of $\sigma$, not the algorithms themselves. This value can be increased up to $99.9\%$ by taking $\sigma=L/5$ at the expense of a lower overall nonlinearity.

The difference between the algorithms is more drastic for short crystals. While the arbitrary-small domains method always achieves an almost optimal result, the annealing algorithm provides less separable JSAs than the algorithm proposed by Dixon et al. Dixon et al. (2013) for crystals having less than $\sim 60$ domains: this might be due to a sub-optimal initial condition for the annealing process which makes the algorithm converge to a local minimum.

The unmodified Tambasco et al. algorithm Tambasco et al. (2016) is the only one that does not yield purities that monotonically increase with crystal length. For crystals of length $80\ell_{c}$ and $100\ell_{c}$, the purity drops when compared with slightly shorter crystals. This can be explained by considering that the two-domain-block structure of the algorithm doesn’t allow enough flexibility to track the target amplitude closely due to the rough discretisation when the total number of domains is small. Numerical artefacts in the Fourier reconstruction prevented a reliable value for the purity for a $50\ell_{c}$ long crystal to be calculated.

It is finally worth mentioning that the symmetrisation of our algorithms is a straightforward procedure whenever the target phase matching function is symmetric (if it is not the case, then it is not possible to have the same spectral response in both directions): either the annealing and the domain engineering with arbitrary small domains can be applied to the first half of the crystal and reproduced specularly on the second half obtaining a symmetric crystal. This could be necessary in experiments which involve a double pump configuration, where a symmetric crystal may be required for having spectrally indistinguishable photons in both directions.

IV Conclusion

Scalable, fault-tolerant quantum photonics will require photons of the highest purities created with near-unity heralding efficiencies. Spectral filtering severely compromises heralding efficiencies Meyer-Scott et al. (2017), increasing the number of multiplexed PDC sources required for approximating nearly-deterministic single-photon sources Christ and Silberhorn (2012); Francis-Jones and Mosley (2014); Bonneau et al. (2015). At high pump powers narrowband filtering can also introduce undesirable mixing in other degrees of freedom Brańczyk et al. (2010).

Domain engineering is therefore crucial for photonic architectures relying on parametric downconversion. The two sub-coherence length methods we introduced are fast, can readily be implemented commercially, and create higher-purity photons than previous algorithms in particular for short nonlinear crystals matched to femtosecond pump lasers Weston et al. (2016).

It will be a challenge to apply our methods to integrated photonics, where photons are typically created via four-wave mixing. Group-velocity matched four-wave mixing photon sources have been demonstrated but the observed photon purity was again limited by spurious frequency correlations caused by the sinc-shaped biphoton amplitudes Halder et al. (2009); Francis-Jones et al. (2016). Techniques for Gaussian apodization in four-wave-mixing have been suggested Fang et al. (2013) but they require manipulation of the effective material dispersion which complicates simultaneous group-velocity-matching. An ongoing effort aims to introduce techniques which offer a free parameter similar to periodic poling in this platform to which we hope to map our domain-engineering algorithm.

We expect that domain engineering with arbitrarily small domains will have a range of interesting applications in ’classical’ nonlinear optics for which domain engineering was originally designed Arbore et al. (1997); Imeshev et al. (2000); Shur et al. (2001); Zhang and Gu (2001); Reid (2003); Kornaszewski et al. (2008), in particular because we can now keep track of both real and imaginary target functions which might offer advantages for phase-sensitive applications.

Acknowledgements.

Acknowledgments The authors acknowledge Luke G. Helt for providing valuable feedback on the manuscript. This work was supported by the Engineering and Physical Sciences Research Council (grant number EP/N002962/1). F.G. acknowledges studentship funding from EPSRC under grant no. EP/L015110/1. Research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Research and Innovation.

Appendix A: Algorithms

Annealing

Define target phase matching function $\phi(\Delta k)$ and initial domain configuration. 2. 2.

Group together all the adjacent domains in bigger blocks having the same orientation. 3. 3.

Define algorithm temperature T, temperature step $\Delta T$ and energy threshold $E_{t}$ (used for accepting a given configuration as the optimal one). 4. 4.

Compute the energy for the initial configuration and store it in the variable $E_{\text{min}}$. 5. 5.

Apply a random perturbation to the blocks’ width of up to $1\%$ of their current width. 6. 6.

Calculate the new phase matching function. 7. 7.

Compute the corresponding energy $E$. 8. 8.

If $E<E_{\text{min}}$, accept the current block configuration as the best configuration, update $E_{\text{min}}=E$ and go to step 11. 9. 9.

Else if $E>E_{\text{min}}$, accept the new configuration with probability $\exp(-E/T)$. 10. 10.

If the new configuration is not accepted, return to the best configuration an decrement $T$ by $\Delta T$. 11. 11.

Iterate 5 to 10 until $E_{\text{min}}<E_{t}$ or $T=0$.

A block diagram of the algorithm is shown in fig. 4.

Arbitrary-small domains method

Define complex target amplitude profile. 2. 2.

Define parameters $w$, $\Lambda$, the number of domains $N$. 3. 3.

Initialize an empty list $S=\{\}$ for storing the domain orientation. 4. 4.

Define $m=1$. 5. 5.

Create two trial lists that are identical except for the last element: $S_{\textsc{up}}=S+\{\textsc{up}\}$ and $S_{\textsc{down}}=S+\{\textsc{down}\}$ 6. 6.

Compute cost functions for the two trial lists: $e_{\textsc{up}}=e_{m}(S_{\textsc{up}})$ and $e_{\textsc{down}}=e_{m}(S_{\textsc{down}})$. 7. 7.

If $e_{\textsc{up}}<e_{\textsc{down}}$, update $S=S_{\textsc{up}}$. 8. 8.

Else if $e_{\textsc{up}}>e_{\textsc{down}}$, update $S=S_{\textsc{down}}$. 9. 9.

Update $m=m+1$. 10. 10.

Iterate Steps 4 to 9 until $m=N+1$.