Quantum advantages for transportation tasks: projectiles, rockets and quantum backflow

Consider a scenario where a quantum particle is initially prepared in some bounded region of space and left to propagate freely. After some time, we verify if the particle has reached some distant target region. We find that there exist"ultrafast"("ultraslow") quantum states, whose probability of arrival is greater (smaller) than that of any classical particle prepared in the same region with the same momentum distribution. For both projectiles and rockets, we prove that the quantum advantage, quantified by the difference between the quantum and optimal classical arrival probabilities, is limited by the Bracken-Melloy constant $c_{bm}$, originally introduced to study the phenomenon of quantum backflow. In this regard, we substantiate the $29$-year-old conjecture that $c_{bm}\approx 0.038$ by proving the bounds $0.0315\leq c_{bm}\leq 0.072$. Finally, we show that, in a modified projectile scenario where the initial position distribution of the particle is also fixed, the quantum advantage can reach $0.1262$.


I. INTRODUCTION
Much of current research in quantum theory focuses on the exploitation of quantum effects in communication and computation.Nevertheless, quantum systems are originally found to be advantageous for mechanical tasks.A paradigmatic example is the tunneling effect [1]: A quantum particle can be detected in regions of space that are classically forbidden by energy considerations.Another noteworthy example is quantum backflow: A free quantum particle with positive momentum can be observed to propagate backwards.Quantum backflow was first identified by Allcock in the context of the timeof-arrival problem [2], and later isolated by Bracken and Melloy [3].More recent examples of quantum advantage in mechanical systems can be found in [4] and [5].
The advantages that quantum mechanical systems might offer for transportation, understood as the quick dispatch of massive particles through free space, are, however, unexplored.Some effort has been paid to investigate the properties of a hypothetical quantum time-ofarrival operator [6] in connection with quantum backflow.Perhaps due to its foundational character, this research program has not produced so far any concrete task where quantum mechanical systems have the upper hand.
In this work, we prove the advantage of quantum mechanical systems over their classical counterparts in a practical transportation task, which we call the projectile scenario.Consider a situation where a non-relativistic one-dimensional quantum particle (a projectile) is prepared in some bounded region of space B and left to propagate freely.After some time ∆T , we measure if the particle is in some distant target region R.For a fixed initial quantum state ρ with spatial support in B, we compare the probability of detection in R with that of a classical particle, initially prepared in B with the same momentum distribution as ρ.
We find that there exist what one might call ultra-fast states (ultra-slow states), whose probability of detection in R at time ∆T is strictly greater (smaller) than that of any classical particle.A natural figure of merit for quantum advantage in the ultra-fast regime is the difference between the quantum and the maximum classical probabilities of arrival.Likewise, in the ultra-slow regime one can consider the difference between the minimum classical and the quantum probabilities of arrival.We find that the maximum quantum advantage in either case does not depend on the distance between the preparation and target regions, but only on the parameter α := M |B| 2 /∆T .For finite values of α, the maximum quantum-classical gap can be computed up to precision δ by diagonalizing an N × N matrix, with N = O (log (1/δ)).
We prove that the maximum quantum advantage, achieved in the limit α → ∞, equals the Bracken-Melloy constant [3], which was numerically estimated to have the value c bm ≈ 0.0384517 [7,8].This conjectured value was, however, not computed with any rigorous error bounds.In fact, until now there was no reason to believe that c bm was smaller than 1.In this regard, we argue that 0.0315 ≤ c bm ≤ 0.0725, hence providing the first upper bound on c bm .
As we show, the appearance of c bm is not a coincidence: through simple metaplectic transformations we connect the quantum projectile problem with a variety of scenarios related to and generalizing quantum backflow, including quantum backflow itself.All such effects are therefore manifestations of the same mathematical phenomenon, seen through different coordinate systems.In the light of the recent interest in experimentally demonstrating quantum backflow [9][10][11][12][13], we argue that projectile scenarios are more experimentally friendly and operationally interesting.
To arrive at a transportation task with a quantum advantage beyond the Bracken-Melloy constant, we consider a scenario in which several projectiles are sequentially released, namely, a quantum rocket.However, it turns out that c bm also limits the advantage of a quantum rocket over a classical analog with the same lift-off zone, combustion chamber size and rocket and fuel mo-mentum distributions.
Nevertheless, we show that a superior quantum advantage can actually be attained in a variant of the projectile scenario where the quantum projectile is compared with a classical particle having the same position and momentum distributions.
The paper is structured as follows: in section II we introduce and solve the projectile scenario; the connection between quantum projectiles and other examples of quantum advantage in mechanical systems is explained in section III.In section IV we provide a simple model for quantum rockets and use it to prove that the classicalquantum gap in such artifacts is also limited by the Bracken-Melloy constant.In section V, we add a natural constraint to the projectile scenario so that the Bracken-Melloy limit can be superseded.Finally, in section VI we present our conclusions.We also provide some Appendices in which the lengthier computations are made more explicit.

II. CLASSICAL VS. QUANTUM PROJECTILES
Our starting point is a classical projectile of mass M , prepared at time t = 0 in the region [0, L].At time t = ∆T > 0, we observe whether the projectile has reached region [a, ∞), with a > L (see Figure 1).If we ignore where exactly in [0, L] the projectile was prepared, then the probability of finding it in [a, ∞) at time ∆T is, at most, Prob (p ≥ M (a − L)/∆T ), where p denotes the projectile's linear momentum.This corresponds to a configuration where the projectile was prepared at x = L at time t = 0. Similarly, the probability to find the projectile in [a, ∞) at time ∆T is, at least, Prob (p ≥ M a/∆T ), which corresponds to an initial preparation at x = 0. Now, let us assume that the projectile is, in fact, a quantum mechanical system.Let S(R) denote the set of quantum states with spatial support in R ⊂ R. We will omit the parentheses whenever R is an interval, and thus denote by ρ ∈ S[0, L] the initial quantum state of the projectile.While the projectile is freely propagating, its dynamics are governed by the kinetic Hamiltonian H = P 2 /2M , where P denotes the projectile's linear momentum operator.The probability to find the quantum projectile in region [a, ∞) after time ∆T can be found by simple application of the Born rule: it is tr U ρU † Θ(X − a) , where U := e −iH∆T and Θ is the Heaviside step function.Note that we work in units where ℏ = 1.
If, after time ∆T , the quantum projectile is found in [a, ∞) with probability greater than any classical particle initially prepared in [0, L] with the same momentum distribution, we say that the quantum projectile is ultra-fast.If, on the contrary, the projectile is detected with probability lower than the classical minimum, we say that the projectile is ultra-slow.To gauge how ultra-fast or ultra-slow a quantum projectile in state ρ is, we consider the difference between the quantum and optimal classical probabilities of arrival.
Let us deal with the ultrafast case first.As we saw in the first paragraph of this section, a classical projectile with momentum distribution ν(p)dp will be detected in [a, ∞) at time ∆T with probability at most Prob (p ≥ M (a − L)/∆T ).The probability of this event is to be evaluated on the distribution ν(p)dp.Since we have assumed ν(p)dp to be the same as the momentum distribution of a quantum particle in state ρ, this implies Thus the quantum advantage, if it exists, is given by tr(ρΩ F (M, a, ∆T )), with where, in the first term of the right-hand side, we made use of the identity [14] U † XU = X + ∆T P/M .We wish to find the largest advantage achievable with a quantum projectile.That is, we are interested in the quantity We next exploit this observation to prove that φ F is just a function of α := M L 2 /∆T .In particular, φ F does not depend on a, the location of the target region: remarkably, quantum projectiles are equally advantageous no matter how large the flight distance.
Let σ : R 2 → R 2 be an affine linear transformation and consider the vector of operators (X, P ).If σ is metaplectic, namely [σ(X, P ) 1 , σ(X, P ) 2 ] = [X, P ] = i, then, as we show in Appendix A, there exists a unitary U σ such that Now, consider the unitary V associated to the metaplectic map For α = M L 2 /∆T , it follows that V Ω F (M, a, ∆T )V † = Θ(X + P ) − Θ(P ) =: Ω, Hence, φ F is just a function of α.We call the righthand side of the above equation the standard projectile problem, or standard problem for short.Note that the standard problem corresponds to determining the maximum quantum advantage of an ultrafast projectile of mass M = 1, prepared in the region [− √ α, 0], to be found in region [0, ∞) after time ∆T = 1.
So far we have only considered ultrafast projectiles.For the ultraslow case, the story is pretty much the same, but opposite: namely, we are now interested in not finding the particle in the target region [a, ∞) after time ∆T has elapsed.The optimal classical strategy is now to concentrate all the mass at point x = 0.In this case, the probability that a classical projectile, prepared at time t = 0 in [0, L] with the same momentum distribution as the quantum state ρ, reaches the target region at time t = ∆T is given by and so the quantum advantage, if it exists, of not finding the particle in the target region is quantified by tr(ρΩ S (M, a, ∆T )), with The maximum quantum advantage is thus tr(ρΩ S (M, a, ∆T )).
As it turns out, φ S = φ, and so the functions φ F , φ S are identical.Indeed, consider the transformation Since [σ(X, P ) 1 , σ(X, P ) 2 ] = −i, this map does not define a unitary transformation over the set of quantum states.Rather, it defines an anti-unitary transformation U σ , as explained in Appendix A. Now, the argument above relating linear optimizations over subsets of quantum states also extends to anti-unitary transformations.The reader can verify that, applying U σ to the standard problem with α = M L 2 /∆T , one ends up with the definition of φ S , and, therefore, φ S (M, L, a, ∆T ) = φ M L 2 /∆T .
In section II A, we will prove that φ(α) > 0 for all α > 0, i.e., there exist ultrafast and ultraslow quantum states in any projectile scenario.From eq. (3) it is clear that φ(α) is a non-decreasing function.Moreover, as shown in section III, its limiting (supremum) value φ(∞) corresponds to the Bracken-Melloy constant c bm [15], conjectured to have the value 0.0384517.We conclude that quantum projectiles can exhibit a limited advantage with respect to their classical counterparts.
We finish this section by introducing yet another projectile scenario.As before, we wish the quantum projectile to have a larger probability of arrival, but this time we award some advantage to the classical projectile: namely, we compare the probability to detect the quantum projectile in the region [a, ∞) with the maximum probability of detecting the classical one in the larger region [a − b, ∞), with b > 0. This problem can be reduced, via the transformation (2), to an optimization of ⟨Θ(X + P ) − Θ(P + β)⟩ ρ over ρ ∈ S[− √ α, 0], with α = M L 2 /∆T , β = b M/∆T .We denote this problem the extended standard problem, with solution φ(α, β).Clearly, φ(α, β) is non-increasing in β and φ(α, 0) = φ(α).Obviously, lim β→∞ φ(α, β) = 0, and so one cannot reduce the extended standard problem to the standard problem.From the formulation of the standard problem (3), one can immediately deduce that φ is a non-decreasing function of α ∈ [0, ∞), with φ(0) = 0 and φ(α) ≤ 1.It remains to see that φ(α) ̸ = 0 for some α.To do this, we need to study the spectrum of Ω := Θ(X + P ) − Θ(P ) restricted to the space S[− √ α, 0].In Appendix B we prove that, in position representation, Let K(x, y) be the kernel of this integral operator.If α > 0, then we can choose z ∈ (− √ α, 0) such that K(0, z) = K(z, 0) * ̸ = 0. Since K(0, 0) = 0, by the determinant criterion it follows that the 2 × 2 matrix {K(x, y)} x,y=0,z is not negative semidefinite.In particular, it has a positive eigenvalue λ, with eigenvector (c 0 , c z ) T .Now, consider the ket where χ C denotes the characteristic function of C ⊂ R. For small enough ε, We conclude that φ(α) > 0 for all α > 0, so ultrafast and ultraslow states exist in all projectile scenarios.
The problem of computing φ(α) for different values of α is more convoluted.Note that the kernel K(x, y) is analytic in x, y; hence, for x, y ∈ [− √ α, 0], we can approximate it up to arbitrary precision by a polynomial on x and y of sufficiently high degree.When we replace K(x, y) by its N th order Taylor expansion, we arrive at a new operator Ω N , which can be shown to be close in operator norm to Ω, restricted to the subspace of wave functions defined in [− √ α, 0].In turn, Ω N only has support on the finite-dimensional subspace spanned by vectors of the form [− √ α,0] dxx k |x⟩, where k runs from 0 to the degree in x of the kernel of Ω N .Hence Ω N can be exactly diagonalized.In Appendix B this argument is developed to conclude that, for finite α, we can compute φ(α) to any precision δ we want by diagonalizing a matrix of size N ≲ max(α, log(1/δ)).In the same appendix, the reader can also find the following (tight) linear upper bound for φ(α): The function φ(α) is plotted for α ∈ [0, 100] in Figure 2. As it can be appreciated, φ(α) roughly looks like a concave function, but not quite: at regular intervals, the slope of the function becomes very small.Such 'steps' seem to decrease in amplitude as α grows, and, actually, for α ≫ 1, the function appears to be well approximated by the ansatz r + sα −1/2 .To grasp the maximum quantum advantage, we need to study the limiting case α = ∞.The problem thus consists in determining the spectrum of Ω, restricted to the space L 2 (−∞, 0].To study this case, it is convenient to switch to the Wigner function representation.
The Wigner function of a quantum state ρ is For convenience, we recall the properties of Wigner functions in Appendix A. The most important one for us is the fact that Wigner functions behave nicely under metaplectic transformations in phase space.Namely, for any metaplectic transformation U σ , it holds that Furthermore, for any bounded measurable function f : R → R and a, b, c ∈ R, we have that tr(ρf (aX where some care has to go into the precise meaning of the integral whenever the integrand is not Lebesgue integrable.Finally, note that, if ρ has a convex support R in either position or momentum, then the support of its Wigner function W ρ (x, p) corresponding to that variable is also contained in R.
Now, for any state ρ, we have, by eq. ( 7), that The last factor on the integrand will vanish everywhere, except in the regions Λ + = {x + p ≥ 0, p ≤ 0}, where it equals 1, and Λ − = {x + p ≤ 0, p ≥ 0}, where it equals −1.However, if ρ ∈ S(−∞, 0], then W ρ (x, p) = 0, for x > 0. Since (x, p) ∈ Λ + implies x ≥ 0, it follows that the first region does not contribute to the integration above.Hence, The problem of integrating Wigner functions over wedges (without any further constraints) was studied by Werner [16] in the context of time-of-arrival operators.The idea is that all wedges can be taken to each other via a metaplectic transformation, and therefore it suffices to study the wedge [0, where we have used that S(−∞, 0] is the space of states that satisfy the condition tr(ρΘ(−X)) = 1.Werner considers the operator corresponding to integrating Wigner functions over the quadrant x, p ≥ 0, and determines its spectrum to be [−0.155940,1.007678].Therefore, φ(∞) ≤ 0.155940.This bound, however, does not take into consideration the constraint tr(ρΘ(X + P )) = 1.
To account for it, we add to Werner's operator a linear combination of operators corresponding to integrating Wigner functions over hyperbolic regions in the quadrant x, p ≤ 0. Since our Wigner functions vanish in that quadrant, the infimum of the spectrum of the new operator (which can also be determined with the techniques in [16]) also provides an upper bound for φ(∞).We numerically find the bound φ(∞) ≤ 0.0725, see Appendix D.
In addition, via variational methods, we show that φ(∞) ≥ 0.0315.This figure is obtained by optimizing linear combinations of the average values of the operators Ω, Θ(X) over density matrices with support on the first N + 1 number states {|n⟩ : n = 0, ..., N }, i.e., (X + iP ) |n⟩ = √ 2n |n − 1⟩, see Appendix C for details.A plot of the Wigner function of a quantum state approximately in S(−∞, 0] and approximately achieving this value can be found in Figure 3 (left).
In the next section, we will show that φ(∞) = c bm , the Bracken-Melloy constant [3], which is conjectured to have the value 0.0384517 [7,8].Our bounds 0.0315 ≤ c bm ≤ 0.0725 therefore support this widespread belief.

III. CONNECTION WITH OTHER QUANTUM MECHANICAL EFFECTS
As we have seen, the ultrafast (ultraslow) projectile problem is equivalent to the standard problem, since a unitary (anti-unitary) transformation takes us from the latter to the former.We next see that the standard projectile problem is similarly connected to the most extreme manifestation of other quantum mechanical effects.The exact correspondences are summarized in Table I.The question of understanding the relation between some of these effects was raised in [17] and partially answered in [18].Our results answer the challenge posed in [17] from a different point of view, namely, that of (anti-)unitary equivalence, and extend the connection to other mechanical effects.
Let us start with the phenomenon of quantum backflow [2,3,8,19].Consider a pure state that only has positive momentum and that is evolving freely.In position representation, we can write it as for some function ϕ such that The probability flux at the origin is therefore and thus the integrated flux at the origin from time 0 to time ∆T is Note the similarity with eq. ( 5).Guided by classical intuition, one would expect this integrated flux to be nonnegative, since the particle is only moving to the right.However, for some quantum states ϕ(x, t), this magnitude can be negative: in that case, we speak of quantum backflow.
Alternatively, we can interpret quantum backflow as a decrease in the probability of detecting a particle with positive momentum in the region [0, ∞).This is so because, by the continuity equation the integrated flux satisfies: where |ψ⟩ = dxψ(x, 0) |x⟩ and U = e −i P 2 2M ∆T .Call P[0, ∞) the space of all states with positive momentum support.From all the above it follows that the maximum amount of backflow is given by where we used the identity Θ(z) = 1 − Θ(−z).The number c bm , known in the literature as the Bracken-Melloy constant [3], is thus the solution a problem of the form sup ρ∈S tr(ρA), for some space of states S and some operator A. In fact, this problem can be obtained from the standard problem with α = ∞ via the anti-metaplectic Most of the optimization problems considered in this paper are of the form maxρ∈S tr(ρΩ), for some operator Ω and some set of states S.This table contains the definitions of each problem and the reversible transformations mapping the standard problem to any other.S(R) denotes the set of states with position support in R ⊂ R, and P(R) denotes the set of states with momentum support in R ⊂ R. We use the shorthand σ(x) := σ(x, p)1 and σ(p) := σ(x, p)2, and omit parentheses whenever R is an interval.
Going through the literature on quantum backflow, one finds that c bm is conjectured to have the value 0.0384517 [7,8].A figure of 0.038452 is obtained in [8] by fitting many points of (an approximation to) the graph of φ(α) with the ansatz r − sα −1/2 and, a figure of 0.0384517 is obtained in [7], by fitting such points to a degree 3 polynomial over α −1/2 .To our knowledge, prior to our work there were no rigorous, non-trivial upper bounds on c bm , and the best lower bound fell 41% short of the conjectured value of the constant [20].Our results in the preceding section hence give mathematical support to the conjecture c bm ≈ 0.0384517.
In Table II we present another set of quantum effects that are mathematically equivalent, not to the standard problem, but to the extended standard problem with α = ∞, which we express, via the transformation σ(x, p) = (x − β, p + β), as an optimization of Ω over the set of states S(−∞, β]. One of these effects is a variant of quantum backflow in which the particle evolves in the presence of a constant force [21].That is, with the Hamiltonian given by H = P 2 /2M − F X. In [17] Goussev proves that this effect is at the same time equivalent to something he calls quantum reentry.Quantum reentry is an effect that consists in preparing a particle in S(−∞, 0], letting it evolve and then measuring a negative probability flow in some point l ≥ 0. That is, the quantity under consideration is − t2 t1 dtj(l, t) for some t 2 > t 1 > 0, which can again be easily transformed to the semi-infinite standard problem, as also shown in Table II.In particular, the maximum probability transfer in both these effects is the same.
Finally, we note that the extended standard problem is equivalent to computing the maximum expression of quantum backflow when the initial momentum is in the region [−γ, ∞) for some γ ∈ R, as shown in Table II.Thus, when the initial momentum is in this region, the probability "backflow" acts as if there were a constant force acting on the system, since these two problems are again equivalent.This seems to have gone unnoticed by Bracken, who studied the former effect in [15], despite having studied the latter in [21] together with Melloy.

IV. CLASSICAL VS. QUANTUM ROCKETS
The low value of c bm constitutes a severe obstruction to any practical application of quantum systems for transportation tasks.How to overcome this limit?A tempting idea is to consider scenarios where a transiting quantum projectile launches a second quantum projectile.Iterating this procedure, we arrive at the notion of a quantum rocket, i.e., a quantum mechanical system that, from time to time, throws away some fuel mass in the direction opposite to the intended motion.Since this rocket scenario encompasses the quantum projectile scenario, its maximum quantum advantage is lower-bounded by the Bracken-Melloy constant.Furthermore, one would imagine that, should we prepare the fuel in the right quantum state, the limited quantum advantage present in quantum projectiles could be somehow bootstrapped, hence increasing the overall advantage of the quantum rocket with respect to a classical rocket whose fuel combustion has an identical momentum distribution.
Unfortunately, this is not the case, at least for a large class of quantum rockets.Consider a minimal model for a quantum rocket, where, at time t, the rocket itself is regarded as a 1-dimensional particle of mass M (t) and zero spin.The state of the rocket at time t is therefore specified through a trace-class positive semidefinite operator ρ(t) : L 2 (R) → L 2 (R).For most of its flight, the rocket will be propagated by the kinetic Hamilto- TABLE II.Some of the problems which are (anti-)metaplectically equivalent to the semi-infinite standard problem, with the same notation as in Table I.In the last row, the normalization factor of the metaplectic transformation is C := (t2 − t1)/t2.nian H = P 2 R /2M (t).At times 0 = t 1 < t 2 < ... < t N , though, the rocket's free evolution is interrupted: namely, at time t j the rocket burns and releases a predetermined amount of fuel m j instantaneously, thus decreasing its overall mass by the same amount.
To model the instantaneous combustion of fuel of mass m < M , we consider a completely positive tracepreserving (CPTP) map Υ that, acting on the rocket's state ρ(t), returns a density matrix representing the joint state of the fuel F and that of the rest of the rocket R, whose mass is now M − m, see Figure 4.
Call X F , P F (X R , P R ) the absolute position and momentum operators of the fuel (the rest of the rocket), and let X CM , P CM (X REL , P REL ) denote the canonical variables of the center of mass (the relative coordinates between systems F and R), with: Let U M,m be the (symplectic) unitary that switches between the R, F and CM, REL representations and define ω CM,REL ≡ U M,m Υ(ρ)U † M,m .Since Υ is an internal and instantaneous operation, it cannot modify the rocket's center of mass degree of freedom.This means that tr REL (ω) = ρ.For ρ = |ψ⟩⟨ψ|, this last relation implies that ω = |ψ⟩⟨ψ| ⊗ σ ψ , for some quantum state σ ψ .
It must be noted that σ should have been prepared in the rocket's combustion chamber.If we assume that the combustion chamber is centered in the rocket's center of mass and has length λ, then σ must have spatial support in [−λ/2, λ/2].
In describing the overall flight of the rocket, we assume that, at time t j , the quantum rocket, with mass M j , will release a mass m j of fuel in state σ j ∈ S[−λ/2, λ/2] (in the relative frame of reference).Hence, the mass and state of the rocket will be instantaneously updated to M j+1 = M j − m j , ρ → tr F (Υ(ρ; σ j , M j , m j )).
We consider the probability to find the rocket at time t N +1 > t N in the region [a, ∞).This is to be compared with the maximum probability that an analog classical rocket arrives at the same region in time t N +1 .Like in the projectile scenario, this classical rocket is assumed to have, at time t 1 , the same initial mass, initial momentum distribution and initial spatial support as the quantum one.At time t j , this classical rocket will burn a mass m j of fuel, and the phase space distribution of the classical fuel in the fuel's reference frame relative to the rocket is demanded to have the same momentum distribution and spatial support as σ j .
In these conditions, in Appendix E we show that the difference between the quantum and classical arrival probabilities is also limited by c bm .This no-go result crucially relies on eq. ( 9), which expresses the assumption that the fuel's interaction with the rocket is instantaneous.Physically, this corresponds to a configuration where the combustion chamber is open on both sides, i.e., the fuel is allowed to exit the rocket, not only against the rocket's direction of motion, but also towards it.Assumption eq. ( 9) allows us to map the computation of the rocket's maximum quantum advantage to the standard problem (with further state constraints) through a metaplectic transformation.

V. A TRANSPORTATION SCENARIO WITH A QUANTUM ADVANTAGE THAT SUPERSEDES THE BRACKEN-MELLOY CONSTANT
In view of the last result, it would be reasonable not to expect significant gaps between the arrival probabilities of quantum and classical particles.As it turns out, though, a simple variation of the way we compare classical and quantum projectiles is enough to find quantum advantages for transportation way beyond the Bracken-Melloy constant.Note that there exist known variations of the quantum backflow problem that achieve quantum advantages greater than the limit set by Bracken and Melloy [18,[22][23][24].Those effects are, however, unrelated to transportation tasks.
In Section II, we compared the behavior of a quantum projectile (or a rocket) with respect to that of a classical one with the same momentum distribution and the same spatial support at time t = 0. Could the quantum advantage be amplified if we demanded further constraints on the initial position distribution µ(x)dx of the classical projectile, besides its support?In the extreme case, we could demand µ(x)dx to coincide with the position distribution of the quantum projectile.
Consider thus the following problem: let ρ denote the density matrix of a particle of mass M , and let µ(x)dx, ν(p)dp be its position and momentum distributions at time t = 0.As before, we let the projectile evolve freely for time ∆T and then check whether the projectile is in [a, ∞); call p q (ρ) the corresponding probability.How much does p q (ρ) differ from the maximum arrival probability of an analog classical particle, with initial position and momentum distributions µ(x)dx, ν(p)dp?
The maximum classical probability of arrival is where W (x, p) represents the probability distribution of the classical particle in phase space at time t = 0.The maximum quantum-to-classical advantage in this projectile scenario is therefore Φ ⋆ = sup ρ∈S W(ρ), where W(ρ) := p q (ρ) − p ⋆ c (ρ).This is a nested max-min optimization problem, whose solution can be proven independent of a, M, ∆T [25].
In Appendix F, p ⋆ c (ρ) is shown to equal s(∞), the solution of the system of ordinary differential equations with initial conditions s(−∞) = q(−∞) = 0.Here Θ + (z) is meant to be 1 for z > 0 and 0 otherwise.s(∞) can be computed numerically via, e.g., Euler's explicit method.Since we know how to compute p ⋆ c (ρ), one could, in principle, use gradient ascent methods to find the maximum of W(ρ), over all quantum states with, say, support on the space spanned by the first N number basis vectors.That is, we could parametrize any such state ρ as ρ = N m,n=0 ρ m,n |m⟩ ⟨n| and then follow the gradient of W(ρ) with respect to the variables ρ m,n .Unfortunately, W is a concave function, so the method is not guaranteed to converge to the absolute maximum.Moreover, we empirically observe that, starting from a random state, projected gradient methods typically converge to very suboptimal To find a suitable starting point for gradient ascent, we considered the following approach: suppose that there existed a linear operator Z such that p ⋆ c (ρ) ≤ tr(Zρ), (12) for all states ρ.Then we could maximize the value over all density matrices with support on the first N number states.The result would provide us with a lower bound on Φ ⋆ .In addition, if the maximizer ρ ⋆ satisfied W Z (ρ ⋆ ) > 0, then that state would be a good starting point for gradient ascent.Now, how to identify an operator Z satisfying (12)?Let f, g : R → R be two functions such that for all x, p.Then, for any distribution W (x, p) in phase space with marginals µ(x), ν(p), It follows that the operator Z = f (X) + g(P ) fulfills condition (12).In fact, the dual of problem (10) is the maximum of the right-hand side of eq. ( 15) over all such functions f, g.Take M = ∆T = 1, a = 0. We observe that the functions f = g = Θ satisfy (14), and hence, the supremum of the spectrum of the operator Ω = Θ(X + P ) − Θ(P ) − Θ(X) provides us with a lower bound for Φ ⋆ , as tr(ρΩ) ≥ Φ * .
If we truncate this operator in the number basis, we are looking at the maximum eigenvalue of the matrix (M (N ) nm : n, m = 0, ..., N ), with For N = 170, the maximum eigenvalue of this matrix is 0.1113: the reader can find a plot of the Wigner function of the corresponding eigenvector in Figure 3 (right).Taking N = 1700, we obtain the tighter bound Φ ⋆ ≥ 0.1228.The maximum quantum advantage in this projectile scenario is therefore substantially greater than the conjectured value of c bm , or even its upper bound 0.0725, derived in section II A.
Applying gradient methods on those states to improve their W value proved to be tricky, though.Call ρ ⋆ the state corresponding to the eigenvector of M (N ) .We observe that, even for low values of N (say, N = 30), we need to use a very small step size in eq. ( 11) to estimate p ⋆ c (ρ ⋆ ) precisely.When we do so, we find that p ⋆ c (ρ ⋆ ) ≈ tr{ρ ⋆ (Θ(X)+Θ(P ))}: that is, for such quantum states, our upper bound (12) on p ⋆ c is (approximately) tight.Around the eigenvectors of M (N ) , the gradient of W explodes, possibly because the function is not everywhere differentiable.Using random perturbations of ρ ⋆ as a seed, projected gradient methods only produced states with a objective value slightly smaller than W(ρ ⋆ ).From all the above, it is thus natural to conjecture that the obtained value of 0.1228 is (close to) a local maximum of W, at least among quantum states with support in {|n⟩ : n = 0, ..., 1700}.
On the other hand, note that after a suitable metaplectic transformation the problem sup ρ tr(ρΩ) becomes sup ρ tr ρ Ω , where with X k := cos(2πk/3)X + sin(2πk/3)P .The operator 2 k=0 sgn(X k )/3 is the one studied by Tsirelson in [4].The best known bounds for its spectrum are given in [5].Using Equation (D20) in [5], one obtains that Φ * ≥ −0.5 + 1.5 × √ 0.17491 = 0.1262.In particular, this shows how unreliable the numerical estimation of these quantities is, even after using a basis with 1700 number states, and thus the importance of getting good upper bounds as well as lower bounds.

VI. CONCLUSION
In this letter, we have investigated how the dynamics of quantum and classical projectiles differ, using the probability of arrival at a distant region of space as a figure of merit.We found that non-relativistic quantum particles can arrive at a distant region with higher or lower probability than any classical particle with the same initial spatial support and momentum distribution.Curiously enough, the maximum gap between quantum and classical probabilities is independent of the distance to the arrival region, and just depends on the mass M and spatial support L of the projectile and its flying time ∆T through the single parameter α = M L 2 /∆T .
The discrepancy between the quantum and classical arrival probabilities is, however, limited by the Bracken-Melloy constant c bm ≈ 0.0384517.As we showed, the maximum quantum advantage of rockets with an open combustion chamber is also bounded by this value.Our no-go result does not apply, however, to rockets with a 1-side closed combustion chamber, which just allows the fuel to exit the rocket opposite to its direction of motion.Whether such rocket models are also limited by c bm , or on the contrary, they can achieve arrival probabilities much higher than classical is an interesting topic for future research.
In a similar direction, we showed that considerable quantum-classical gaps of at least 0.1262 can be observed if we demand classical projectiles to reproduce the initial position distribution of the quantum projectile.It is an open problem whether this figure is indeed close to the maximum quantum advantage, and whether this effect can be exploited for real transportation tasks.

VII. APPENDIX
In this Appendix we perform some of the lengthy computations that give the claims of the main text.It is organized as follows.In section A, we review the relevant properties of the Wigner function.In section B we give an explicit formula for φ(α) and show how to numerically approximate it for finite α.In section C we compute a numerical lower bound for φ(∞).In section D we compute a numerical upper bound for φ(∞).Section E proves that there is no advantage in our model of the quantum rocket with respect to a single projectile.Finally, in section F we describe the numerical methods used in the restricted projectile scenario.

Appendix A: Notes on the Wigner function
The Wigner function of a quantum state ρ is This is a partial Fourier transform on the function ⟨x − y/2| ρ |x + y/2⟩, and as such is defined with the usual density arguments from the states ρ such that ⟨x| ρ |y⟩ is a Schwarz function of two variables.We now study the action of metaplectic transformations on the Wigner function.Suppose then that σ : R 2 → R 2 is affine-linear, and [σ(X, P ) 1 , σ(X, P ) 2 ] = [X, P ].Then, calling σ the linear part of σ, we must conclude that σ ∈ SL 2 (R).It is well known that SL 2 (R) = Sp 2 (R), which gives rise to the name metaplectic that we have used in the main text.Furthermore, the KAN decomposition of SL 2 (R) is That is, every matrix decomposes as a product of a rotation, a dilation and a translation.These all correspond to time evolutions of quadratic Hamiltonians (P 2 + X 2 , XP +P X and P 2 or X 2 , respectively).Finally, the affine part of the map can be realized by time-evolving with the Hamiltonians X and P .These six Hamiltonians thus give rise to the unitaries U σ mentioned in the main text.
On the other hand, since all such Hamiltonians are at most quadratic in momentum and position, the time evolution of the Wigner function must satisfy the Liouville equation [26,27], i.e., it must evolve classically in phase space.Therefore, If, on the other hand, [σ(X, P ) 1 , σ(X, P ) 2 ] = −[X, P ], then there is an extra matrix 1 0 0 −1 in the KAN decomposition of the linear part of σ.This operator corresponds to the antiunitary map ρ → ρ * , where * denotes complex conjugation in the position basis.Indeed, a short computation shows that Given an operator Ω, we define its Wigner function as With these choices of normalization, a short computation shows that tr(ρΩ) = R 2 dxdpW ρ (x, p)W Ω (x, p).
In the main text we are primarily concerned with operators of the form f (aX + bP + c) for some bounded measurable function f .We now prove that tr(ρf (aX Now, take ρ to be a coherent state, i.e., ρ = |α⟩⟨α|, with It follows that On the other hand, the Wigner function of a coherent state is known to be [28] with r 2 = x 2 + p 2 .Cancelling the factor e −|α| 2 in both sides of (A2) and expanding the remaining exponential in (A5) as a power series in α, ᾱ, we can compare the coefficients multiplying α m ᾱn on both sides of the resulting equation, thus obtaining Note that in [4] Tsirelson provides the complex conjugated formula for the same quantity.This mistake does not, however, invalidate the main result of [4], namely, the computation of the spectrum of a given linear operator.This is so because the spectra of a self-adjoint operator and its complex conjugate in a given basis coincide.
Next, we invoke (A6) to derive the matrix elements O nm (ϕ) := ⟨n| Θ(cos(ϕ)X + sin(ϕ)P ) |m⟩ and show that We will use this expression in Appendices C, F to lower bound the maximum quantum advantage in the standard and restricted projectile scenarios.

(A8)
We can evaluate the right-hand side of the above equation by changing to polar coordinates.The result is with where, in the last step, we changed the sum variable k → m + n − k so that a comparison with eq.(1.5) in [4] can be made.
As it turns out, the final expression for w nm can be written in terms of the generalized hypergeometric function p F q .Thanks to such an identity, we were able to compute w nm accurately for large values of m, n.
Appendix B: Properties of φ(α) First, we will calculate the kernel of Ω in position representation.Our starting point is the identity where the integral must be understood as a Cauchy principal value.We start from the easily verifiable identity: Applying the identity to eq. (B5), we find that On the other hand, when averaged over elements of S(α), A 0 has support on a two-dimensional subspace, namely, the span of the vectors The maximum eigenvalue of A 0 is therefore the result of solving the generalized eigenvalue problem min{λ : λG − F ≥ 0} with 2 × 2 matrices F, G given by The horizontal line is −0.0725.This spectrum was computed by numerically integrating with Mathematica, taking ε = 0.001 rather than a limit.
can be paired, i.e., dS = N µ(x)dx.In that case, after removing those, the remaining entries of vector ⃗ x are in the interval (x + dx, ∞).Also, the number of unpaired elements of ⃗ y greater than a − (x + dx) are Q + dQ, with dQ = N dx(ν(a − x) − µ(x)).If Q = 0, then there are two possibilities: (1) ν(a − x) ≥ µ(x), in which case N µ(x)dx entries can be paired, and so dS = N µ(x)dx, dQ = N dx(ν(a − x) − µ(x)); (2) ν(a − x) < µ(x), in which case just N ν(a − x)dx can be paired, and so dS = N ν(a − x)dx, dQ = 0. Defining q ≡ Q N , s ≡ S N , we hence have that the functions q(x), s(x) follow the system of differential equations Call (s(x), q(x)) the solution of the system of ordinary differential equations (F5) with the boundary condition s(−∞) = q(−∞) = 0. From all the above it follows that the solution of Problem 1 is p ⋆ c (µ, ν) = s(∞).
2. The gradient of p ⋆ c (µ, ν) Suppose that the distributions µ, ν depend on one parameter λ, i.e., µ = µ(x; λ), ν = ν(x; λ).We wonder how much p ⋆ c (λ) = p ⋆ c (µ(•; λ), ν(•; λ)) differs from p ⋆ c (λ + δλ), with δλ ≪ 1.Let us assume that the roots of q(x; λ) can be expressed as where An increment of λ will thus have two effects on p ⋆ c (λ).On one hand, the functions f + , f − will respectively change by the amounts ∂ ∂λ f + δλ, ∂ ∂λ f − δλ.On the other hand, the set of points x where q(x, λ+δλ) vanishes will change.Assuming that the kernel of q(•, λ + δλ) is of the form From eqs. (F5), and, assuming that µ, ν are smooth, we have that ν Indeed, if this quantity were negative, then q would have remained zero; and, if it were positive, then there would exist x − i < x < x + i such that ν(a−x; λ)−µ(x; λ) = 0, and q would have lifted itself from zero at x instead of x + i .By (F7), this implies that the last term of eq.(F8) vanishes.

Maximizing W(ρ)
Since the maximum quantum advantage is independent of the parameters M, ∆T, a, from now on we take M = ∆T = 1, a = 0 in W(ρ) that is W(ρ) ≡ ⟨Θ(X + P )⟩ ρ − p ⋆ c (ρ). (F15) We will only perform projected gradient ascend in the subspace H N = span{|n⟩ : n = 0, . . ., N }, noting that we can get better achievable lower bounds with increasing number of iterations and increasing N .For a learning rate ϵ, each iteration updates the state according to To solve this program, we used the MATLAB package YALMIP [32] in combination with the semidefinite programming solver MOSEK [33].
Given a set of states S and an operator A, we have, for any unitary U , that sup ρ∈S tr(ρA) = sup ρ∈U SU † tr ρU AU † .

FIG. 1 .
FIG. 1. Projectile scenario.A projectile is prepared at time t = 0 in [0, L] and, at time t = ∆T , we verify that it has reached region [a, ∞).Maximum quantum advantage in probability of arrival as compared to a classical particle is found to be the Bracken-Melloy constant, 0.0315 ≤ c bm ≈ 0.0384517 ≤ 0.0725.

FIG. 3 .
FIG. 3. Wigner functions of (left) near-optimal state for the projectile scenario and (right) conjectured-optimal state for the constrained projectile scenario.Both states are obtained by truncating to the harmonic oscillator energy level N = 170.The left state is the eigenstate of [Θ(−X)]170[(Θ(X + P ) − Θ(P ))]170[Θ(−X)]170 with eigenvalue 0.0331, where [C]N denotes the restriction of the operator C to the subspace spanned by the first N + 1 number states.The right state is the eigenstate of [Θ(X + P ) − Θ(X) − Θ(P )]170 with eigenvalue 0.1113.
f is bounded and measurable, we have that (as a tempered distribution) it has an inverse Fourier transform f , and we may write f (x) = R dt f (t)e ixt .Via functional calculus, we thus havef (aX + bP + c) = R dt f (t)e it(aX+bP +c) = R dt f (t)e itbP e itaX e −it 2 ab/2 e itc ,where we have used the Baker-Campbell-Haussdorf formula e i(ξX+ζP ) = e iζP e iξX e −i ξζ 2 .A straightforward computation now shows that⟨x − y/2| f (aX+bP +c) |x + y/2⟩ = R dt f (t)e it(ax+bp+c) , from which we conclude the result.Since it will be useful soon, we next compute the Wigner function of the operator |m⟩⟨n|, in number basis, i.e., (X + iP ) |n⟩ = √ 2n |n − 1⟩.First note that, by linearity of the Wigner function, for any state ρ = m,n ρ mn |m⟩⟨n|, we have that W ρ (x, p) = ρ mn W |m⟩⟨n| (x, p).(A2)

i 2 (y 2 −x 2 ) 4 X 2 , 2 − x 2 )
the final expression, we invoked the Baker-Campbell-Haussdorf formula e i(ξX+ζP ) = e iζP e iξX e −i ξζ 2 in the first line; the resolution of the identity I = dx |x⟩⟨x| = dx |y⟩⟨y|, in the second one; the relation ⟨x|p⟩ = e ipx / √ 2π (assuming that the bra is an element of the position basis; and the ket, of momentum basis), in the third one; and the relation dpe ips = 2πδ(s), in the fourth one.Using the same techniques, one finds that sign(be further reduced to a real kernel by conjugating it with the unitary e i |x⟩⟨y| , (B5) where sinc(z) := sin(z)/z.a. Bounding φ(α)

ρ
k+1 = P(ρ k + ϵ∇ ρ W(ρ)), (F16) with the projection P ensuring a valid density matrix.It remains to compute various quantities above.Firstly, for any matrix M , the projection P(M ) = argmin Z ∥Z − M ∥ 2 to the set of density matrix can be cast as a semidefinite program[31