Observation of the quantum Gouy phase

Controlling the evolution of a photonic quantum states is crucial for most quantum information processing and metrology tasks. Because of its importance, many mechanisms of quantum state evolution have been tested in detail and are well understood. However, the fundamental phase anomaly of evolving waves called the Gouy phase has not been studied in the context of elementary quantum states of light such as photon number states. Here we outline a simple method for calculating the quantum state evolution upon propagation and demonstrate experimentally how this quantum Gouy phase affects two-photon quantum states. Our results show that the increased phase sensitivity of multi-photon states also extends to this fundamental phase anomaly and has to be taken into account to fully understand the state evolution. We further demonstrate how the Gouy phase can be used as a tool for manipulating quantum states of any bosonic system in future quantum technologies, outline a possible application in quantum-enhanced sensing, and dispel a common misconception related to the nature of the increased phase sensitivity of multi-photon quantum states.

Controlling the evolution of a photonic quantum states is crucial for most quantum information processing and metrology tasks. Because of its importance, many mechanisms of quantum state evolution have been tested in detail and are well understood. However, the fundamental phase anomaly of evolving waves called the Gouy phase has not been studied in the context of elementary quantum states of light such as photon number states. Here we outline a simple method for calculating the quantum state evolution upon propagation and demonstrate experimentally how this quantum Gouy phase affects two-photon quantum states. Our results show that the increased phase sensitivity of multi-photon states also extends to this fundamental phase anomaly and has to be taken into account to fully understand the state evolution. We further demonstrate how the Gouy phase can be used as a tool for manipulating quantum states of any bosonic system in future quantum technologies, outline a possible application in quantum-enhanced sensing, and dispel a common misconception related to the nature of the increased phase sensitivity of multi-photon quantum states.
The wave dynamics dictating the evolution of quantum states is of utmost importance in both fundamental studies of quantum systems and quantum technological applications. For photons, the evolution of their spatial structure has been the key in a plethora of promising techniques for quantum communication [1,2], information processing [3,4], simulation [5], and metrology [6]. One particular feature of a converging wave travelling through its focus is the acquisition of an additional phase shift when compared to a collimated beam or a plane wave traveling the same distance. This effect, which is known as the Gouy phase, was first observed and described by Gouy more than a century ago [7,8]. Although the phenomenon is well established and can be described through methods in physical optics [9,10], the Gouy phase continues to be the topic of studies discussing its underlying physical origin by linking it to properties such as the geometry of the focus, geometric phases, and the uncertainty principle [9,[11][12][13][14][15][16][17][18]. In addition to the continued interest aiming at providing an intuition for the phenomenon, this phase anomaly is often harnessed to realize novel tools in optics [19][20][21][22].
Despite the Gouy phase being a general wave phenomenon, studies investigating its role in quantum state evolution have been limited to a few matter wave studies [23][24][25][26][27] and spatially separated photon pairs [28,29]. While these demonstrations utilize (locally) single quantum systems and, thus, observe the effect known for classical light waves, more complex quantum states consisting of multiple identical quantum systems, i.e., bosonic systems with multiple excitations, have not been studied before. We term the specific phase acquired by such quantum states the quantum Gouy phase.
In general, any phase accrued by a mode of a photonic quantum system leads to a photon-number dependent phase for the quantum state. This means that whereas a single photon or a classical field would acquire a phase φ upon propagation, when N -photons occupy the same mode (|N ) the quantum state is left with N times the same phase, i.e. exp(iN φ) |N [30]. This increased phase sensitivity of photon number states is utilized in so-called N00N states that have garnered popularity due to their potential in pushing the sensitivity of measurements to what is considered the absolute physical limit [31]. N00N states can be compactly expressed for two orthogonal modes p and p as Hence, the enhancement in measurement sensitivity is enabled by the phase difference between the two components being N times the phase difference between the underlying modes. More importantly, using such a N00N state configuration allows the study of the speed-up of the quantum Gouy phase compared to the classical case.
In the present work, we describe theoretically how an N -photon number state evolves upon propagation and verify experimentally the speed-up of the quantum Gouy phase with two-photon N00N states through interference in the transverse structure of a bi-photon. We further show that the quantum Gouy phase speed-up can be applied to super-resolving longitudinal displacement measurements using the quantum Fisher information (QFI) formalism and solidifying its link to the uncertainty interpretation of the Gouy phase [12]. Finally we show that our results for N-photon states cannot be simulated by classical light with a λ/N wavelength, demonstrating that the often-used effective de Broglie wavelength approach for multi-photon states, although useful in specific cases [32][33][34], is not always accurate. As such, our work brings the fundamental wave feature of the Gouy phase to the quantum domain, thereby opening the path to its utilization in quantum technological applications FIG. 1. Observing the quantum Gouy phase through a changing on-axis interference along the propagation direction. a) Conceptual image of the observation scheme. The image displays the intensity structure of a superposition of a radial mode with p = 4 and a Gaussian reference (p = 0) at different distances from the focus. The inset shows the intensity of the field on the optical axis. b) Intensity of a classical light beam prepared in the same superposition as in a), with a Gaussian waist w0 = 25 µm and z0 = 0. c) Spatially varying two-photon probability for a two-photon N00N state prepared in the same radial modes as in a) and b). To make this structure visible, we post-select for cases where the two photons exist in the same position using the projection P (x, y, z) = | Ψ(z)|2 x,y | 2 . For b) and c) the intensities/probabilities are calculated on a plane cutting through the optical axis (see a) for reference).
through its unique quantum state manipulation properties.

PROBING THE QUANTUM GOUY PHASE
To observe the quantum Gouy phase of N -photons, an interferometric measurement scheme can be used. We chose to use the transverse-spatial modes of paraxial light beams as the different "arms" of the interferometric scheme, where one mode acts as the required reference "arm". More specifically, we used Laguerre-Gaussian (LG) modes, which are a family of orthogonal solutions to the paraxial wave equation in cylindrical coordinates [35]. In the case of a classical monochromatic field, the Gouy phase of these modes evolves as [35] where z is the propagation distance, k is the wavenumber, is an integer giving the number of orbital angular momentum quanta per photon, p is a positive integer defining the radial transverse structure of the field, w 0 is the beam waist defining the transverse extent of the beam at its focus, and z 0 gives the position of the beam focus along the optical axis. Since the Gouy phase depends on the mode order S = 2p + | | + 1, its anomalous phase behaviour can be observed through the change of the transverse structure during propagation when the light is in a superposition of spatial modes of different mode orders [36]. For radial modes, which are LG modes with = 0, this change results in a varying intensity along the optical axis, as can be seen in Fig. 1a). Thus, to probe the quantum Gouy phase and contrast it to its classical counterpart, we study the superposition of a Gaussian reference mode (p = 0) and different higher order radial modes in both the classical domain and the aforementioned quantum setting, i.e., a N00N state superposition. By measuring the change in intensity and two-photon detection rate, respectively, observed in a single mode fiber (SMF) scanned through the focus, we are able to directly observe the speed-up of the quantum Gouy phase.

THEORETICAL EVOLUTION UPON PROPAGATION
In our measurement scheme, we expect the propagation to result in a photon number dependent Gouy phase when the state |N p is translated through a focus. To verify these expectations theoretically, we start with N photons occupying a monochromatic paraxial mode at a position z = 0, with a complex field structure u p (ρ, 0). To translate the mode along the optical axis, we apply the translation operator e iPzz/ to the mode in the angular spectrum representation, in which the quantized mode of light can be expressed aŝ where F p (κ, 0) represents the normalized complex amplitude of the plane wave mode with transverse wave vector κ, andâ † (κ) is the corresponding operator density [37,38]. After applying the translation operator, the mode takes the form which is identical to the initial mode being propagated by z using the angular spectrum method (ASM) [9,39].
We thus see that the quantized mode evolves identically to a classical light field, i.e., the propagated LG mode has an identical spatial structure u p (ρ, z) only differing by the propagation-related change in wavefront curvature and beam radius. Due to the beam evolving according to the ASM, we can extract the Gouy phase evolution by defining a new modeb † p (z) which has the structure of the field after translation, without the accumulated Gouy phase, i.e., u p (ρ, z)e −ikz e iΦ G (z) . Using this new mode, we can express the mode after propagation as a single mode with a phasê We can then simply state the Gouy phase evolution of an N -photon Fock state as which explicitly contains the photon number dependent Gouy phase evolution. For a detailed derivation, see the Supplementary.

EXPERIMENT
In the experiment, we first prepared laser light in a superposition of the Gaussian reference mode and one of the higher-order radial modes. The structuring of the laser beam was done with a single hologram on a spatial light modulator (SLM), using a holographic method commonly known as mode carving [40]. After structuring, the beam was imaged one focal distance away from a 75 mm lens which performs an optical Fourier transform on the transverse structure while focusing [39]. Since the transverse structure and its Fourier transform are identical for LG modes, the beam structure at the focus was identical to the structure carved at the SLM, up to a phase factor of π between the superposed LG modes which needed to be accounted for with odd values of the radial index [19,20]. To measure the Gouy phase induced change in the interference along the optical axis, we placed an SMF at the focus and moved it longitudinally using a stage with a computer controlled piezo actuator. The laser source was a continuous-wave diode laser operating at 810 nm and the SLM used for Two photons with Gaussian beam profiles were sent on to separate sections of an SLM where they were independently structured into orthogonal superpositions of radial modes. Then these photons were probabilistically overlapped using a beamsplitter, after which they bunched into a radial mode N00N state [6]. Finally, this two-photon N00N state was focused down to a 25 µm Gaussian beam waist and coupled into a SMF (with a mode field diameter of 5 µm) that was scanned through the focus (from behind the focus towards the lens). The twophotons were then probabilistically split into two single photon avalanche diodes (SPAD) and we post-selected on both of the detectors detecting a photon at the same time using a coincidence counter (CC). For more details, see the main text and Supplementary. structuring the light was wavefront corrected using the method described in Ref. [41]. Furthermore, to get the generated modes as close as possible to the correct transverse structure at the wanted beam radius, we employed an additional Gaussian correction in the mode carving that minimized any effect of the initial Gaussian beam structure in the carved mode (see Supplementary).
For a classical field, we can extract the theoretically expected measurement results simply by calculating the overlap of the Gaussian eigenmode of the SMF and the normalized transverse structure of the scalar field u tot (ρ, z) = 1 √ 2 (u 0p (ρ, z) − e iθ u 0p (ρ, z)). Thus, for laser light, the amount of laser power coupled in to the fiber is proportional to where A j (z) refer to the overlap between the normalized radial mode j, at a distance z from its focus, and the normalized Gaussian eigenmode of the fiber. To see the Gouy phase dependence of the detection probability, the above equation can then be stated as where the term φ(z) is an extra phase contribution from the curvature of the wavefront acquired upon propagation. However, since the wavefront curvature is very small near the optical axis, the only significant contribution to the phase of the overlaps A j (z) comes from the Gouy phase difference ∆Φ G (z). Thus, scanning the fiber through the focus results in a signal which oscillates as cos 2(p − p) arctan 2(z−z0) kw 2 0 underneath some envelope function caused by the z-dependence of the overlap functions.
For the measurements, we kept the reference mode, i.e., a Gaussian mode with radial index p = 0, fixed and varied the index p of the probe mode between 1 and 4, which lead to four different measurement scenarios with differing Gouy-phase contributions. The measured data can be found on the top row of Fig. 3. The measurements follow the probability introduced above very well, which we verified by fitting curves matching Eq. (7) to the data. In each fit, we fixed the mode field diameter of our fiber to the 5 µm specified by the manufacturer and only had 4 fitting parameters: an overall scaling factor of the function, the beam waist w 0 , focal position z 0 , and the z-independent phase offset θ. The average adjusted R-squared value of the fits was 0.986, meaning that the data corresponds well to the theoretical model.
After first verifying the methods viability using a laser and showing the effect of the Gouy phase on a classical interference pattern along the optical axis, we extended the measurement scheme to observe the quantum Gouy phase. Following the same general idea, we now generated different two-photon N00N states between a ref-erence Gaussian mode (p = 0) and higher order radial modes, and studied the two-photon interference pattern along the optical axis. To prepare such a N00N state, we first generated photon pairs through spontaneous parametric down-conversion (see Supplementary for more information) and then shaped each of the two photons individually into a well-defined superposition of the wanted radial modes using two holograms performing two different mode carvings. Once each of the photons was structured, we directed the photons into the same beam path using a beamsplitter. As demonstrated in [6], once in the same beam path, indistinguishable photons bunch into the wanted spatial mode N00N state given in Eq. (1). A simplified sketch of the two-photon experimental setup can be seen in Fig. 2.
To calculate the N -photon coincidence probability, we project the radial mode N00N state |Ψ(z) onto the state where all of the photons have been coupled successfully into the SMF P = | Ψ(z)|N SM F | 2 . Assuming that we produce perfectly balanced N00N states of radial modes with a phase offset θ, the N -photon detection probability can be reduced to the form As before, we can express this coincidence probability as which is similar to the detection probability of the classical field, leading to an oscillating interference underneath some envelope function. However, in the above equation we see the photon number dependent scaling for both the frequency of the oscillation as well as the envelope term. Note that a probability curve with half the amplitude but same shape can also be observed for photon pairs prepared similarly without bunching. Thus, to verify that we generate radial mode N00N states in our experiment, we prepared the two photons in the corresponding radial mode superpositions and showed that the probability of coupling both of the photons into the SMF roughly doubles when the photons are made indistinguishable in time, which is a clear signature of bunching (see Supplementary for the measured data). For detailed derivations of the detection probabilities, see the Supplementary. For the N00N state measurements, we used the same set of radial modes in superposition with the reference Gaussian mode leading to the data shown on the bottom row of Fig. 3. As before, the data follows very well the theoretically expected curves, verifying the abovepresented equations and their described behaviours. Fits of Eq. (8) to the data, with the same parameters as in the classical case, resulted in an average adjusted R-squared value of 0.951. The slight imperfections in the data can all be accounted for by imperfections in the alignment, imaging, the SMF eigenmode, spatial mode generation, and errors in the stage position. Besides the errors in the stage positions, all of these can be effectively categorized as contaminations of our state space by modes not included in the theoretical analysis. Hence, our results demonstrate that the quantum Gouy phase leads to a speed up in the accumulated phase upon propagation and also modulates the underlying envelope function. As we will discuss next, both features shed new light on the fundamental understanding of the Gouy phase, as well as hint at quantum enhanced metrology applications.

QUANTUM FISHER INFORMATION
As the quantum Gouy phase evolves faster with a larger number of photons, one application could be super-sensitive measurements of longitudinal displacement. This prospect can be investigated by calculating the QFI achieved through translation, which is of the form [31,[42][43][44] When calculating this variance for the radial mode N00N state |Ψ(z) we get the QFI where ∆ 2 k z | i and k z i are the variance and the average of k z for the mode i, respectively, calculated using the angular spectrum of the corresponding mode. It is worth noting that the QFI does not depend on z, since the angular spectrum of a mode only acquires a phase structure upon translation. From Eq. (10) we can see that the second term of the QFI has Heisenberg scaling. As we show in the Supplementary, this term relates to the Gouy phase difference between modes p and p . Hence, radial mode N00N states along with their quantum Gouy phase properties should be able to enhance the sensitivity of longitudinal displacement measurements. However, although these states provide benefits such as intrinsic interferometric stability when translating the mode along z, the spatial extent of the modes changes, making it challenging to device a real measurement capable of saturating the QFI at any z. The form of Eq. (10) also shows that it could be possible to engineer different spatially structured quantum states to measure different physical parameters. Due to the form of the quantum Fisher information, the key feature that needs to be optimized in such state engineering should be maximizing variance of a specific momentum of the quantum state. For example, this would mean maximizing the variance in orbital angular momentum for rotation sensing [6] or linear momentum for sensing the longitudinal position (9). See the Supplementary for the derivations of the QFI and the Fisher information calculated for the projection used in our experiment.

MOMENTUM UNCERTAINTY
In addition to showing the potential for Heisenberg scaling, there is an interesting connection between the QFI and the uncertainty interpretation of the Gouy phase which fundamentally links the potential change in the spread of the transverse momentum to the evolution of the Gouy phase [12]. Feng and Winful also noted that a larger momentum spread of higher-order modes results in a bigger Gouy phase shift [12]. Since the Gouy phase is increased by the photon number N , which is accompanied by a photon number dependent momentum spread, as can be seen in Eq. (10), our results make a further connection between the quantum Gouy phase and its uncertainty interpretation. Similarly to Ref. [12], one can further link this behaviour to a tighter spatial confinement of the photons which can be made visible, e.g. by measuring the spatial extent of the N -photon state as shown in Fig. 1c).

DE BROGLIE WAVELENGTH OF LIGHT
Finally, our results show that the behaviour of a two-photon N00N state cannot be replicated simply by switching to a classical field with half the wavelength. The difference is clear if we note that the Gouy phase has a nonlinear dependence on the wavenumber which means that simply ascribing an effective de Broglie wavelength λ/N to the N -photon state does not produce the correct quantum Gouy phase. This is in contrast to the phase accrued by a non-converging field upon propagation and arguments discussed in such a context [32][33][34]. In order to investigate this fundamental difference in more detail, we have plotted in Fig. 4 the measured data for two radial mode N00N states, along with overlap curves calculated for classical 405 nm modes with two different mode orders and waists. From these comparisons we see that the effect is not reproduced by a simple switching of the wavelength or doubling of the mode order.
Based on the comparison in Fig. 4 and Eq. (6) the only exact description of the N -photon Fock state evolution seems to be that it evolves as the underlying mode, taken to the power of N . Although, doubling the mode order and halving the wavelength seems to replicate quite well the shown two-photon behaviour. Since the state evolves as the mode taken to poower N , this evolution of the N -photon quantum state results in a more rapid phase change and tighter confinement of the N -photon. Both of these features have been taken advantage of in different studies and experiments. Either in the form of N00N-state super-resolution measurements [30,45] or in increasing the confinement [46].

CONCLUSION
In summary, here we have verified in theory and experiment that the increased phase sensitivity of multiphoton quantum states also extends to the fundamental phase anomaly of converging waves called the Gouy phase. We have shown through single-path interferometric measurements along the optical axis, that two-photon N00N states experience twice the Gouy phase when traveling through a focus. Since the Gouy phase is a fundamental feature of converging waves, our results should apply broadly to quantum states of any bosonic system. Moreover, as the Gouy phase is an important factor in systems such as optical cavities [46,47], and a powerful tool in various applications such as mode sorters and mode converters [19][20][21], our results can be widely utilized in applications in quantum optics and quantum in- . The curves match well near the focus. However, the blue curve does not exhibit the same fringe pattern and the red curve has a larger relative amplitude outside the focal region. Hence, the quantum Gouy phase behaviour cannot be exactly reproduced by simply changing the wavelength and mode order. formation science. In addition to providing a tool for quantum state manipulation, we showed that our results allow Heisenberg-limited scaling in measurements of the longitudinal displacement and, as such, might inspire new superresolution measurement schemes.
Besides these possible technological applications, we have linked the speed-up of the Gouy phase in the quantum domain to an increased spread in the momentum of an N -photon state. Hence, our results show that the uncertainty interpretation of the phase anomaly [12] holds true in the quantum domain. Finally, due to the nonlinear relation between the Gouy phase and the wavenumber, our results unambiguously demonstrate that an Nphoton state cannot be rigorously modelled by using a classical field with a wavelength λ/N . However, our results suggest that an additional N -fold increase in the mode order can reproduce the effect of the quantum Gouy phase when the beam Rayleigh lengths are matched. This hints at a possible link between an N -photon state and the N th harmonic of a classical field, which introduces an increase of the mode order and decrease of the beam waist, in addition to doubling the frequency. Thus, our study not only outlines possible applications using the quantum features of spatially structured photons, it also sheds new light on the fundamental understanding of the Gouy phase, a property intrinsic to all systems described by converging or diverging waves. from the Doctoral School of Tampere upon a translation of z along the optical axiŝ whereâ † p (z) represents the transverse mode at any plane z characterized by the integers and p, ρ contains the transverse coordinates,â † (ρ) is an operator density, and u p (ρ, 0) is the normalized structure of the mode. The operatorP z is the longitudinal component of the linear momentum operator, defined for a monochromatic field in the plane wave basis asP where k z (κ j ) is the longitudinal wave vector for the mode j, which depends on the transverse wave vector κ j as k z (κ j ) = k 2 − ||κ j || 2 , with k being the wave number, common to all of the plane wave modesâ † (κ j ). To operate on our mode with this unitary, we initially express our mode in the plane wave basis aŝ where F p (κ, 0) is the normalized angular spectrum of the mode p at z = 0. Now we can restate Eq. (S2) then using the Baker-Hausdorff lemma, we can finally obtain The form of Eq. (S6) is now identical to a mode with an angular spectrum F p (κ, 0) propagated by a distance z using the angular spectrum method (ASM) [3,4]. Since the ASM accounts for every aspect of paraxial beam propagation, including the Gouy phase and the plane wave phase e −ikz , we can extract the Gouy phase out of the translated mode by initially expressing the translated mode in real space using the complex field distribution where u p (ρ, z) is the normalized transverse field for the mode, calculated classically through the ASM. We can then choose an artificial set of orthogonal spatial modesb † p (z) which have exactly the structure of the p modes at the position z, without the Gouy phase, i.e., u p (ρ, z) = u p (ρ, z)e −ikz e iΦ G . Expressing the position basis mode density in this new mode basisâ † (ρ) = p u * p (ρ, z)b † p (z) the translated mode takes the form Due to the orthonormality of the chosen spatial mode basis [5], the above integral reduces tô Hence, we can then state the evolution of an N-photon Fock state in the spatial mode p as showing that upon translation, any phase accumulated by the mode results in N-times the same phase being accumulated by the N-photon Fock state. Interestingly, from the fact that any arbitrary modeb † z can be chosen, one can infer that any changes in the amplitude of the field during propagation is also magnified N-times. Similarly to the N-times increase in the phase however, in order to observe this N-fold change in the amplitude the measurement needs to be chosen in a manner that displays this change. One example of such a measurement is the post-selection on both of the photons existing in the same spatial position which manifests as an increased confinement in Fig. 1c) of the main text. Similarly, the same can be seen in the measured data in Fig. 3 as a narrowing of the envelope within which the fringes are observed. As a result, these changes could be summarized as the state experiencing any change in the mode to the power of N, which is in conflict with the effective de Broglie wavelength interpretation whenever the amplitude changes or a phase is acquired that has a nonlinear dependence on the wavenumber.

DERIVATION OF THE MEASUREMENT PROBABILITY
To model the expected measurement signal we calculate the probability of N photons coupling into a single mode fiber (SMF) from a radial mode N00N state. We start off with the radial mode N00N state at position z where θ is a constant phase offset between the two terms. By decomposing the eigenmode of the SMF into a superposition of LG modes at a distance z from the beam waist we can calculate the probability of N photons coupling into the SMF using (S13) Note that in the above equations, using the normalized transverse structures, the overlaps are defined as which also includes the Gouy phase of each mode (note that for A p (z), = 0). From the overlap calculation in Eq. (S13), only the states with N photons in either mode 0p or 0p survive and the above equation simplifies to (S15) To get a more intuitive expression of this equation, we can further state is small since the SMF only probes the phase difference close to the optical axis where the phase of the overlap is mostly dependent on the Gouy phase and not the changing wavefront curvature. By then defining the Gouy phase difference as ∆Φ G (z) = Φ G − Φ G (Φ G and Φ G correspond to Gouy phases acquired by LG modes of different mode order S = 2p + 1 and S = 2p + 1) we can write the probability as (S16) From the form of this equation, it is clear that the probability amplitude oscillates according to cos(2N (p − p) × atan(z/z R )), meaning that the oscillations get less frequent the further away we go from z = 0. In addition to this, the oscillation happens inside an envelope function defined by the overlaps |A p | and |A p |.
To derive a similar expected signal for the case of a classical monochromatic field with the normalized transverse structure u tot (ρ, z) = 1 √ 2 (u 0p (ρ, z) − e iθ u 0p (ρ, z)) being coupled into the SMF, we can calculate the power of the light in the SMF as an overlap between the structure of the incident field and the normalized eigenmode of the SMF where the overlaps A i (z) correspond to the same overlaps calculated in Eq. (S14) and the same assumptions can be made about the overlaps near the optical axis.

OVERLAP BETWEEN THE EIGENMODE OF A FIBER AND A LAGUERRE-GAUSSIAN MODE
To speed up the data processing, we analytically derived the overlap of a monochromatic p-mode and a Gaussian mode corresponding to the eigenmode of the SMF where the normalized p-mode is defined according to [6] u 0p (r, ϕ, z) = 2 π where r is the radial coordinate, ϕ is the azimuthal coordinate, w(z) = w 0 1 + [(z − z 0 )/z R ] 2 is the beam radius, 2 is the Rayleigh length, and k is the wave number. The normalized eigenmode of the fiber can be similarly defined as For simplicity we have set z 0 = 0 and insert the above definitions into equation (S18) If we then define B(z) = w(z) , and x(z) = 2r 2 w 2 (z) , we can simplify the above equation to We can then use the identity ∞ 0 y n exp [−ay] dy = n! a n+1 , which holds when Re(a) > 0 and n ∈ N, to arrive at the expression which we then used for calculating the overlaps in equations (S16) and (S17).

DERIVATION OF THE QUANTUM FISHER INFORMATION
The quantum Fisher information (QFI) carried by a probe state |ψ(x) about a parameter x, encoded in the state by the unitary evolutionÛ (x) is given by [7][8][9] whereĤ is the generator of the unitaryÛ and its variance ∆ 2Ĥ = Ĥ 2 ψ − Ĥ 2 ψ is taken with respect to the input state. In our case, where z is the longitudinal displacement, the translation operatorÛ = e iPzz/ is the unitary of interest, yielding Thus, to obtain the QFI we need an expression for the variance of the longitudinal momentum operator, which we provide here for the radial N00N state (S11).
We start by using equations (S3) and (S11) to obtain where we used and is the classical average of the longitudinal wave vector k z for the mode p. Similarly, we obtain that The resulting QFI is, therefore where ∆ 2 k z p = k 2 z − k z 2 , and equivalently for p .
To obtain further insight on the QFI, we swap the modes p and p to the Hermite-Gaussian (HG) modes HG mn and HG m n , respectively. In the HG basis, the average and variance of the longitudinal wave vector have the simple expressions resulting in where S = m + n is the mode order, S 2 = m 2 + n 2 (with equivalent expressions for S and S 2 ), and ∆Φ G (z) = (S − S ) tan −1 (z/z r ) is the Gouy phase difference between the modes. The QFI is z-independent, consistent with the self-similarity of the angular spectrum upon propagation. It is worth noting that the second term of the QFI remains unaltered in the LG basis, while a closed form for the first term can be obtained by decomposing the LG modes in the HG basis [10].
We note that the QFI is comprised of two fundamentally distinct terms. The first, standard-quantum limited (proportional to N ), increases monotonically with the indices of the modes p and p , and thus with dimension of the state space. We hypothesize that this term carries information about the full distribution of photons in the transverse plane, such that it's retrieval cannot be fully achieved by interferometric measurements alone. On the other hand, the second term displays Heisenberg-limited scaling (proportional to N 2 ) and increases with the slope of the Gouy phase difference at the focus Equation (S31), therefore, suggests that measurements of the quantum Gouy phase with N00N states can reach quantum-optimal sensitivity in the estimation of longitudinal displacements.
To support our claim, we also calculate the classical Fisher information using the modeled measurement probability (S16) without the wavefront curvature. Assuming |A p | = |A p | = A independent of z, around the region of interest z = 0, we obtain [11] where ] is the QFI term proportional to N 2 , P 1 is given by (S16) and P 2 = 1 − P 1 . The maximum value of the Fisher information is obtained at the focus, and it's given by The above equation shows that the proposed measurement setup is indeed sensitive to the Heisenberg-limited part of the QFI, but does not capture any of the information contained in the SQL part. It is worth noting that equations (S33) and (S34) are valid in the limit where the SMF mode is much smaller than the input modes, such that the contribution from the wavefront curvature to the measurement probability is negligible. In this case, the high losses yield A 1 and the Fisher information is far from reaching its quantum bound.

EXPERIMENTAL DETAILS
The experimental setup consists of a photon pair source and the spatial mode manipulation part (see Fig. S1 for a detailed schematic). In the photon pair source we use a 405 nm continuous-wave pump laser focused down to a 12 mm long, periodically poled nonlinear crystal made out of potassium titanyl phosphate (ppKTP). The pump laser had a slightly astigmatic focus with an approximately 67 µm Gaussian beam waist in the crystal. In the crystal, some of the 405 nm photons go through spontaneous parametric downconversion (SPDC) which is type 0 phase matched. The frequencies of the two emerging 810 nm photons were made degenerate by tuning the phase matching by controlling the temperature of the crystal. After the SPDC, the pump laser is filtered out using two bandpass filters (BPF), where the first one has a 10 nm bandpass and the second one a 3 nm bandpass around 810 nm. The first filter is used to remove the pump laser and the second filter is used to tune the spectral properties of the photon pair. To split each pair of photons, we use their momentum anti-correlations by placing a lens one focal length away from the crystal after which we place a D-shaped mirror one focal length behind the lens. The lens then performs an optical Fourier transform and the momentum anti-correlations of the photons are used to split each pair. After splitting the photons, one of them is sent through a delay stage which controls the temporal indistinguishability of the two photons. The photons are then coupled into separate SMF's placed on coupling stages (xyz) using a lens and a microscope objective.
FIG. S1. A detailed drawing of the experimental system. The photon pair source (left) generates photon pairs which are split and coupled into SMF's. The source also includes a delay stage which can tune the temporal distinguishability of the photon pair. The photons are sent to the spatial mode manipulation setup using the SMF's after which they are sent onto two separate regions of a spatial light modulator (SLM). After the photons have been strucutred at the SLM they are imaged with 4f systems to roughly one focal length away from the 75 mm lens which is used to focus the photons onto a SMF placed on a stage. Within these identical 4f systems we place a beamsplitter (BS) which is used to probabilistically combine the photons into the same beam path. The half wave plates (HWP) and quarter wave plates (QWP) are used to facilitate optimum efficiency of the polarization dependent SPDC process and the SLM.
The photons then exit the SMF's after which they are collimated and sent through waveplates which align their polarizations such that the polarization sensitive SLM operates at optimum efficiency. The SLM used was a Holoeye Pluto 2 and it was wavefront corrected using the method described in [12]. The two photons are then modulated on two separate regions of the SLM, using mode carving which is a technique where amplitude and phase modulation can be performed on a single phase-only hologram [13]. In the process, a transverse field structure |u A (x, y)|e −iΦ(x,y) is carved onto the first diffraction order of the hologram, out of a larger Gaussian input beam. Quite often the effects of the initial Gaussian structure of the input beam are removed by making the Gaussian large enough to have an effectively flat amplitude distribution in the area of the hologram. However, to produce the best possible modes at the beam radius required for our measurements, we removed the structure of the initial Gaussian by generating a hologram for the field |u A (x,y)| u0(x,y;win) e −iΦ(x,y) using the same mode carving technique. This effectively removes the structure of the incident Gaussian u 0 (x, y; w in ) in the first diffraction order, and we only need to measure the radius w in of the beam and do not have to worry about making the incident Gaussian much larger than the hologram. However, we note that the incident Gaussian still needs to be slightly larger than the hologram for optimal performance. For more information on this added Gaussian correction see the supplementary of [14].
At the SLM, we structure one of the photons in the p-mode superposition 1 √ 2 [u 00 (ρ, 0) − e iθ1 u 0p (ρ, 0)] and the other one in the superposition 1 √ 2 [u 00 (ρ, 0) + e iθ2 u 0p (ρ, 0)]. Here the phases θ 1 and θ 2 correspond to small corrections that need to be made due to slight imperfections in the imaging. Due to these same imperfections, we employed additional lens terms on one of the holograms to match the focal points of the two beams as closely as possible. After structuring, both photons go into identical 4f imaging systems and they are also sent into the same beam path, probabilistically, using a 50:50 beamsplitter. Once in the same beam path the photons can bunch into the radial mode N00N state, as long as they are indistinguishable in every degree of freedom, which we verified in our experiments as shown below. After this, the photons are focused down using a 75 mm focusing lens on to a SMF (Thorlabs 780HP FC/PC) placed on a translation stage that scans the SMF through the focus using a computer controlled piezo actuator (Thorlabs PIA13). The final SMF is placed on a coupling stage (xyz) connected to a mount which can control the tip and tilt of the SMF. After the photons are coupled into the SMF, we split them using a 50:50 fiber beamsplitter and sent them into separate single-photon avalanche diodes (SPAD; laser components COUNT T) from which we post-selected on coincident detections of two photons using a coincidence counter (CC; IDQ ID900). The accidental coincidences were removed from all of the data. We calculated the accidental rates as R 1 R 2 τ , when R i correspond to the measured single photon rates in each detector and τ = 1 ns is the coincidence window used.
In the measurements we initially set the stage as far away from the focus as possible and moved it closer to the lens, scanning the focus in discrete steps. When calculating the distance between subsequent measurement points, we used the typical step size provided for the piezo actuator (20 nm per piezo step). However, the manufacturer states that this step size can vary up to 20 % depending on the component variance, change of direction, and application conditions. Thus we tried to keep the conditions of the lab consistent while only scanning the piezo in one direction for all of the measurement.
For the measurements with a laser, an 810 nm continuous-wave laser was used and we only required one input hologram. We also replaced the coincident detection scheme with a power meter before and after the setup. The results were then calculated as the power measured after the last SMF normalized by the power measured before the laser was coupled out of the input SMF. The reference power was measured continuously by splitting the input light field using a fiber beamsplitter. Additionally, when structuring the laser field into the superposition 1 √ 2 [u 00 (ρ, 0) − u 0p (ρ, 0)] at the measurement fiber, we had to take into account the differing Gouy phases of odd and even order modes. Hence, when p was an odd integer we had to generate the field 1 √ 2 [u 00 (ρ, 0) + e iθ1 u 0p (ρ, 0)] instead of 1 √ 2 [u 00 (ρ, 0) − e iθ1 u 0p (ρ, 0)] to compensate for the π phase difference in Gouy phase when performing the optical Fourier transform of the field.
FIG. S2. Measured photon bunching curve with two photons in orthogonal radial mode superpositions. Accidental coincidences have been removed from the data and the errorbars (±1 standard deviations) are calculated from 21 repetitions of the measurement. The change in stage position should be multiplied by two to arrive at the effective path length change (see Fig. S1 for the configuration of the mirrors on the delay stage). The results show that the rate at which the photon pairs couple into the SMF roughly doubles when the photons bunch into a radial mode N00N state.
To verify photon bunching in to radial mode N00N states, we set the final SMF close to the focus of the field and generated a radial mode N00N state 1 √ 2 [|2, 0 0,2 + |0, 2 0,2 ] by individually structuring the photons into the superpositions 1 √ 2 [u 00 (ρ, 0) − iu 02 (ρ, 0)] and 1 √ 2 [u 00 (ρ, 0) + iu 02 (ρ, 0)], as described earlier. We then scan the delay stage in the source to vary the distinguishability of the photons. From the results of the scan (Fig. S2) we see that, as expected, the rate of photon pairs coupling into the SMF roughly doubles when the photons are made indistinguishable.