Factorization with a logarithmic energy spectrum

In the field of cryptography [2], the complexity of prime number decomposition is crucial for the security of codes. However, the Shor algorithm [3] that takes advantage of entanglement of quantum systems has endangered this security reviving an interest in factorization algorithms and creating the thriving field of quantum information [4]. So far, several experiments [5] have implemented the Shor algorithm but in the meantime other approaches such as Gauss sum factorization [6, 7] have been proposed, and realized experimentally.

At first glance our protocol to factor N = p × q is straightforward and simple: we consider two atoms moving in a potential with logarithmic energy spectrum [8]. Starting from the two-particle ground state, we excite them into a two-particle state of total energy E(N) = E₀ ln N. A measurement of the single-particle energies might result in E₀ ln p or E₀ ln q from which we can determine the factors p and q. In this case the problem would be solved.

Unfortunately, in general the situation is much more complicated. Instead of finding the factors p and q one of the following three scenarios may take place: (i) the atoms may still be found in the two-particle ground state. (ii) Due to the identity N = p × q = 1 × N the two-particle state with energy E(N) is at least two-fold degenerate and we could find the integers 1 and N instead of p and q. (iii) Even worse is the possibility that the excitation may bring the atoms into a two-particle state with a different energy E(N') ≠ E(N). This process is quite probable especially when the density of energy levels around E(N) is large. It is needless to say that integer N' may factor in a completely different way.

In this paper we propose a method that guarantees that after the excitation and the subsequent measurement the atoms are found with a probability of nearly unity in the desired states with single-particle energies E₀ ln p and E₀ ln q. This remarkable feature is a consequence of the fact that only two of all possible two-particle states participate in the quantum dynamics caused by the excitation. The creation of a two-level system consisting of the ground state and the factor state, that is, the state representing the two factors, results from a resonant excitation of the factor state. The other states with this excitation energy have dramatically different transition matrix elements. Indeed, the Rabi oscillations in this two-level system ensure that at half the Rabi period, we end up with the factor state with almost unity probability. Moreover, we can even find an approximate but analytical expression for this time in terms of the parameters of the system. We emphasize that this formula only contains N and is therefore independent of the individual factors. This property allows us to determine the moment at which we have to turn off the excitation and to carry out the energy measurement.

It is interesting to compare and contrast our proposal with the analysis reported in [9], which also addresses the possibilities of factoring contained in logarithmic energy spectra. However, the authors of [9] focus almost exclusively on the limit of a large number of atoms when a description in terms of partition functions is particularly suited. In contrast, we restrict ourselves to two particles. Another major difference between these two approaches lies in the fact that we rely on an explicitly time-dependent interaction, whereas the treatment of [9] by its very nature of using thermodynamical concepts was restricted to a time-independent situation.

This paper is organized as follows. In section 2 we discuss the single-particle states in a potential with a logarithmic energy spectrum and design a time-dependent interaction between two identical bosonic atoms in such a potential which induces a transition preferentially into a two-particle energy state characterized by the factors p and q. Moreover, we introduce in section 3 the two-particle Schrödinger equation governing these transitions.

Section 4 presents our factorization scheme and both, a numerical and an approximate but analytical solution, give insight into the time-dependent probabilities of the two-particle states. Indeed, a numerical example based on typical atomic and trap parameters and neglecting decoherence shows that with our method a correct factorization with more than 90% probability is possible in a single run. However, due to experimental limitations our method is not able to factor arbitrarily large numbers and we carefully discuss these restrictions. A short discussion of the case of more than two prime factors is also addressed. We conclude by summarizing, in section 5, our main results and providing an outlook on possible future developments.

In order to focus on the main ideas while keeping the paper self-contained, we put the details into several appendices. In appendix A we discuss the realization of the one-dimensional (1D) trap and the required effective 1D interaction between the bosonic atoms. In appendix B we show how the dynamics of the two-particle system can be reduced to two levels. Appendix C finally contains an asymptotic treatment of the transition matrix elements needed in our calculations based on the Wentzel–Kramers–Brillouin (WKB) technique [10].

In a recent paper [8], a potential V =V (x) was constructed such that the motion of a non-relativistic particle of mass μ along the x-axis is characterized by a logarithmic energy spectrum

$$\begin{equation} E_\ell=\hbar\omega_0\ln(\ell+1) \end{equation} \tag{ 1 }$$

with $\ell = 0, 1, 2, \ldots$. Here we have introduced the frequency

$$\begin{equation} \omega_0\equiv\frac{E_1-E_0}{\hbar\,\ln2}=\frac{E_1}{\hbar\,\ln2} \end{equation} \tag{ 2 }$$

and the energy is not measured with respect to the minimum of V but to the energy E₀ = 0 of the ground state.

Given the energy spectrum (1) the time-independent Schrödinger equation

$$\begin{equation} {\hat H}_0\,\varphi_\ell(x)\equiv \left(-\frac{\hbar^2}{2\mu}\,\frac{\mathrm{d}^2}{\mathrm{d}x^2}+V(x)\right)\varphi_\ell(x)=E_\ell \varphi_\ell(x) \end{equation} \tag{ 3 }$$

was used to determine iteratively the potential V =V (x) as well as the energy wave functions φ_ℓ = φ_ℓ(x).

In figure 1 we show V together with φ_ℓ for 0 ⩽ ℓ ⩽ 9. In appendix A.1 we describe a method to realize the 1D potential V in a highly elongated cigar-shaped trap.

Zoom In Zoom Out Reset image size

Next we turn to the case of two particles in such a potential. Moreover, for our factorization procedure we need transitions between these two-particle energy eigenstates; that is, we need to introduce a time-dependent interaction between the atoms.

We are considering identical bosonic particles described by the symmetric two-particle energy state

$$\begin{equation} |{uv}\rangle_{\rm B} \equiv \frac{1}{\sqrt{2}}\,\left(|{uv}\rangle+|{vu}\rangle\right),\quad u\ne v, \end{equation} \tag{ 4 }$$

$$\begin{equation} |{uu}\rangle_{\rm B}\equiv|{uu}\rangle. \end{equation} \tag{ 5 }$$

Because these states are stationary for non-interacting particles, we have to introduce an interaction ${\hat H}_{\mathrm {W}}$ which we model by a 1D Fermi contact potential [11, 12]

$$\begin{equation} {\hat H}_{\rm{1D}}\equiv\gamma\,\delta(x_1-x_2) \end{equation} \tag{ 6 }$$

of strength γ, where x₁ and x₂ are the positions of the atoms. We emphasize that quantum gates based on the time-independent Hamiltonian (6) have been suggested [13] and realized experimentally [14].

In order to induce transitions between the states, we choose a sinusoidal interaction

$$\begin{equation} {\hat H}_{\rm{W}}\equiv\gamma\,\sin(\omega_{{\rm ext}}t)\,\delta(x_1-x_2) =\sin(\omega_{{\rm ext}}t)\,{\hat H}_{\rm 1D}, \end{equation} \tag{ 7 }$$

where the external frequency ω_ext may be adjusted in such a way as to cause the desired transitions. Moreover, our control over the interaction allows us to switch it on and off.

In appendix A.2 we show that the interaction (7) can be realized experimentally for bosonic atoms in the neighbourhood of a Feshbach resonance [15]. The dependence of the 1D interaction strength γ on both, the s-wave scattering length and the three-dimensional (3D) geometry of the trap, is given in appendix A.3.

The matrix elements of the interaction Hamiltonian given by (6) with the real-valued single-particle eigenfunctions φ_j read

$$\begin{equation} \langle{rs}|{\hat H}_{\rm 1D}|{uv}\rangle\equiv \gamma \, \int \mathrm{d}x\,\varphi_r(x)\varphi_s(x) \varphi_u(x)\, \varphi_v(x)\equiv \gamma \kappa_{rs,uv}. \end{equation} \tag{ 8 }$$

For bosons, a simple calculation using (4) and (5) gives

$$\begin{equation} \left.{_{\rm B} \langle{rs}|{\hat H}_{\rm 1D} |{uv}\rangle_{\rm B} }\right. =2\,\langle{rs}|{\hat H}_{\rm 1D} |{uv}\rangle, \end{equation} \tag{ 9 }$$

$$\begin{equation} \left.{_{\rm B} \langle{rr}|{\hat H}_{\rm 1D}|{uv}\rangle_{\rm B} }\right. =\sqrt{2}\,\langle{rr}|{\hat H}_{\rm 1D} |{uv}\rangle, \quad \left.{_{\rm B} \langle{rs}|{\hat H}_{\rm 1D}|{uu}\rangle_{\rm B} }\right. =\sqrt{2}\,\langle{rs}|{\hat H}_{\rm 1D} |{uu}\rangle, \end{equation} \tag{ 10 }$$

$$\begin{equation} \left.{_{\rm B} \langle{rr}|{\hat H}_{\rm 1D}|{uu}\rangle_{\rm B} }\right. =\,\langle{rr}|{\hat H}_{\rm 1D} |{uu}\rangle, \end{equation} \tag{ 11 }$$

where we have assumed r ≠ s and u ≠ v. These matrix elements form the backbone of the quantum dynamics described in the next section.

Suppose that the integer N = pq is a product of two primes p and q and our goal is to find p and q. Our algorithm starts from two non-interacting atoms in the ground state |0,0〉 of our potential with logarithmic energy spectrum (1). At time t = 0 the time-dependent interaction ${\hat H}_{\mathrm {W}}$ given by (7) is switched on and the frequency ω_ext of the harmonic time dependence (7) is chosen to be resonant with the transition into a state of the non-interacting two-particle system with total energy ℏω₀ ln N leading to the condition

$$\begin{equation} \omega_{\rm{ext}}=\omega_0 \ln N= \omega_0(\ln p + \ln q). \end{equation} \tag{ 16 }$$

At some time T we switch off the interaction and a subsequent measurement of the energy of one particle results in ℏω₀ ln p'. If the quotient q' = N/p' is an integer not equal to N or 1 we have found both prime factors p = p' and q = q'. We shall show that there exists a time when the probability of finding the two-particle state which contains the factors dominates all the other ones. In addition, we shall give an approximate but analytical expression for this time. We emphasize that in this algorithm we only make a single measurement on a single particle to find both factors.

In appendix B we show that the complete system of equations (14) can be approximated by the two equations (B.7) which couple the probability amplitudes for the two-particle ground state |0,0〉 and the factor state |p − 1,q − 1〉. In the present section, we derive an approximate but analytic expression for the time T where the probability |Ψ_p−1,q−1(T)|² for the two atoms to be in the factor state is large.

We introduce in (B.7) the dimensionless time τ ≡ ω₀t as well as the dimensionless matrix element

$$\begin{equation} {\tilde \kappa}_{rs,uv}\equiv \frac{\gamma}{2\hbar\omega_0}\,\kappa_{rs,uv} \end{equation} \tag{ 17 }$$

defined by (8) and arrive at the system

$$\begin{eqnarray} &&{\dot {\bar \Psi}}_{0,0}(\tau) = -{\tilde \kappa}_{00,p-1,q-1} {\bar \Psi}_{p-1,q-1}(\tau), \nonumber\\ &&{\dot {\bar \Psi}}_{p-1,q-1}(\tau) = {\tilde \kappa}_{00,p-1\,\,q-1} {\bar \Psi_{0,0}}(\tau) \end{eqnarray} \tag{ 18 }$$

for the approximate probability amplitudes ${\bar \Psi }_{0,0}$ and ${\bar \Psi }_{p-1,q-1}$ For the initial conditions ${\bar \Psi }_{0,0}(0)=1$ and ${\bar \Psi }_{p-1,q-1}(0)=0$, we immediately obtain the solutions

$$\begin{equation} {\bar \Psi}_{0,0}(\tau)=\cos \left(\Omega \tau\right) \mbox{ and } {\bar \Psi}_{p-1,q-1}(\tau)=\sin\left( \Omega \tau \right), \end{equation} \tag{ 19 }$$

which describe Rabi oscillations between the ground state and the factor state with the dimensionless frequency

$$\begin{equation} \Omega\equiv {\tilde \kappa}_{00,p-1\,\,q-1}. \end{equation} \tag{ 20 }$$

The corresponding probabilities $|{\bar \Psi }_{0,0}(t)|$ and $|{\bar \Psi }_{31-1,41-1}(t)|$ for the example of N = 1271 = 31 × 41 discussed in section 4.3.1 are displayed in figure 2 by dashed lines. These curves show that the probability $|{\bar \Psi }_{31-1,41-1}(t)|$ reaches unity for the time t = π/(2Ωω₀). We study a more realistic example in section 4.3.1 below and find that if we would switch off the interaction at this time the system is found in the factor state with a high probability. Unfortunately, Ω still contains via the matrix element ${\tilde \kappa }_{00,p-1\,\,q-1}$ the unknown factors p and q. We overcome this dilemma by the relation

$$\begin{equation} \Omega \approx \frac{\gamma\alpha}{\hbar \omega_0} \,\,\frac{ 1}{\pi (2N\ln N)^{1/2}} \end{equation} \tag{ 21 }$$

following from the definitions (17) and (20) of the dimensionless matrix element ${\tilde \kappa }_{00,p-1\,\,q-1}$ and of Ω, respectively, the expression (C.6) for α and the asymptotic expansion (C.22) of ${\tilde \kappa }_{00,p-1\,\,q-1}$. Note again that due to the bosonic nature (10) of our atoms we have multiplied (C.22) by a factor of 2^1/2 assuming that p ≠ q.

Zoom In Zoom Out Reset image size

We emphasize that in the expression (21) the factors p and q do not appear anymore. Hence, we can approximate the optimal time T_opt by

$$\begin{equation} T_{\rm opt}\approx \frac{\hbar}{\gamma\alpha}\,\,\frac{\pi^2}{\sqrt{2}}\left( N \ln N \right)^{1/2} =\frac{\pi}{4\sqrt{2}} \,\frac{1}{a_{\rm max}\alpha}\,\frac{2\pi}{\omega_\perp} \left(N\,\ln N\right)^{1/2}, \end{equation} \tag{ 22 }$$

where in the last step we have recalled the expression (A.19) for γ.

In the preceding section we have neglected all off-resonant transitions. To test the consequences of this approximation, we include now a finite number M of the off-resonant transitions

$$\begin{equation} |{0,0}\rangle\rightarrow |{p'-1,q'-1}\rangle \end{equation} \tag{ 23 }$$

nearest to the factor state with p'q' = N ± 1,N ± 2,... in (18). According to (C.19) the matrix elements of the transitions (23) have magnitudes comparable to (B.6) if M ≪ N. We arrive at the (M + 2)-dimensional system

$$\begin{eqnarray} &&{\dot {\Psi}}_{0,0}(\tau) = -{\tilde \kappa}_{00,p-1\,\,q-1}\,\Psi_{p-1,q-1}(\tau) -\sum_{k=1}^M {\tilde \kappa}_{00,P_k-1\,\,Q_k-1}\, \mathrm{e}^{\mathrm{i}\Delta \omega_k\tau } \,\Psi_{P_k-1,Q_k-1}(\tau),\nonumber\\ &&{\dot {\Psi}}_{p-1,q-1}(\tau) = {\tilde \kappa}_{p-1\,\,q-1,00}\,\Psi_{0,0}(\tau),\nonumber\\ &&{\dot {\Psi}}_{P_k-1\,\,Q_k-1}(\tau) = {\tilde \kappa}_{P_k-1\,\,Q_k-1,00}~\mathrm{e}^{-\mathrm{i}\Delta \omega_k \tau}\, \Psi_{0,0}(\tau) \end{eqnarray} \tag{ 24 }$$

with the dimensionless frequency Δω_k ≡ ln N − ln (P_k Q_k) and k = 1,...,M.

One might wonder why only in the first of equations (24) additional terms corresponding to off-resonant transitions appear. The explanation for this approximation is as follows. Firstly, we consider 'upward' transitions from the factor state into states with energy ≈ℏω₀ ln N². As shown in appendix C their matrix elements scale with N⁻¹ and therefore we neglect them. Secondly, we consider 'downward' transitions from the factor state into the ground state and states nearby. The nearest state |0,1〉 has an energy of ℏω₀ ln 2. Consequently, the term

$$\begin{eqnarray*} &&\kappa_{p-1\,\,q-1,0\,1} \,\mathrm{e}^{-\mathrm{i} \tau \ln 2 }\Psi_{0,1}(\tau) \end{eqnarray*}$$

oscillates much faster than the oscillating terms in the first equation and following the arguments of appendix B is therefore neglected.

The linear system (24) has time-dependent coefficients but with the substitution

$$\begin{eqnarray*} &&\mathrm{e}^{\mathrm{i}\Delta \omega_k \tau}\,\Psi_{P_k-1,Q_k-1}(\tau) \equiv {\tilde \Psi}_{P_k-1,Q_k-1}(\tau) \end{eqnarray*}$$

it reduces to a linear system with time-independent coefficients which can be solved by a standard numerical procedure.

4.3.1. Example

We consider the case N = 1271 = 31 × 41 and solve (24) subject to the initial condition

$$\begin{equation} \Psi_{00}(0)=1 \quad {\rm and }\quad \Psi_{rs}(0)=0 \end{equation} \tag{ 25 }$$

for all r + s > 0 taking into account the M = 8 nearest off-resonant transitions. There is no transition into the two-particle state with N' = 1270 because all matrix elements e.g. ${\tilde \kappa }_{00,10-1\,\,127-1}$ vanish by the symmetry of the integrand in (8). But for N' = 1272 four transitions contribute and are taken into account as well as three transitions for N' = 1269 and a single transition for N' = 1273. This results in a system of 20 real equations for 20 unknown real functions since the probability amplitudes are complex.

In section A.3 of the appendix, the interaction strength γ = 2ℏω_⊥a_max is expressed by the classical frequency ω_⊥ of a transverse oscillation in the quasi-1D trap and by the amplitude a_max of the oscillating scattering length of the bosonic atoms. As a result the dimensionless matrix element ${\tilde \kappa }$ given by (17) translates into

$$\begin{eqnarray*} &&{\tilde \kappa}_{rs,uv}=\frac{\omega_\perp}{\omega_0}\,a_{\rm max}\,\kappa_{rs,uv}. \end{eqnarray*}$$

For our example we consider two $\rm \left .^{85}Rb\right .$ atoms with scattering length a₀ = 127 nm [12] in a quasi-1D trap of length 0.2 cm, aspect ratio ω_⊥/ω₀ = 10, with longitudinal frequency ω₀ = 2π × 1000 Hz and the amplitude a_max = 1 nm of the oscillating scattering length. These parameters together with the asymptotic expressions (C.19) for the matrix elements κ_00,uv derived in appendix C complete the system (24) of equations. Figure 2 shows the time dependence of the occupation |Ψ_0,0(t)|² of the ground state |0,0〉 starting from unity and |Ψ_{31−1,41−1}(t)|² of the factor state |31 − 1,41 − 1〉 starting from zero due to the initial conditions (25), both in solid lines. A Rabi oscillation between these states is clearly seen. Moreover, figure 2 shows that the probabilities |Ψ_{Pk − 1,Qk − 1}(t)|² for the off-resonant transitions are small. They are displayed in figure 3 once more on a different scale.

Zoom In Zoom Out Reset image size

If we switch off the interaction at T_opt indicated by a vertical line in figure 2, the probability that we find the two atoms in the state |31 − 1,41 − 1〉 is 93%, whereas for |0,0〉 it is only 1.7%. The remaining probability is distributed over the states resulting from off-resonant transitions.

We recall that in (B.2) we have derived the condition ωT/N ≫ 2 for the validity of the secular approximation leading to the system (18). Using the parameters of the present example, we find that ω₀T_opt/N = 9.0, which is larger than two. This is consistent with the result that the transition probabilities calculated from (18) and (24), respectively, and displayed in figure 2 agree reasonably.

The present example is typical in the following sense. The asymptotic condition N ≫ 1 is clearly fulfilled. But the number to be factored is also limited by the experimental conditions as shown in section 4.4. We conclude that the parameters of the example allow us to find the factors of any integer N composed of two prime factors chosen from a range of about 10²–10⁴ with a high probability. In section 4.5 we sketch the cases of a different number of prime factors.

4.3.2. Detuned external frequency

Because the energy levels in the logarithmic spectrum (1) become closer and closer spaced as the quantum number ℓ increases, it is obvious that the external frequency must meet the resonance condition ω_ext = ω₀ ln N with a high accuracy. Here we demonstrate that even a small detuning of the external frequency changes the probability of finding the factor state at t = T_opt drastically.

We return to our example of section 4.3.1. At N = 1271 we find a spacing ΔΩ = ω₀Δω = 5.94 s⁻¹ between the levels N and N + 1. This has to be contrasted with the resonant frequency ω_ext = ω₀ ln N = 4.5 × 10⁴ s⁻¹.

In figure 4 the probability of finding the factor state at time t = T_opt is shown as a function of the detuning δω of the external frequency ω_ext = ω₀ ln N + δω as calculated from the system (24). It can be seen that a precision of significantly less than 1 s⁻¹ of the external frequency ω_ext is required if a single run of the experiment results in the factors with a probability of more than 90%.

Zoom In Zoom Out Reset image size

We now analyse the restrictions on the largest integer N that can be factored by our method which are imposed by the experimental limitations and the model. The most severe condition is dictated by decoherence.

4.4.1. Trap parameters

We begin by considering the properties of the trap and identify two parameters which limit the number N to be factored by the protocol described in this paper. The first one is the aspect ratio ω_⊥/ω₀ which due to assumption (A.4) to have a 1D trap takes the form

$$\begin{equation} N< \mathrm{e}^{\omega_\perp/\omega_0}. \end{equation} \tag{ 26 }$$

Since this condition involves an exponential it is not the limiting factor.

The length L of a trap is limited by technical reasons and represents the second parameter constraining the number N. According to appendix C we estimate the distance L_N between the two classical turning points x_N of a particle with energy E_N by

$$\begin{equation} L_N=2|x_N|\approx \sqrt{2\pi}\,N\,\alpha^{-1} \end{equation} \tag{ 27 }$$

which yields with the expression (C.6) for α the inequality

$$\begin{equation} N< \frac{\alpha L}{\sqrt{2\pi}}, \end{equation} \tag{ 28 }$$

which is an upper bound on N. This condition is more serious than (26) since it scales only linearly.

4.4.2. Off-resonant transitions

The third parameter limiting N is the strength γ of the Fermi contact interaction (6). Indeed, the frequency ω_ext is chosen to be resonant with the transition |0,0〉 → |p − 1,q − 1〉 where pq = N. We now ask under what conditions transitions into neighbouring states |P − 1,Q − 1〉 with PQ = N ± 1,N ± 2,... have small probabilities.

Consider, for example, the transition |0,0〉 → |P − 1,Q − 1〉 with PQ = N' ≡ N − 1 driven by a frequency ω_ext = ω₀ ln N, which is off-resonant by a detuning

$$\begin{eqnarray*} &&\Delta\omega\equiv\omega_0\,(\ln N-\ln N')\approx \frac{\omega_0}{N}. \end{eqnarray*}$$

Since Δω decays inversely with N we recognize that with increasing N neglecting off-resonant transitions becomes more and more a problem. If the transition into |P − 1,Q − 1〉 with PQ = N − 1 were the only one possible for the system, we can apply the theory of Rabi oscillations presented in many textbooks [10]. According to this theory, the amplitude of the Rabi oscillation of the transition probability approaches zero as the ratio ℏΔω/|W_if| of the energy ℏΔω corresponding to the detuning and the transition matrix element W_if of the perturbation causing the transition approaches infinity, that is,

$$\begin{eqnarray*} &&1\ll \frac{\hbar \Delta\omega}{|W_{\mathrm{if}}|}. \end{eqnarray*}$$

Of course, there are many more transitions of our two-particle system but the frequency of the off-resonant transition discussed above is the one closest to the frequency ω_ext of the perturbation. We now give a rough estimate of how the number N to be factored is limited by the parameters of the trap as well as the interaction strength γ. With the help of the expression (A.16) for γ and the WKB formula (C.22) for the matrix element $W_{\mathrm {if}}=\langle {0,0}|{\hat H}_{\rm {1D}}|{p-1,q-1}\rangle$, we obtain the inequality

$$\begin{equation} \frac{N}{\ln N}<\left(\frac{\pi}{2\sqrt{2}}\,\frac{\omega_0}{\omega_\perp}\frac{1}{a_{\rm max} \alpha}\right)^2. \end{equation} \tag{ 29 }$$

We note that the parameters chosen for the example of section 4.3.1 fulfil the conditions (26)–(28).

4.4.3. Decoherence

Up to now we have not addressed the issue of decoherence. Indeed, it is completely neglected in the system of equations (24). Therefore, decoherence must be negligible at least for the duration T_opt of the interaction.

With the help of (22) we cast (29) into the form

$$\begin{eqnarray*} &&N<\frac{1}{\pi}\,\omega_0 T_{{\rm opt}}\leqslant \frac{1}{\pi}\,\omega_0 T_{\rm dec}, \end{eqnarray*}$$

which shows that the number to be factored is limited not only by the trap geometry but, most crucially, by the time T_dec the system can be maintained free of decoherence. In our example in section 4.3.1 we find that T_opt = 1.83 s.

Energy levels within neutral atoms in optical lattices are reported [17] to have long decoherence times typically in the range of seconds. Instead, in our model the controlled interaction (7) causes transitions between vibrational states of the atoms in the trap. A theoretical work [18] very similar to ours on vibrating neutral atoms in 1D optical microtraps estimates times up to 1 s, which hopefully may improve as technology develops.

Finally, we briefly discuss three cases where the number of prime factors is different from two. To perform at least an order of magnitude evaluation of the probabilities of finding the prime factors, we assume all of the factors to be of approximately the same size.

We start with the case of three different prime factors N = pqr. Here we must apply the scheme of section 4.1 repeatedly as we now demonstrate. Acting as before with an external frequency ω₀ ln N onto the two-particle ground state, it is highly probable that one of the excited states |p − 1,qr − 1〉, |q − 1,pr − 1〉 or |r − 1,pq − 1〉 gets populated, each with a probability of roughly 1/3. Any single-particle energy measurement can lead to one of the six integers contained in these states and if an energy ℏω₀ ln ℓ₁ is measured, then the integers ℓ₁ and ℓ₂ = N/ℓ₁ are calculated. To find all three factors p, q, r a second run is needed. If here any of the four remaining integers results all three factors can be determined. But if one of the two integers ℓ₁ or ℓ₂ already known reappears, which happens with a probability of about 1/3, then we ran in the same two-particle state as before and a third run is needed and so on. We estimate the probability of full success with n repeated runs to be of the order 1 − 1/3ⁿ⁻¹.

Secondly, consider the case of only one factor, i.e. N is prime, but assume that this is not known to us and we try to apply the protocol of section 4.1 assuming erroneously N would be a factor of two primes. But provided that N is odd only a transition into the state |00,0 N − 1〉 is possible. As discussed in appendix C the Rabi frequency κ_00,0 N−1 given by (C.24) is about a factor of N^−1/2 smaller than that calculated from (C.22). At the time t = T_opt given by (22) the system would with high probability still be in the ground state, and from a measurement resulting in the single-particle ground state it can be concluded that N is prime.

Our last example demonstrates that even with a moderate number of prime factors our protocol becomes impractical because it needs a number of repeated runs. Consider the factorization of N = 1275 = 3 × 5² × 17. From these prime factors five different pairs of integers can be formed whose product equals N = 1275, such as e.g. 15 and 85. We leave it to the reader to convince himself that at least five runs are needed to find the four prime factors and the probability of geting all of them with just five runs has an order of magnitude of about 4%. This means that typically many more runs are needed to find the prime factors.

We consider a rotationally symmetric potential

$$\begin{equation} V_{\rm{3D}}(x,y,z)\equiv V(x)+\frac{1}{2}\mu\omega_\perp^2\,(y^2+z^2), \end{equation} \tag{ A.1 }$$

which is harmonic in the y–z-direction, while in the x-direction the potential V =V (x) gives rise to the logarithmic spectrum (1).

We realize V_3D by an optical potential formed by a holographic mask which can create a desired intensity distribution I = I(r) of an electromagnetic field in position space [24]. If the frequency of this radiation is far-detuned from any transition frequency of the atom, the resulting potential V_3D = V_3D(r) is proportional to I(r) [25]. Hence, it is possible with today's technology to tailor the potential and, in particular, create one that gives rise to a logarithmic energy spectrum.

The energy eigenvalues of a particle of mass μ moving in the potential V_3D defined by (A.1) read

$$\begin{equation} E_{\ell mn}=\hbar \omega_0 \ln(\ell+1) +\hbar \omega_\perp (m+n+1) \end{equation} \tag{ A.2 }$$

with the eigenfunctions

$$\begin{equation} \psi_{\ell mn}(x,y,z)=\varphi_\ell (x)\, u_m(y)\, u_n(z) \end{equation} \tag{ A.3 }$$

with ℓ,m,n = 0,1,2,.... Here u_m = u_m(y) denote the well-known energy wave functions of a harmonic oscillator with mass μ and eigenfrequency ω_⊥, and φ_ℓ(x) are the eigenfunctions of the Hamiltonian ${\hat H}_0$ discussed in section 4.3.1.

Our quasi-1D model is defined as usual [11]: the transverse quantum numbers m and n of the particle are confined to those of the ground state, that is, to m = n = 0. As long as the longitudinal energy E_ℓ is below the energy ℏω_⊥, that is,

$$\begin{equation} \hbar\omega_0 \ln (\ell+1) < \hbar\omega_\perp, \end{equation} \tag{ A.4 }$$

no transition to an excited transverse state is possible and the energy of the particle is characterized by the longitudinal quantum number ℓ and the trap can be considered to be longitudinal.

To model a short-range interaction in three dimensions, one cannot resort to a delta function because it does not influence the dynamics in a scattering process. For pure s-wave scattering, however, the Schrödinger equation can be solved for a pseudo potential acting on a spherically symmetric wave function ψ(r) as

$$\begin{equation} V_{{\rm pseudo}}({\bf r})\psi(r)=\gamma_{{\rm 3D}} \delta^{(3)}({\bf r})\frac{\partial}{\partial r}\left(r \psi(r)\right), \end{equation} \tag{ A.5 }$$

which mimics a hard-core potential with an interaction constant

$$\begin{equation} \gamma_{\mathrm{3D}}\equiv \frac{4\pi\hbar^2 a}{\mu} \end{equation} \tag{ A.6 }$$

proportional to the s-wave scattering length a [26]. Acting with the pseudo potential (A.5) on a function ψ_reg(r) regular at the origin, it simply follows that

$$\begin{equation} V_{{\rm pseudo}}({\bf r})\psi_{\rm{reg}}({\bf r})=\gamma_{{\rm 3D}} \delta^{(3)}({\bf r})\frac{\partial}{\partial r}\left(r\psi_{{\rm reg}}({\bf r})\right)= \gamma_{{\rm 3D}} \delta^{(3)}({\bf r})\psi_{{\rm reg}}({\bf r}) \end{equation} \tag{ A.7 }$$

(note that δ(x)f(x) = δ(x)f(0)), i.e. here the pseudo potential (A.5) may be substituted by a contact potential.

In this paper we do not solve a Schrödinger equation for the scattering of the two atoms, but instead in our perturbative approach we only calculate matrix elements of the interaction potential with the unperturbed wave functions which are certainly regular everywhere in space. Therefore we are allowed to use the contact potential

$$\begin{equation} {\hat H}_{{\rm 3D}}\equiv \gamma_{{\rm 3D}}\delta^{(3)}({\bf r}_1-{\bf r}_2) \end{equation} \tag{ A.8 }$$

in our further calculations. This contact potential is valid for very long wavelengths of the interacting cold atoms.

If the atoms exhibit a Feshbach resonance [15], it is possible to tune the interaction (A.8) between the atoms via the scattering length

$$\begin{equation} a(B)=a_0\left(1+\frac{\Delta B}{B-B_0}\right), \end{equation} \tag{ A.9 }$$

which near the resonance B = B₀ depends [12] on a magnetic induction B as shown in figure A.1. The parameters B₀ and ΔB are determined by the microscopic structure of the atoms [12] but we treat them here simply as phenomenological parameters.

Zoom In Zoom Out Reset image size

For a time-dependent magnetic induction B of the form

$$\begin{equation} B(t) =B_0-\Delta B \,\left(1+\frac{a_{\rm{max}}}{a_0}\,\sin (\omega_{\rm{ext}}t)\right) \end{equation} \tag{ A.10 }$$

with a_max ≪ a₀ the scattering length,

$$\begin{equation} a(t)\approx a_{\rm{max}} \sin (\omega_{\rm{ext}} t) \end{equation} \tag{ A.11 }$$

varies sinusoidally with t.

Hence, we have realized the two-atom interaction (7) by a time-dependent magnetic induction B(t) oscillating with frequency ω_ext around B = B₀ − ΔB where the interaction (7) vanishes.

We now calculate by following [27] the matrix elements ${\cal M} \equiv \langle {r00,s00}| {\hat H}_{{\mathrm { 3D}}} |{u00,v00}\rangle$ of the 3D contact potential (A.8) with the unperturbed wave functions 〈xyz|ℓ00〉 ≡ ψ_ℓ00(r) defined by (A.3) using the ground state wave function

$$\begin{equation} u_0(y)=\left(\alpha_\perp^2/\pi\right)^{1/4}\,\exp(-\alpha_\perp^2y^2/2) \end{equation} \tag{ A.12 }$$

with α²_⊥ = μω_⊥/ℏ.

From the definition

$$\begin{equation} {\cal M}=\gamma_{\rm{3D}}\int \mathrm{d}^3 r_1\, \mathrm{d}^3 r_2\, \psi_{r00}({\bf r}_1) \psi_{s00}({\bf r}_2)\delta^{(3)}({\bf r}_1-{\bf r}_2) \psi_{u00}({\bf r}_1) \psi_{v00}({\bf r}_2) \end{equation} \tag{ A.13 }$$

of ${\cal M}$ and from the expression (A.3) for ψ_ℓmn, we find the formula

$$\begin{equation} {\cal M}=\gamma_{\rm{3D}} \left[\int \mathrm{d}y \, u_0(y)^4\right]^2\,\int \mathrm{d}x \,\varphi_r(x) \varphi_s(x)\varphi_u(x)\varphi_v(x). \end{equation} \tag{ A.14 }$$

After performing the integration of the ground state wave functions (A.12) over y, we find that

$$\begin{equation} {\cal M}=\gamma_{\rm{3D}}\,\frac{\alpha_\perp^2}{2\pi}\,\int \mathrm{d}x\, \varphi_r(x) \varphi_s(x)\varphi_u(x)\varphi_v(x). \end{equation} \tag{ A.15 }$$

Comparison with the 1D matrix element

$$\begin{equation} \langle{rs}| {\hat H}_{{\rm 1D}}|{uv}\rangle=\gamma\,\int \mathrm{d}x\,\varphi_r(x)\varphi_s(x)\varphi_u(x)\varphi_v(x) \end{equation} \tag{ A.16 }$$

gives the expression

$$\begin{equation} \gamma= \frac{\alpha_\perp^2}{2\pi}\,\gamma_{\rm{3D}}=2\hbar\omega_\perp a, \end{equation} \tag{ A.17 }$$

where we have used (A.6).

When we recall the relation (A.11) for the time-dependent scattering length, we find that

$$\begin{equation} 2\hbar\omega_\perp a(t)=2\hbar\omega_\perp a_{\rm max} \sin(\omega_{{\rm ext}}t) \, \, =\gamma \sin(\omega_{\rm ext}t), \end{equation} \tag{ A.18 }$$

which defines the effective 1D coupling constant

$$\begin{equation} \gamma\equiv 2\hbar\omega_\perp a_{\rm max}. \end{equation} \tag{ A.19 }$$

This relation may also be found in [11, 28], but in a different form and with a different derivation.

So far we have considered distinguishable particles. However, in the bosonic case these matrix elements have to be multiplied by the appropriate factors summarized in (9)–(11).

We start with the well-known WKB approximation [10]

$$\begin{equation} \varphi_\ell^{\rm{(WKB)}}(x)\equiv 2\left(\frac{\mu}{T_\ell p_\ell(x)}\right)^{1/2}\, \cos \left[\frac{1}{\hbar}\,\int_x^{x_\ell} \mathrm{d}x'\,p_\ell(x')-\frac{\pi}{4}\right] \end{equation} \tag{ C.1 }$$

of the ℓth single-particle energy wave function valid inside the oscillatory regime and far away from both classical turning points ± x_ℓ defined by the condition

$$\begin{equation} E_\ell=V(\pm\, x_\ell). \end{equation} \tag{ C.2 }$$

When we insert the Rydberg–Klein–Rees potential [29] applied to the logarithmic spectrum calculated in [8] into the classical momentum

$$\begin{equation} p_\ell(x)=\sqrt{2\mu}\,\sqrt{E_\ell-V(x)}, \end{equation} \tag{ C.3 }$$

we find the Bohr–Sommerfeld–Kramers quantization rule

$$\begin{equation} S\equiv 2 \int_{-x_\ell}^{x_\ell} \mathrm{d}x\, p_\ell(x) =2 \pi\hbar \left(\ell+\frac{1}{2}\right). \end{equation} \tag{ C.4 }$$

Having in mind that we need the WKB wave functions (C.1) for large indices ℓ only, we use the approximate asymptotic potential [8]

$$\begin{equation} V(x)\sim\hbar\omega_0\,\ln\left(\sqrt{\frac{2}{\pi}}\alpha|x|\right) \end{equation} \tag{ C.5 }$$

with

$$\begin{equation} \alpha\equiv \sqrt{\frac{\mu\omega_0}{\hbar}}, \end{equation} \tag{ C.6 }$$

to determine the classical turning points

$$\begin{equation} x_\ell=\pm \sqrt{\frac{\pi}{2}}\,(\ell+1)\,\alpha^{-1}. \end{equation} \tag{ C.7 }$$

The classical period

$$\begin{equation} T_\ell = \frac{\partial S}{\partial E_\ell}=\frac{\partial S}{\partial \ell} \frac{\partial \ell}{\partial E_\ell} =2\pi\,(\ell+1)\,\omega_0^{-1} \end{equation} \tag{ C.8 }$$

is calculated using the action given by (C.4) as well as the energy spectrum (1).

To obtain the scaling properties of our factorization scheme and to completely avoid numerical evaluations, we observe that in the matrix elements

$$\begin{equation} \kappa_{00,p-1\,\,q-1} = \int_{-\infty}^{\infty} \mathrm{d}x\, \varphi_0(x)^2\,\varphi_{p-1}(x)\varphi_{q-1}(x) \end{equation} \tag{ C.9 }$$

used in section 4.2 the square of the ground state wave function φ₀ appears. Because of its limited range shown in figure C.1 the wave functions φ_ℓ(x) with ℓ = p − 1 and ℓ = q − 1 are needed only for values of x close to the origin.

Zoom In Zoom Out Reset image size

This fact allows us to apply the approximation

$$\begin{equation} p_\ell(x)\approx p_\ell(0)=\sqrt{2\mu}\,\sqrt{E_\ell-V(0)}. \end{equation} \tag{ C.10 }$$

Because the numerical treatment of (3) gave V (0) = −0.488 ℏω₀ [8], we substitute

$$\begin{equation} V(0)/(\hbar\omega_0)\approx -1/2 \end{equation} \tag{ C.11 }$$

into (C.10) and with the help of the relation

$$\begin{equation} \frac{1}{\hbar}\int_x^{x_\ell} \mathrm{d}x\, p_\ell(x) = \frac{1}{\hbar}\int_0^{x_\ell} \mathrm{d}x\, p_\ell(x)-\frac{1}{\hbar}\int_0^x \mathrm{d}x\, p_\ell(x), \end{equation} \tag{ C.12 }$$

which reduces, with the help of the quantization rule (C.4), to

$$\begin{equation} \frac{1}{\hbar}\int_x^{x_\ell} \mathrm{d}x\, p_\ell(x)\approx\frac{\pi}{2}\left(\ell+\frac{1}{2}\right)-\frac{1}{\hbar}\, p_\ell(0)\,x, \end{equation} \tag{ C.13 }$$

we arrive at the simplified WKB wave function

$$\begin{equation} \fl \varphi^{\rm{(WKB)}}_\ell(x)\approx \frac{\sqrt{2\, \alpha}\,}{(2\ln (\ell+1)+1)^{1/4}{\sqrt{\pi(\ell+1)}}} \cos\left(\sqrt{2\ln (\ell+1) +1}\,\alpha x-\frac{\ell\pi}{2}\right). \end{equation} \tag{ C.14 }$$

We note that the wave functions with odd index ℓ vanish at the origin.

Figure C.1 compares this approximation to the numerically obtained wave functions for parameters arising in our example of section 4.3.1. We emphasize that indeed in the neighbourhood of the origin, the two wave functions are in remarkable agreement.

In order to find an approximate analytical expression for the matrix element (A.16), we fit a normalized Gaussian

$$\begin{equation} G(x)\equiv\left( \frac{1}{\pi \sigma^2}\right)^{1/4} \mathrm{e}^{-x^2/2\sigma^2} \end{equation} \tag{ C.15 }$$

to the numerically obtained [8] ground state wave function φ₀(x) and find that

$$\begin{equation} \sigma =1.22\,\alpha^{-1} \approx \sqrt{3/2}\,\alpha^{-1}. \end{equation} \tag{ C.16 }$$

Again in figure C.1 we compare the numerically exact expression and the approximation and find that there is excellent agreement.

With the approximations (C.14) and (C.15) the matrix element (C.9) for the transition |00〉 → |p − 1,q − 1〉 takes the form

$$\begin{equation} \kappa_{00,p-1\,\,q-1} \approx \int^{\infty}_{-\infty}\mathrm{d}x\,G(x)^2\,\varphi_{p-1}^{\rm{(WKB)}}(x)\varphi_{q-1}^{\rm{(WKB)}}(x) \end{equation} \tag{ C.17 }$$

and becomes a standard integral.

Indeed, with the abbreviation

$$\begin{equation} w(p)\equiv(2\ln p+1)^{1/2} \end{equation} \tag{ C.18 }$$

and for p and q both even, or both odd, we find the expression

$$\begin{eqnarray} \kappa_{00,p-1\,\,q-1} & \approx & \frac{\alpha}{\pi(p q)^{1/2}} \frac{{\mathrm{e}}^{-3/8(w(p)-w(q))^2 }+(-1)^{p+1}\,\mathrm{e}^{-3/8(w(p)+w(q))^2 }}{(w(p) w(q))^{1/2}} \nonumber\\ \end{eqnarray} \tag{ C.19 }$$

which reduces with the help of the relation

$$\begin{eqnarray*} &&w(p)^2+w(q)^2=2\ln(p q)+2=2\,\ln N + 2 \end{eqnarray*}$$

to

$$\begin{eqnarray*} &&\kappa_{00,p-1\,\,q-1} = \frac{\alpha}{\pi }\,\,\frac{\mathrm{e}^{3/4 (w(p) w(q)-1)} +(-1)^{p+1}\,\mathrm{e}^{-3/4(w(p)\cdot w(q)+1) }}{N^{5/4}(w(p) w(q))^{1/2}} \end{eqnarray*}$$

otherwise it vanishes. Hence, the last line shows that the matrix element depends only on the product

$$\begin{equation} w(p) w(q)= \left( (\ln N+1)^2-(\ln (p/q))^2\right)^{1/2}. \end{equation} \tag{ C.20 }$$

We assume that

$$\begin{equation} |\ln(p/q) |\ll \ln N \end{equation} \tag{ C.21 }$$

and approximate (C.20) by

$$\begin{eqnarray*} &&w(p) w(q) \approx \ln N +1. \end{eqnarray*}$$

As a consequence, we arrive at the asymptotic expression

$$\begin{equation} \kappa_{00,p-1\,\,q-1} \approx \frac{\alpha}{\pi}\,\frac{1}{ (N\ln N)^{1/2}}\propto N^{-1/2} \end{equation} \tag{ C.22 }$$

no matter whether p and q are both even or both odd.

For N odd the matrix element for the transition |0,0〉 → |N − 1,0〉 given by

$$\begin{equation} \kappa_{00,0\,\,N-1} \approx \int^{\infty}_{-\infty}\mathrm{d}x\,G(x)^3\,\varphi_{N-1}^{\rm{(WKB)}}(x) \end{equation} \tag{ C.23 }$$

yields after integration the asymptotic behaviour

$$\begin{equation} \kappa_{00,0\,\,N-1} \approx\frac{2\alpha}{3\pi}\left(\frac{6\pi}{e}\right)^{1/4} \frac{1}{N\,w(N)^{1/2}} \propto N^{-1} \end{equation} \tag{ C.24 }$$

while for N even it vanishes.

For the matrix element of the transition κ_{p−1,q−1 P−1,Q−1} with PQ = N² it is sufficient here to give a crude estimate which contains four WKB wave functions. With the help of (C.14) it follows that

$$\begin{equation} \left|\kappa_{p-1,q-1\,\,P-1,Q-1}\right| < \frac{\mbox{const.}}{N}. \end{equation} \tag{ C.25 }$$

Again we note that we are considering distinguishable particles here. In the bosonic case the matrix elements must be multiplied by the appropriate factors given by (9)–(11).

Finally, we use the example 7 × 181 = 1267 to illustrate how much our factorization protocol suffers if the factors p and q deviate considerably from the assumption (C.21). Here the approximation (C.22) gives a Rabi frequency which is 70% larger than the frequency (C.19) calculated from factors 7 and 181. As a consequence, the interaction is switched off at a time T_opt given by (22), which is too early and the probability of finding the factor state |7 − 1,181 − 1〉 is only 53%. Three runs are therefore needed to find the factors with a probability of about 90%.