Higher dimensional quantum communication in a curved spacetime: an efficient simulation of the propagation of the wavefront of a photon

TL;DR

This paper presents an efficient numerical simulation method for the propagation of complex photon wavefronts in inhomogeneous media, with applications to higher-dimensional quantum communication and potential satellite QKD systems.

Contribution

It introduces a novel simulation approach for complex wavefronts in inhomogeneous media, enabling analysis of high-dimensional quantum communication channels in curved spacetime environments.

Findings

Simulated the propagation of optical vortices in holographic elements.

Analyzed an OAM volume-holographic sorter.

Proposed optical analog models for wavefront evolution in Earth's curved spacetime.

Abstract

A photon with a modulated wavefront can produce a quantum communication channel in a larger Hilbert space. For example, higher dimensional quantum key distribution (HD-QKD) can encode information in the transverse linear momentum (LM) or orbital angular momentum (OAM) modes of a photon. This is markedly different than using the intrinsic polarization of a photon. HD-QKD has advantages for free space QKD since it can increase the communication channel\~Os tolerance to bit error rate (BER) while maintaining or increasing the channels bandwidth. We describe an efficient numerical simulation of the propagation photon with an arbitrary complex wavefront in a material with an isotropic but inhomogeneous index of refraction. We simulate the waveform propagation of an optical vortex in a volume holographic element in the paraxial approximation using an operator splitting method. We use this code to analyze an OAM volume-holographic sorter. Furthermore, there are analogue models of the evolution of a wavefront in the curved spacetime environs of the Earth that can be constructed using an optical medium with a given index of refraction. This can lead to a work-bench realization of a satellite HD-QKD system.

Equations

$$-\frac{n^2}{c^2}{E}_{,tt}=-{\nabla }^{2}E-2\nabla \left(E\cdot \nabla \left(\log (n)\right)\right).$$

$$n=\frac{\sqrt{{\mu }_0ϵ}}{\sqrt{{\mu }_0{ϵ}_0}}=\frac{c}{v},$$

$$E(r,t)=u(r)T(t),$$

$${\nabla }^{2}u+\underset{\approx 0}{\underbrace{2\nabla \left(u\cdot \nabla \left(\log (n))\right)\right)}}+n^2{k}_0^{2}u=0,$$

$$T(t)\propto e^{i{\omega }_{0}t}.$$

$${\nabla }^{2}u+k^{2}u=0,$$

$$\left(\frac{{\partial }^2}{\partial z^2}+{\nabla }_T^2\right)u+k^{2}u=0,$$

$${\nabla }_T^2:=\left(\frac{{\partial }^2}{\partial x^2}\right)+\left(\frac{{\partial }^2}{\partial y^2}\right)$$

$$u=\psi e^{i{k}_{0}z}.$$

$$\left(\underset{\approx 0}{\underbrace{\frac{{\partial }^2\psi }{\partial z^2}}}+2i{k}_0\frac{\partial \psi }{\partial z}-k_0^2\psi +{\nabla }_T^2\psi \right)+k^2\psi =0.$$

$$\left|\frac{{\partial }^2\psi }{\partial z^2}\right|\ll \left|k_0\frac{\partial \psi }{\partial z}\right|.$$

$$2i{k}_0\frac{\partial \psi }{\partial z}=-{\nabla }_T^2\psi +\left(k_0^2-k^2\right)\psi =0.$$

$$\boxed{\begin{array}[]{lrl}\hbox{Paraxial Optics:}&-\frac{1}{i}\frac{\partial\psi}{\partial z}&=-\frac{1}{2k_{0}}\nabla_{T}^{2}\psi+V\psi,\\ &\\ \hbox{Quantum Mechanics:}&-\frac{\hbar}{i}\frac{\partial\psi}{\partial t}&=-\frac{\hbar^{2}}{2m}\nabla^{2}\psi+V\psi,\\ \end{array}}$$

$$V=\frac{1}{2}\left(k_0^2-k^2\right)=\frac{1}{2}{k}_0^2\left(1-n^2\right),$$

$$\begin{array}[]{rll}\hbox{QM}&}{&\hbox{PO}\\ \hline\cr\hbar&\longrightarrow&1\\ m&\longrightarrow&k_{0}\\ t&\longrightarrow&z\\ \end{array}$$

$$-\frac{\hslash }{i}\frac{\partial }{\partial t}=H\psi =(T+V)\psi ,$$

$$T:=-\frac{{\hslash }^2}{2m}{\nabla }^2,$$

$$\psi (r,t)=\mathcal{U}\left(t,t_0\right)\psi (r,t_0)⟹\mathcal{U}\left(t,t_0\right)=e^{-\frac{i}{\hslash }{\int }_{t_0}^{t}Hdt}.$$

$$e^{A+B}=e^A{e}^B{e}^{-[A,B]/2},$$

$$e^{H(t-t_0)}=e^{A/2+B+A/2}\approx e^{A/2}{e}^B{e}^{A/2},$$

$$A$$

$$B$$

$$\psi \left(r,t_0+\Delta t\right)$$

$$\approx$$

$$n(x,y,x)=n_0+\Delta n\cos \left((k_s-k_r)\cdot \rho +\ell \varphi \right),$$

$$u(x,y,z_0,t_0)=\psi (x,y,z_0,t_0)e^{i{k}_{0}z}=e^{-{\rho }^2/{\sigma }^2+i\ell \varphi }{e}^{i{k}_{0}z}.$$

$$d{s}^2={\Omega }^2(\mathbf{r})c_0^{2}d{t}^2-{\Phi }^{-2}(\mathbf{r})d{\mathbf{r}}^2,$$

$$c(\mathbf{r})=|d\mathbf{r}/dt|=c_0\Phi (\mathbf{r})\Omega (\mathbf{r}).$$

$$n(\mathbf{r})=\Phi (\mathbf{r}{)}^{-1}\Omega (\mathbf{r}{)}^{-1}.$$

$$d^2\mathbf{r}/d{A}^2=\nabla \left(\frac{1}{2}n(\mathbf{r}{)}^2{c}_0^2\right)$$

$$d^2\mathbf{r}/d{A}^2=\nabla \left(\frac{1}{2}{n}^2(\mathbf{r})V^2(\mathbf{r})\right)$$

Full text

1. INTRODUCTION

We are interested in numerical tools to efficiently and accurately simulate a photon wavefront propagating in a mildly curved ($M/r\ll 1$ ) spacetime geometry for the secure distribution of a one time only key from a sender Alice ($A$) to a receiver Bob ($B$), as illustrated in Fig. 1. To secure such a key requires a quantum key distribution (QKD) system, and involves state preparation, state propagation, and state detection. This QKD scenario has been exhaustively studied in the literature and is replete with security proofs for numerous protocols, e.g. [5, 9, 8]. These security proofs have been extended in many cases to higher-dimensional state spaces [11], and all of the protocols have been or are currently being demonstrated successfully on earth, and in free-space links to satellites.[24, 10].

Conventional realizations of QKD today involve transmitting heavily-attenuated laser pulses from Alice to Bob and encoding qubit information in each packet by utilizing the spin of the photon. This allows Alice and Bob, who are suitably authenticated, the possibility to establish and share an arbitrarily-secure one-time only key between them. Here they have access to a two-dimensional state space and can therefore form three MUBs each with two orthogonal polarization states. Such a six-state QKD scheme [6] has limited bandwidth and optical fidelity constraints. These constraints can be ameliorated by extending the QKD to higher-dimensional state space [11].

The potential of extending photon-based QKD to higher dimensions was first introduced in a delayed choice experiment using a photons orbital angular (OAM) in 1983 by Wheeler and Miller, and was popularized in 1992 when Allen et al. showed that Laguerre-Gaussian light beams possessed a quantized orbital angular momentum (OAM) of $l\hbar$ per photon [19, 3]. This opened up an arbitrarily high dimensional quantum space to a single photon [2]. This was further developed by Mair et al. when they showed that pairs of OAM photons can be entangled using parametric down conversion, and then when Molina-Terriza et al. introduced a scheme to prepare photons in multidimensional vector states of OAM [18, 22, 21] . These early discoveries began the approach of OAM QKD.

There are advantages and disadvantages of using the extrinsic wavefront modulation of a photon to communicate quantum information as opposed to the intrinsic polarization states of a photon. The main advantage is that higher-dimensional QKD can, in principle, increase bandwidth and be more tolerant to bit error rate (BER) [11]. It also offers a higher degree of entanglement. However, there are disadvantages as well. First, such modulated photon wavefronts can be degraded in propagation due to variations in index of refraction, where a photon’s polarization is more robust [23, 1] Second, an apertured beam diverges. For OAM photons, this divergence is proportional to the square root of the angular momentum, and substantially larger apertures will be needed. Finally, one needs to efficiently sort OAM photons. Given these difficulties, it has been pointed out that photon wavefronts that are a superposition of a certain set of transverse linear momentum states can be used as a signal for higher-dimensional QKD, and these states can be feasibly and efficiently sorted using volume holographic elements [20].

We focus here on the efficient and accurate simulation of the detection and propagation of single photons with arbitrary complex wavefronts in a curved spacetime geometry. In particular we simulate OAM modes or optical vortices. In Sec. 2 we review the paraxial wave equation and its similarity to the Schrödinger equation. In this section, we also discuss the split operator method to solve the paraxial wave equation. This is the numerical approach we will use in our simulations. It is efficient and accurate as it utilizes the fast Fourier transform (FFT). In Sec. 3 we simulate an optical vortex sorter using our split operator code. It involves the propagation of a wavefront through a isotropic but inhomogeneous optical medium. In Sec. 4 we review an optical-mechanical analogue of a conformally flat spacetime. In this way we can simulate the propagation of our wavefronts in a curved spacetime through an equivalent optical medium with varying index of refraction. In Sec. 2 we solve the propagation of a modulated wavefront in a gravitational field using the optical-mechanical analogue. The idea of simulating a gravitational field by an effective index of refraction has been around for a long time,[7] and have recently come to the fore in quantum optics and quantum computing technologies.We solve the paraxial wave equation for this state for an equivalent index of refraction, and in Sec. 7 we discuss the results.

2. THE PARAXIAL WAVE EQUATION IN FLAT SPACETIME AND THE OPERATOR SPLITTING METHOD

The vector equation for the electric field in a flat spacetime in some inertial observers frame is obtained by the covariant Maxwell’s equations where,

[TABLE]

Here, we have assumed that the source-free medium is linear and isotropic, and that it is non magnetic, $\mu=\mu_{0}=const.$ However the medium may be inhomogeneous $\epsilon=\epsilon(r)$, as well as the index of refraction,

[TABLE]

where we use geometric units where the speed of light, $c=1$. Separating variables in Eq. 1 with

[TABLE]

we obtain the usual modified Helmholtz equation for the spatial electric field,

[TABLE]

where $k_{0}=\omega_{0}/c$ is the vacuum wave number and $\omega_{0}$ is the vacuum angular frequency constant from the separation of variables,

[TABLE]

The under-braced term in Eq. 4 is ordinarily non-zero for an inhomogeneous medium. However, we can neglect this term for our purposes since our apertures are large with respect to our wavelength, and since we are dealing with slowly-varying and small index of refraction variations in a thick dielectric diffraction grating (volume hologram). If we neglect this term then we obtain the usual Helmholtz, or scalar wave equation,

[TABLE]

where the wave number $k=nk_{0}$ depends on the inhomogeneous isotropic index of refraction of the medium we are propagating through. This equation can be solved efficiently and accurately if we make one more assumption that applied to our analysis in this manuscript. In particular, we make the slowly-varying envelope approximation, where we assume the electric field is propagating along the $z$ axis in our observers laboratory frame and the remaining field amplitude is a slowly varying function, $\psi(x,y,z)$. Here we rewritie Eq. 6

[TABLE]

where $\nabla_{T}$ is the transverse Laplacian,

[TABLE]

and

[TABLE]

Inserting, Eq. 8 into Eq. 7 yields

[TABLE]

At this point, we assume that the wave vector is close to the optical axis, $z$, so

[TABLE]

This is referred to as the paraxial approximation in wave optics [25], resulting in a parabolic equation,

[TABLE]

This paraxial equation for optics is a 2-dimensional Schrödingers equation for a quantum particle in a potential $V$,

[TABLE]

where the effective optical potential depends on the index of refraction of the medium,

[TABLE]

and we make the identifications,

[TABLE]

Fortunately, there are highly efficient numerical approaches to solve time-dependent Schrödinger-type equations. We know of no better numerical algorithm for our purposes than the split operator method [16, 15, 17]. This method was developed initially to solve the paraxial optics equation, and later it was applied to the Schrödinger equation. This particular approach uses the fast Fourier transform (FFT) , is second or higher order (or higher) convergent, and preserves normalization throughout the evolution. We briefly outline this algorithm.

We find it convenient to express the Hamiltonian operator in Eq. 12 in terms of the sum of two operators,

[TABLE]

where the momentum or “slip operator,”

[TABLE]

and the “kick operator,” $V$. We seek a unitary transformation operator that is a solution of Eq. 15, therefore

[TABLE]

This evolution operator can be subtile when there is time dependence (we will not consider this here), and care must be taken because of non-commutitivity of exponentiated operators. One must utilize the Campbell-Baker-Housdorff identity where for two operators $A$ and $B$,

[TABLE]

is exact only if $[A,[A,B]]=[B,[A,B]]=0$. One can use this identity to remove lower-order error terms in the unitary evolution operator to obtain a second-order accurate operator,

[TABLE]

. with

[TABLE]

If we consider a small time step where $t\rightarrow t_{0}+\delta t$ then the evolution from Eqs.17-21, then

[TABLE]

These “kick” and “slip” operators can be determined numerically by a succession of FFT, and inverse FFT in the following five steps:

1. (1)

Multiply to define, $\psi_{1}(r,t_{0})=e^{-i\frac{V}{2\hbar}\Delta t}\psi(r,t_{0})$, 2. (2)

Fourier transform to get $\phi_{1}(p,t_{0})={\mathcal{F}}\left(\psi_{1}(r,t_{0})\right)$, 3. (3)

Multiply to generate $\phi_{2}(p,t_{0})=e^{-i\frac{T(p)}{\hbar}\Delta t}\phi_{1}(p,t_{0})$, 4. (4)

Inverse Fourier transform to obtain $\psi_{2}(r,t_{0})={\mathcal{F}}^{-1}\left(\phi_{1}(p,t_{0})\right)$,  finally 5. (5)

Multiply to define, $\psi(r,t_{0}+\Delta t)=e^{-i\frac{V}{2\hbar}\Delta t}\psi_{2}(r,t_{0})$.

One need not worry about the normalization factors that accompany the FFT since there is always an even number of FFT’a and inverse FFT’s in the split operator algorithm. By making use of the transformations of Eq. 14, this split operator method applies equally well to the paraxial wave equation and its evolution, and this optical application is the primary focus of this manuscript.

Before we consider optical wavefront propagation in a curved spacetime geometry, we will use the split operator method to analyze the propagation of an optical vortex through a volume hologram with inhomogeneous index of refraction in the next section. And in particular we will numerically calculate the efficiency of a volume-holographic optical-vortex sorter.

3. THIS VORTEX THAT VORTEX: AN OPTICAL VORTEX SORTER

It is sufficient here to thoroughly analyze a single optical vortex grating and to exhibit its ability to differentiate a photon with one value of OAM from another. Once this has been demonstrated, the higher-dimensional sorters can be constructed by literally stacking or sequencing such gratings, one from each basis state in a given MUB. Another approach will be incoherently multiplex each of the gratings into a single holographic substrate. In this section,we examine a holographic emulsion produced by the coherent interference of a unit-amplitude plane wave and the unit amplitude OAM wave front corresponding to $\ell$. Specifically, the hologram’s face is oriented so that its normal is along the $z$-axis of our observer. The hologram’s face or aperture is rectangular with dimensions $L_{x}$ by $L_{y}$ and the thickness of the hologram is $L_{z}$. The bulk index of refraction is $n_{0}$, but this is modulated by an interference pattern of amplitude $\Delta n$. The hologram is generated by the interference pattern of two waves of angular frequency $\omega_{0}$, (1) an incident reference plane wave whose wave vector, $k_{r}$, is in the $x$-$z$ plane and tilted at an angle $\theta_{r}$ from the $z$-axis, and (2) a modulated signal plane wave of wave vector, $k_{s}$, also in the $x$-$z$ plane but titled an angle $\theta_{s}$ from the $z$-axis. The modulation is an optical vortex with orbital angular momentum $\ell$. The interference of these two waves gives rise to a modulation in the index of refraction in the volume hologram,

[TABLE]

where, $\rho^{2}=x^{2}+y^{2}$ is the radius and $\phi$ is the azimuthal angle in the observer’s cylindrical coordinates, $\{\rho,\phi,z\}$. This is illustrated in Fig. 2.

We can examine the sorting efficiency of this hologram by (1) sending in an identical signal wave with a Gaussian envelope into the hologram and measuring how much of the signal is directed toward the detector in the far-field behind the hologram at an angle $\theta_{r}$ in the $x$-$z$-plane, and compare this with (2) the signal intensity for a beam with a miss-matched value for its orbital angular momentum $\ell_{s}\neq\ell$.

The initial amplitude of the optical wavefront in Eq. 8 i is a Gaussian beam of width $\sigma$ that is modulated by an optical vortex with orbital angular momentum, $\ell$,

[TABLE]

4. EFFECTIVE INDEX OF REFRACTION FOR A BLACK HOLE SPACETIME

Metrics in general relativity (GR) that can be written in isotropic coordinates lend themselves to a simple exact formulation such that the effects of the curved spacetime are encapsulated in a spatially varying index of refraction $n(\mathbf{r})$ [12]. A metric in isotropic coordinates has the form

[TABLE]

where $\mathbf{r}=(x,y,z)=\{x^{i}\}_{i\in\{1,2,3\}}$ are Cartesian coordinates. Thus, a metric in isotropic coordinates is conformally flat. The isotropic coordinate speed of light $c(\mathbf{r})$ is determined by the condition that the geodesics are null ($ds=0)$ leading to

[TABLE]

Thus the effective index of refraction for light in the gravitational field is then

[TABLE]

The index of refraction may be used to formulate a Newtonian-like “F=ma” formulation of geometric optics in which the equation for the path of the ray assumes the form of a Newton’s law of motion [14, 13]

[TABLE]

where $dA=dt/n(\mathbf{r})^{2}$ is the so called optical action [12]. The lefthand side of this equation has the form of an acceleration, while the righthand side has the form of a force, the gradient of a potential-energy function $U(\mathbf{r})=-(\frac{1}{2}\,n(\mathbf{r})^{2}\,c_{0}^{2})$. Using the optical action, one can also show that motion of massive particles is given by a similar equation,

[TABLE]

where

[TABLE]

Thus, in the massless limit, $V\rightarrow c_{0}$ this formula reduces to the one for light. Thus the latter formulas holds for both massive and massless particles. While in this work we will be primarily interested in massless particles (light), in this section we will carry out a single formalism that is valid for both massive and massless particles.

The optical-mechanical analogy states that for light, the normal to surfaces of constant phase is the wavevector, while for massive particles, the normal to surfaces of constant action is the particle’s momentum. This can put on a firm footing by considering a variational principle for the particle’s geodesics via

[TABLE]

using the line element $ds$ which contains the metric coefficients defining the spacetime. This can be cast into a Lagrangian form

[TABLE]

where

[TABLE]

Here, $V^{2}=\sum_{i=1}^{3}(V^{i})^{2}$ with $V^{i}=dx^{i}/dt$ the components of the coordinate velocity. The canonical momenta are given by

[TABLE]

and the effective Hamiltonian (energy) when expressed in terms of the momenta is given by

[TABLE]

which is constant of the motion along the geodesic. Lastly, the magnitude of the coordinate velocity is given by

[TABLE]

which generalizes the usual phase velocity for light $v=c_{0}/n$.

To bring in quantum mechanics, we follow [12] de Broglie and define the constant coordinate frequency $\omega$ via $H\equiv\hbar\,\omega$. The measurable proper frequency $\tilde{\omega}$ along the particles geodesic is given by using the isotropic metric to obtain $\tilde{\omega}(\mathbf{r})=\omega/\Omega(\mathbf{r})$. By using the de Broglie relationship $p(\mathbf{r})\equiv\hbar\,k(\mathbf{r})$ we obtain the coordinate wavelength $\lambda=\hbar\,c_{0}^{2}/(n^{2}(\mathbf{r})\,H\,V(\mathbf{r}))$ with measurable proper wavelength $\tilde{\lambda}(\mathbf{r})=\lambda/\Phi(\mathbf{r})$. One can then define an index of refraction $N(\mathbf{r})$ for both massive and massless particles that satisfies $\lambda(\mathbf{r})\,N(\mathbf{r})=$ constant with $N(\mathbf{r})\equiv n\,[1-m^{2}\,c_{0}^{4}\,\Omega^{2}(\mathbf{r})/H^{2}]^{1/2}$, which reduces to $N(\mathbf{r})\rightarrow n(\mathbf{r})$ for light by setting $m=0$. We obtain gravitation redshifts by using $N(\mathbf{r}_{1})\,\lambda(\mathbf{r}_{1})=N(\mathbf{r}_{2})\,\lambda(\mathbf{r}_{2})$. Converting to proper wavelengths we have $N(\mathbf{r}_{1})\,\tilde{\lambda}(\mathbf{r}_{1})\,\Phi(\mathbf{r}_{1})=N(\mathbf{r}_{2})\,\tilde{\lambda}(\mathbf{r}_{2})\,\Phi(\mathbf{r}_{2})$. We can define the phase velocity $v_{p}=\omega/k=H/p$, or equivalently $v_{p}(\mathbf{r})=c_{0}/N(\mathbf{r})$, which shows that $N(\mathbf{r})$ is playing the proper role as an index of refraction. By substituting the de Broglie definitions into the Hamiltonian we obtain $\hbar^{2}\,\omega^{2}=m^{2}\,c_{0}^{2}\,\Omega^{2}(\mathbf{r})+c_{0}^{2}\,\hbar^{2}\,k^{2}(\mathbf{r})/n^{2}(\mathbf{r})$. Upon differentiating this expression with respect to $k$ we obtain the group velocity $v_{p}(\mathbf{r})=d\omega/dk=c_{0}N(\mathbf{r})/n^{2}(\mathbf{r})$, and with a little algebra this yields $v_{g}(\mathbf{r})=V(\mathbf{r})$, the magnitude of the coordinate velocity.

Finally, we can first quantize the above Hamiltonian by converting it to a relativistic Klein-Gordon equation for the wavefunction $\psi(x)$ for a spinless massive particle with the substitution of $p_{\mu}\rightarrow-i\,\hbar\,\nabla_{\mu}$ into the curved-spacetime mass-shell constraint

[TABLE]

where $\nabla_{\mu}$ is the usual covariant derivative,

[TABLE]

with Christoffel connection

[TABLE]

Thus, the Klein Gordon equation

[TABLE]

becomes

[TABLE]

In a static isotropic metric $g_{\mu\nu}(x)=\textrm{diag}(\Omega^{2}\,c_{0}^{2},-\Phi^{-2},-\Phi^{-2},-\Phi^{-2})$, and with the aid of the identity $g^{\mu\nu}\,\Gamma^{\lambda}_{\mu\,\nu}=-(-g)^{-1/2}\,\partial_{\mu}[(-g)^{1/2}\,g^{\mu\lambda}]$ we obtain

[TABLE]

where $\boldsymbol{\nabla}$ is the Cartesian divergence and $\boldsymbol{\nabla}^{2}$ is the Cartesian Laplacian. In Eq. 43 we have defined $\xi(\mathbf{r})=\ln(\Omega(\mathbf{r})/\Phi(\mathbf{r}))$ and $k_{c}^{2}(\mathbf{r})=m\,c_{0}/(\hbar\,\Phi(\mathbf{r}))$ which can be interpreted as the particle’s coordinate Compton wavelength (the particle’s invariant proper Compton wavelength is $h/(m\,c_{0})$). We can remove the cross term $\boldsymbol{\nabla}\xi\cdot\boldsymbol{\nabla}\psi$ arising from the covariant derivative by letting $\psi(\mathbf{r})\equiv f(\mathbf{r})\,\phi(\mathbf{r},t)$ where $f(\mathbf{r})=(\Omega/\Phi)^{1/2}\equiv e^{-\xi/2}$. Then $\phi(\mathbf{r},t)$ satisfies the wave equation $(n^{2}/c_{0}^{2})\,\partial_{t}^{2}\phi-\boldsymbol{\nabla}^{2}\phi+[k_{c}^{2}(\mathbf{r})+\eta(\mathbf{r})]=0$ where $\eta(\mathbf{r})=\frac{1}{2}\,\nabla^{2}\xi(\mathbf{r})+\frac{1}{4}\,|\nabla\xi(\mathbf{r})|^{2}$. Finally, letting $\phi(\mathbf{r},t)\equiv e^{i\,H\,t/\hbar}u(\mathbf{r})$ we arrive at the Helmholz equation valid for both massive and massless particles,

[TABLE]

where we have used $H\equiv\hbar\,\omega_{0}$ and $k_{0}=\omega_{0}/c_{0}=H/(\hbar\,c_{0})$.

For the Schwarzschild metric in isotropic coordinates we have

[TABLE]

where $H^{{}^{\prime}}\equiv H/(m\,c_{0}^{2})$ and the radial coordinate is given by $\rho=r/r_{s}$ with $r_{s}=G\,M/c_{0}^{2}$ (half the Schwarzschild radius $r_{S}=2\,G\,M/c_{0}^{2}$). The constant $H^{{}^{\prime}}$ can be evaluated at any point on the particle’s geodesic where it’s coordinate velocity is known.

It can be shown [12] that $\eta(\mathbf{r})\ll k^{2}(\mathbf{r})$ for all values of $1\leq\rho\leq\infty$ (the valid range of the isotropic scaled radius $\rho$), so that we can drop the former. For massless particles $H^{\prime}\rightarrow\infty$ ($m\rightarrow 0$) and $k=k_{0}\,n$, so that for light we end up with the Helmholz equation for a scalar field with a spatially varying wavevector $k(\mathbf{r})$,

[TABLE]

The above derivations reveal that the complete effects of the gravitational field are now encapsulated in the spatially varying index of refraction $n(\mathbf{r})$ (for light) given by the combination of the metric coefficients $(\Omega(\mathbf{r})\,\Phi(\mathbf{r}))^{-1}$. Thus, for example, the well known bending of light around a massive object can now be seen as the usual bending of light as the photon moves into a more intense index of refraction [12, 4]. The fact the matter waves bend more than light is now attributed to the latter’s slower velocity. The gravitational wave optics of massive and massless particles differ only in a dispersion effect, since the index of refraction for massive objects $N(\mathbf{r})$ depends on the parameter $m\,c_{0}^{2}/H$.

In the next section we will explore the paraxial approximation to Eq.(46) for light and explore the propagation of wavefronts in a gravitational field described by the static Schwarzschild metric.

5. PARAXIAL WAVE EQUATION IN A CURVED SPACETIME

In this section we will assume that both Alice and Bob are at rest in the isotropic coordinates $\{t,x,y,z\}$ of our conformally-flat spacetime of Eq. 26. This assumption we make is not out of necessity, but will avoid complications of boosting into Alice and Bob’s frame at event $A$ and $B$ shown in Fig. 1. The photon trajectory connecting $A$ to $B$ is determined by the wave 4-vector,

[TABLE]

and in the optical-mechanical model we are following in Evans et al. the spatial projection of this null trajectory onto the $t=const$ spacelike hypersurface of Allice is given by their force equation,

[TABLE]

where $dA=dt/n^{2}$ and $\nabla={}^{\!(3)}\nabla$ is the 3-dimensional gradient operator. This trajectory is a spacelike curve ${\mathcal{C}}$ and is illustrated in Fig. 4 by the solid black line connecting $A$ to $B_{\perp}$.

Fortunately , the optical-mechanical Eqns. 46-48 already take into account general relativistic curvature effects within the geometrical optics approximation, and provides us with the Helmholtz equation in the spacelike hypersurface of $A$ and it is written in the Cartesian coordinates of Alice’s laboratory. The derivation of the paraxial wave equation, (Eq. 12) for optics from the Helmholtz equation (Eq. 7) in Sec. 2 can be applied to this optical-mechanical Helmholtz equation (Eq. 46) with one modification. In flat spacetime the a plane-wave with wave number $k_{0}$ will propagate in a straight line in the laboratory with the same frequency. In a curved spacetime the photon will ordinarily be redshifted and its trajectory with tangent 3-vector, $k$, will ordinarily be curved in the laboratory 3-space. For these reasons we would like to prove that we can transform the Helmoltz given in Eq. 46 at some point $P$ in in the $t=const$ spacelike hypersurface (e.g. Alice and Bob’s laboratory) as

[TABLE]

and $\overline{\nabla}_{T}$ is the transverse Laplacian in the barred coordinates,

[TABLE]

This “barred” frame is Fermi-Walker transported along the photon trajectory in the 3-space connecting $A$ to $B_{\perp}$ as illustrated in Fig. 4.

Here $\bar{z}$ is defined at each point $P$ along ${\mathcal{C}}$ to be the proper spacelike distance along ${\mathcal{C}}$ from $A$ to $P$,

[TABLE]

and the $\bar{z}$-axis, $e_{\bar{z}}=\partial_{\bar{z}}$, is the tangent to the curve ${\mathcal{C}}$ at each $P$ and is parallel to the spatial projection of the wave vector,

[TABLE]

the transverse axes $e_{\bar{x}}$ and $e_{\bar{y}}$ are Fermi-Walker transported along ${\mathcal{C}}$,

[TABLE]

Whereas $\bar{z}$ measures the distance along $\mathcal{C}$ from $A$ to $P$, $\bar{x}$ and $\bar{y}$ measures the proper distance perpendicularly away from $\mathcal{C}$ along $e_{\bar{x}}$ and $e_{\bar{y}}$; respectively. The spatial part of $k$ is tangent to the photons trajectory in this photon trajectory frame. This as illustrated in in Fig. 4 by the coordinate axes at event $P$. A necessary condition for the derivation of the paraxial equation in this hypersurface is that we solve Eq. 48 for the curve ${\mathcal{C}}$, and we have the Jacobian representing the rotational operator in space connecting the bared frame to the metric frame

[TABLE]

and its inverse.

Paralleling the development of the paraxial equation in Sec. 2, we assume the field is the product of a slowly varying amplitude $\psi$ that serves to modulate a plane wave train traveling tangent to the wave vector. In our barred coordinate frame on ${\mathcal{C}}$, we let

[TABLE]

Here $u$, $\psi$ and $k$ all are functions of $\bar{x}$, $\bar{y}$ and $\bar{z}$. We begin to transform Eq. 46 by Eq. 53,

[TABLE]

Inserting, Eq. 54 into Eq. 46 will yield the Helmholtz equation in the the barred coordinates. It is rather involved and we will not show it here. Furthermore, this equation will not easily be put into the form we assumed in Eq. 49. In particular, we have mixed derivatives of $psi$ for an arbitrary transformation matrix $\Omega$. We believe we can justify approximations; however, this is beyond the scope of this manuscript and beyond the week-field application we intended. Nevertheless, for weak gravitational fields, e.g. around the Earth, we can work in the isotropic coordinates without going to the bared coordinates. That is we expect minimal bending of the photon trajectories and the paraxial approximation about the $z$-axis should be valid to good approximation. We then assume the transformation matrix is the identity matrix,

[TABLE]

We then rewrite Eq. 46 in the unbarred coordinates by singling out the $z$-axis for propagation,

[TABLE]

and $\nabla_{T}$ is the transverse Laplacian in the barred coordinates,

[TABLE]

Again we parallel the development of the paraxial equation in Sec. 2 but this time for the coordinates $r=\{x,y,z\}$ of the $t=constant$ spacelike hypersurface of our conformally-flat spacetime. We again assume that the field is the product of a slowly varying amplitude $\psi$ that serves to modulate a plane wave train traveling tangent to the wave vector. In our global coordinates we again let

[TABLE]

Here $u$, $\psi$ and $k$ all are functions of $x$, $y$ and $z$. Since we are in isotropic flat coordinates, the $z$-direction is just as good as any direction we could have chosen. Inserting Eq. 8 into Eq. 7, we find

[TABLE]

We would like to prove that we can ignore the three terms in the second line of this equation, and we think that this is reasonable. However, we have not been able to provide a proof. We continue with this assumption for the remainder of this manuscript. We will address the proof in a future paper. Therefore given these assumptions, we obtain the optical-mechanical equivalent of the paraxial equation for a beam propagating primarily down the $z$-axis in our spacetime coordinates,

[TABLE]

where the effective wave number is now

[TABLE]

and the effective potential is

[TABLE]

This paraxial equation still predicts the bending of the photon trajectory and its redshift. Connecting with the previous section we generalized Eq. 8 and wrote,

[TABLE]

where, the surface of constant phase

[TABLE]

has normals in the direction of $\nabla S_{0}(r)=k(r)$ at each point $P$ along the photon’s 3-trajectory connecting Alice to Bob. The wave 3-vector threads crossed each of these phase surfaces at normal angles, and this defines the trajectory of our photon in the 3-space from Alice to Bob as illustrated in the left of Fig. 4. We will further this calculation for stronger gravitational bending effects by using the Fermi-Walker transported coordinates for Alice’s lab in order to construct the paraxial equation from the Helmholtz equation. This will be discussed in a future manuscript.

We are now in a position to simulate a photon’s wavefront evolution in the optical-mechanical analogue using our split operator code to solve Eq. 61. We provide a simple numerical example of the solution in the next section.

6. SIMULATION OF AN OPTICAL VORTEX IN A CURVED SPACETIME

We simulate an optical vortex given by Eq. 25 propagating outward along the $z$-axis in our isotropic coordinates of a Schwarzschild spacetime of mass $M$ given in conformally-flat coordinates. We start the simulation at $z=z_{0}$ and examine the the evolution using our split operator code. We numerically solve Eqs. 61-63.

In order to determine the effective wave number $\kappa$ and effective potential $\mathcal{V}$, we use the wave vector in Eq 46 and Eq. 4. We find,

[TABLE]

and the effective potential

[TABLE]

The effective wave number and the effective potential are graphed in Figs. 5-6.

We used the split operator code to propagate the $\ell=10$ optical vortex with a Gaussian envelope that we used in Sec. 2. Here we set $M=1$,, $k_{0}=2\pi$ and started the wavefront at $z_{0}=25$ and propagating inward toward the origin. We illustrate in Fig. 7 the initial wavefront intensity and the evolved wavefront showing the development of the orbital angular momentum ($\ell=10$) ring.

7. CONCLUSION

In this manuscript we reviewed the paraxial approximation and we outlined the split operator method to solve the paraxial equation. We applied our numerical approach in flat Minkowski spacetime to numerically simulate a volume holographic sorter by propagating both a matched and mismatched optical vortex through the glass. Using this example, we demonstrated the efficiency, accuracy and convergence of the code. We then numerically examined an optical-mechanical analogue model of the propagation of a photon in a curved spacetime. We derived a curved-space paraxial wave equation for this model that was derived from the optical-mechanical Helmholtz equation. While we have not yet formulated an paraxial evolution equation using the Fermi-Walker coordinates that track the trajectory of the photon, we are nevertheless able to construct a Schrödinger-like equation for a weak gravitational fields in the isotropic coordinates. More work is required in order to prove that two of the three terms we neglected are small compared to the other leading-order terms. We nevertheless, used the curved-spacetime paraxial equation to simulate an optical vortex propagating in a Schwarzschild spacetime that was expressed in an isotropic conformally-flat coordinate system. We assumed that our two observers were co-moving synchronous observers to avoid further Lorentz boosts into an arbitrary rest frame, although this can be easily added by a Lorentz boost of our results.

We look forward to reformulating this curved-spacetime paraxial approximation to simulate the scattering of wavefronts by a black hole at moderate impact parameters. This would require us to address the paraxial approximation in the Fermi-Walker coordinates where the wave vector is tangent to the world line of the photon in the spacelike hypersurface.

ACKNOWLEDGMENTS

PMA and WAM would like to acknowledge support of the Air Force Office of Scientific Research (AFOSR). This work was supported under joint iC2S2 grants from the Korean Ministry of Science and Future Planning (MSIP) IITP 2017-0-00266 and the US Air Force Asian Office of Aerospace Research & Development (AOARD) FA2386-17-1-4070. Any opinions, findings, conclusions, or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of AFRL.