Evolution equation for multi-photon states in turbulence

TL;DR

This paper develops a Wigner functional-based evolution equation for multi-photon quantum states propagating through turbulent atmospheres, unifying different approaches and enabling new solution methods.

Contribution

It introduces a Fokker-Planck type evolution equation for multi-photon states in turbulence using Wigner functional theory, extending previous methods.

Findings

The evolution equation is consistent with prior solutions.

The equation is in the form of a Fokker-Planck equation with functional derivatives.

Potential methods for finding additional solutions are discussed.

Abstract

The recently developed Wigner functional theory is used to formulate an evolution equation for arbitrary multi-photon states, propagating through a turbulent atmosphere under arbitrary conditions. The resulting evolution equation, which is obtained from an infinitesimal propagation approach, is in the form of a Fokker-Planck equation for the Wigner functional of the state and therefore incorporates functional derivatives. We show consistency with previously obtained solutions from different approaches and consider possible ways to find additional solutions for this equation.

Equations

$${\nabla }_T^{2}g(\mathbf{x},z)-\text{rmi}2k{\partial }_{z}g(\mathbf{x},z)+2k^2\tilde{n}(\mathbf{x},z)g(\mathbf{x},z)=0.$$

$${\partial }_{z}G(\mathbf{a},z)=\frac{\text{rmi}2{\pi }^2}{k}|\mathbf{a}{|}^{2}G(\mathbf{a},z)-\text{rmi}k\int N(\mathbf{a}-{\mathbf{a}}^{\prime },z)G({\mathbf{a}}^{\prime },z){\text{rmd}}^2{a}^{\prime },$$

$${\partial }_{z}G(\mathbf{a},z)=\text{rmi}\int M(\mathbf{a},{\mathbf{a}}^{\prime },z)G({\mathbf{a}}^{\prime },z){\text{rmd}}^2{a}^{\prime }=\text{rmi}M(z)\diamond G(z),$$

$$M(\mathbf{a},{\mathbf{a}}^{\prime },z)\equiv \pi \lambda |\mathbf{a}{|}^2\delta (\mathbf{a}-{\mathbf{a}}^{\prime })-kN(\mathbf{a}-{\mathbf{a}}^{\prime },z),$$

$$G(\mathbf{a},z+dz)=\int \left[\delta (\mathbf{a}-{\mathbf{a}}^{\prime })+\text{rmi}dzM(\mathbf{a},{\mathbf{a}}^{\prime },z)\right]G({\mathbf{a}}^{\prime },z){\text{rmd}}^2{a}^{\prime }.$$

$$|G(z)⟩=\int |\mathbf{a}⟩G(\mathbf{a},z){\text{rmd}}^{2}a.$$

$$1\equiv \int |\mathbf{a}⟩⟨\mathbf{a}|{\text{rmd}}^{2}a,$$

$$\hat{M}(z)\equiv \int |\mathbf{a}⟩M(\mathbf{a},{\mathbf{a}}^{\prime },z)⟨{\mathbf{a}}^{\prime }|{\text{rmd}}^{2}a{\text{rmd}}^2{a}^{\prime },$$

$$\hat{U}(z,dz)$$

$$\text{mathds}1=\sum _n\frac{1^{\otimes n}}{n!}={\exp }_{\otimes }(1)$$

$$|\psi (z+dz)⟩=\hat{U}(z,dz)|\psi (z)⟩=|\psi (z)⟩+\text{rmi}dz\hat{\mathcal{M}}(z)|\psi (z)⟩.$$

$$\begin{gathered} \text{eqalign}{\partial }_z|\psi (z)⟩=\text{rmi}\hat{\mathcal{M}}(z)|\psi (z)⟩, \\ {\partial }_z⟨\psi (z)|=-\text{rmi}⟨\psi (z)|\hat{\mathcal{M}}(z). \end{gathered}$$

$${\partial }_z\hat{\rho }(z)=\text{rmi}\left[\hat{\mathcal{M}}(z),\hat{\rho }(z)\right].$$

$$\hat{\mathcal{R}}(z;\eta )={\exp }_{\otimes }\left[\int |\mathbf{a}⟩T(\mathbf{a},{\mathbf{a}}^{\prime },z;\eta )⟨{\mathbf{a}}^{\prime }|{\text{rmd}}^{2}a{\text{rmd}}^2{a}^{\prime }\right],$$

$$T(\mathbf{a},{\mathbf{a}}^{\prime },z;\eta )\equiv \exp \left(\text{rmi}\pi \eta \lambda |\mathbf{a}{|}^2\right)\delta (\mathbf{a}-{\mathbf{a}}^{\prime })-\text{rmi}\eta kN(\mathbf{a}-{\mathbf{a}}^{\prime },z).$$

$${\partial }_{\eta }{\hat{\mathcal{R}}(z;\eta )|}_{\eta =0}=\text{rmi}\hat{\mathcal{M}}(z).$$

$${\partial }_z\hat{\rho }(z)={\partial }_{\eta }{\left[\hat{\mathcal{R}}(z;\eta )\hat{\rho }(z){\hat{\mathcal{R}}}^{†}(z;\eta )\right]}_{\eta =0}.$$

$$W_{\hat{\mathcal{R}}}[\alpha ]=\frac{1}{\det \{1+T\}}\exp \left[-2{\alpha }^{\ast }\diamond (1-T)\diamond (1+T{)}^{-1}\diamond \alpha \right],$$

$${\partial }_z{W}_{\hat{\rho }}(z)[\alpha ]$$

$${\partial }_z{W}_{\hat{\rho }}(z)[\alpha ]=$$

$$\begin{gathered} \text{eqalign}{\alpha }_b\to {\left(1-\frac{\text{rmi}\xi }{2}M\right)}^{-1}\diamond {\alpha }_c, \\ {\alpha }_b^{\ast }\to {\alpha }_c^{\ast }\diamond {\left(1+\frac{\text{rmi}\xi }{2}M\right)}^{-1}. \end{gathered}$$

$${\partial }_z{W}_{\hat{\rho }}(z)[\alpha ]={{\partial }_{\xi }\frac{W_{\hat{\rho }}(z)\left[{\left(2+\text{rmi}\xi M\right)}^{-1}\diamond \left(2-\text{rmi}\xi M\right)\diamond \alpha \right]}{\sqrt{\det \left\{1+\frac{1}{4}{\xi }^2{M}^2\right\}}}|}_{\xi =0}.$$

$${{\partial }_{\xi }{\left(\det \left\{1+\frac{1}{4}{\xi }^2{M}^2\right\}\right)}^{-1/2}|}_{\xi =0}=0.$$

$${\partial }_z{W}_{\hat{\rho }}=\text{rmi}{\alpha }^{\ast }\diamond M\diamond \frac{\delta W_{\hat{\rho }}}{\delta {\alpha }^{\ast }}-\text{rmi}\frac{\delta W_{\hat{\rho }}}{\delta \alpha }\diamond M\diamond \alpha ,$$

$$⟨M⟩=\pi \lambda |\mathbf{a}{|}^2\delta (\mathbf{a}-{\mathbf{a}}^{\prime }).$$

$${\partial }_z{W}_{\hat{\rho }}=\text{rmi}\pi \lambda {\alpha }^{\ast }\diamond |\mathbf{a}{|}^{2}1\diamond \frac{\delta W_{\hat{\rho }}}{\delta {\alpha }^{\ast }}-\text{rmi}\pi \lambda \frac{\delta W_{\hat{\rho }}}{\delta \alpha }\diamond |\mathbf{a}{|}^{2}1\diamond \alpha ,$$

$$W_{\hat{\rho }}(z)=W_{\hat{\rho }}(z_0)+\text{rmi}{\int }_{z_0}^z{\alpha }^{\ast }\diamond M(z_1)\diamond \frac{\delta W_{\hat{\rho }}(z_1)}{\delta {\alpha }^{\ast }}-\frac{\delta W_{\hat{\rho }}(z_1)}{\delta \alpha }\diamond M(z_1)\diamond \alpha \text{rmd}{z}_1.$$

$$W_{\hat{\rho }}(z)=$$

$${\int }_{z_0}^z{\int }_{z_0}^{z_1}⟨M({\mathbf{a}}_1,{\mathbf{a}}_2,z_1)M({\mathbf{a}}_3,{\mathbf{a}}_4,z_2)⟩\text{rmd}{z}_2\text{rmd}{z}_1$$

Full text

Keywords: Free-space propagation, scintillation, multi-photon state, Wigner functional, Fokker-Planck equation

1 Introduction

The effect of a turbulent atmosphere on a photonic quantum state during propagation plays an important role in any free-space quantum information system (FS-QIS), such as free-space quantum communication [1, 2], quantum teleportation over free-space [3, 4, 5], and quantum protocols [6] for free-space time-transfer [7]. The random phase modulations induced by the fluctuating refractive index distort the spatial characteristics of the state, which in turn affects the photon-number characteristics of the state.

Current implementations of such FS-QISs focus on either the spatiotemporal degrees of freedom of a few discrete photons [8, 9, 10, 11, 12, 13, 14, 15, 16], or the particle-number degrees of freedom (so-called continuous variables) for a few discrete spatial modes [17, 18, 19, 20, 21, 22, 23]. The implementations that are based on the spatial degrees of freedom often use orbital angular momentum (OAM) modes to encode the information [24, 25, 26, 27]. The polarization degree of freedom has been incorporated in both trends [24, 28, 29].

The natural culmination of these endeavors would be an approach that incorporates both the spatiotemporal degrees of freedom and the particle-number degrees of freedom to their full extent. For such an implementation, it would be necessary to understand the evolution during propagation through a turbulent atmosphere of a photonic state that is arbitrary in terms of both its spatiotemporal degrees of freedom and particle-number degrees of freedom.

Such an analysis needs a special formalism in which neither the spatiotemporal degrees of freedom, nor the particle-number degrees of freedom, are treated as discrete variables — a few discrete photons or a few discrete modes. While the multi-photon nature of the states is represented by the particle-number degrees of freedom, the turbulent medium affects the spatial degrees of freedom of the state, scattering any input mode into an infinite dimensional modal space. This modal scattering process happens continuously, resulting in optical power being scattered back and forth among an infinite number of modes. A proper analysis of the evolution of multi-photon states, propagating through a turbulent atmosphere, therefore requires a formalism that incorporates both the particle-number degrees of freedom and the spatiotemporal degrees of freedom as continuous degrees of freedom.

Such a special formalism has recently been developed in the form of a Wigner functional formalism [30]. It is a generalization of the well-known Wigner function [31], which has been used extensively (see [32] and references therein). The Wigner functional formalism is based on the spatiotemporal quadrature bases that are eigenstates of the momentum dependent quadrature operators [33]. This formalism is an extension of the Moyal formalism for quantum mechanics [34, 35, 36]. The Wigner functionals are defined in a phase space that represents a Hilbert space for all possible states in quantum optics. The operations are represented by a functional version of the star-product of the Moyal formalism.

The effect of turbulence on photonic quantum states has been studied extensively, but, as with the implementations, these studies invariably only focused on one of the two sets of continuous degrees of freedom. Those studies that focus on the spatial degrees of freedom [37, 8, 9, 11, 38, 12, 14], are mostly based on the single-phase screen approximation [39], which assumes weak scintillation [40]. However, we have also considered a multi-phase screen approach (based on infinitesimal propagation) to derive an evolution equation for the state [41, 42, 43, 44, 45], which can address strong scintillation conditions [46]. Since the scintillation process directly affects the spatial degrees of freedom, the studies that focus on the particle-number degrees of freedom need a method to represent the effect of turbulence on the spatial degrees of freedom. They get around this problem by modeling the effect of turbulence simply as a loss mechanism, ignoring any other role that the spatial degrees of freedom may play in the process [47, 48]. Another approach, which employs a Wigner function for the spatiotemporal degrees of freedom, instead of the particle-number degrees of freedom, represents the effect of turbulence as a random position dependent force [49, 50, 51].

In this paper, we’ll use the Wigner functional formalism for an infinitesimal propagation approach to derive an evolution equation for the Wigner functional of an arbitrary photonic state. The result is a functional differential equation in the form of a Fokker-Planck equation, which we call the infinitesimal propagation functional equation (IPFE). This is the main result of the paper.

The IPFE is difficult to solve. We discuss a few approaches for finding solutions. One of these approaches only provide trivial solutions, the vacuum and thermal states. We show that previous results obtained from considering states with a few discrete photons [44, 45] are also solutions of the IPFE. An interesting observation emerges from considering the characteristic functional of the solution, namely that it obeys an identical functional differential equation as the Wigner functional solution.

2 Infinitesimal propagation

2.1 Unitary evolution

Under general scintillation conditions, a photonic quantum state encounters two effects while propagating through the atmosphere. One is the free-space propagation process (without the effect of turbulence). The other is a continuous phase modulation caused by fluctuations in the refractive index of the medium. While the latter involves a modulation in the position domain, the former represents a modulation in the Fourier domain. The full dynamics, which incorporates both these effects, is captured in the classical context, by the equation of motion for paraxial propagation through a turbulent medium

[TABLE]

Here $g(\mathbf{x},z)$ is the classical scalar electromagnetic field, $k=2\pi/\lambda$ is the wavenumber and $\tilde{n}(\mathbf{x},z)$ is the fluctuation in the refractive index of the atmosphere. The refractive index of the atmosphere can thus be represented as $n=1+\tilde{n}$. In the transverse Fourier domain, the equation becomes

[TABLE]

where $\mathbf{a}$ is the two-dimensional transverse spatial frequency vector, $G(\mathbf{a},z)$ is the transverse angular spectrum, which is also a function of the propagation distance $z$, and $N(\mathbf{a},z)$ is the transverse Fourier transform of the refractive index fluctuation. Moreover, $N^{*}(\mathbf{a},z)=N(-\mathbf{a},z)$, because the refractive index fluctuation is a real-valued function. Unless the field represented by the photonic state is strong enough to induce a nonlinear interaction with the atmosphere, the classical description of the propagation process, as given in (1), suffices to provide the basis for the development of the evolution of photonic states propagating through the atmosphere.

The approach that we’ll follow is to derive a dynamical equation from an infinitesimal propagation of the state, as proposed in [41, 43]. If $\hat{U}(z)$ represents that unitary evolution operator for the propagation of a state over a distance $z$ through the atmosphere and $\hat{\rho}$ is the density operator for the input state, then the output state is given by $\hat{\rho}(z)=\hat{U}(z)\hat{\rho}\hat{U}^{{\dagger}}(z)$. For an infinitesimal propagation, one expects to find that $\hat{U}(dz)=\mathbf{1}+\rmi dz\hat{P}_{\Delta}$, where $\hat{P}_{\Delta}$ represents the infinitesimal propagation. It thus leads to an evolution equation of the form $\partial_{z}\hat{\rho}(z)=\rmi[\hat{P}_{\Delta},\hat{\rho}]$. The infinitesimal propagation operator $\hat{P}_{\Delta}$ assumes full knowledge of the medium that causes the scintillation process. Since we normally only know the statistical properties of the medium, we’ll eventually have to apply an ensemble averaging. However, we start by deriving an expression for the infinitesimal propagation operator, using the classical equation in (2).

To define quantum operators for the scintillation process, it is convenient to combine the two terms in (2), so that

[TABLE]

where the kernel is given by

[TABLE]

and where we adopt the notation introduced in [30] to write the integral contraction as a binary operation denoted by $\diamond$. Note that the kernel is Hermitian $M^{*}(\mathbf{a},\mathbf{a}^{\prime},z)=M(\mathbf{a}^{\prime},\mathbf{a},z)$. Represented as an infinitesimal propagation, the equation in (3) becomes

[TABLE]

The transverse angular spectrum is interpreted as a Fourier domain wave function $G(\mathbf{a},z)=\left\langle{\mathbf{a}}|{G(z)}\right\rangle$, where $\left|{\mathbf{a}}\right\rangle$ is an element of the transverse spatial frequency basis, which obeys an orthogonality condition $\left\langle{\mathbf{a}}|{\mathbf{a}^{\prime}}\right\rangle=\delta(\mathbf{a}-\mathbf{a}^{\prime})$. It allows us to define the single-photon states as

[TABLE]

The infinitesimal propagation for a single-photon state is $\left|{G(z+dz)}\right\rangle=\hat{\mathcal{P}}_{1}(z,dz)\left|{G(z)}\right\rangle$, where $\hat{\mathcal{P}}_{1}(z,dz)\equiv\mathbf{1}+\rmi dz\hat{M}(z)$, with

[TABLE]

being the projection operator for single-photon states, and

[TABLE]

representing the infinitesimal propagation process for single-photon states, which also incorporates a projection operation for single-photon states. Since $\hat{M}(z)$ is Hermitian, it follows that $\hat{\mathcal{P}}_{1}^{{\dagger}}(z,dz)=\mathbf{1}-\rmi dz\hat{M}(z)$.

While $\hat{\mathcal{P}}_{1}$ implies a projection operation for single-photon states, the equivalent for $n$ photons, expressed as $\hat{\mathcal{P}}_{n}\equiv\hat{\mathcal{P}}_{1}^{\otimes n}$, would involve a projection operation for $n$-photon states. For general multi-photon states, one needs to sum over all such operators. Therefore, the unitary operator is given by

[TABLE]

where

[TABLE]

is the identity operator for the complete space, and $\exp_{\otimes}(\cdot)$ implies that all products in its expansion are tensor products and the first term in the expansion contains only the vacuum state $\left|{\rm vac}\right\rangle\left\langle{\rm vac}\right|$. The last line in (9) follows under infinitesimal propagation conditions due to the smallness of $dz$. The sub-leading order term serves as the scintillation operator, $\hat{\mathcal{M}}(z)\equiv\hat{M}(z)\otimes\mathds{1}$. Although it does not involve any projection operations, only one photon is affected by the scintillation process.

For an arbitrary pure state, the infinitesimal propagation gives

[TABLE]

In the limit $dz\rightarrow 0$, it leads to the dynamical equations

[TABLE]

An arbitrary (pure or mixed) state, expressed as a density operator $\hat{\rho}(z)$, would thus obey the dynamical equation

[TABLE]

We can now identify the scintillation operator $\hat{\mathcal{M}}(z)$ with the infinitesimal propagation operator $\hat{P}_{\Delta}$, mentioned above. As such, it describes the unitary evolution of the state propagating through the turbulent atmosphere. However, the detail of the medium is only known in a statistical sense. To have a predictive ability, one needs to consider the ensemble average of the process. For that purpose, it is more convenient to consider the state in terms of its Wigner functional.

2.2 Wigner functional approach

To make the expression in (13) more analytically manageable, we define an auxiliary function for the scintillation operator

[TABLE]

where $\eta$ is an auxiliary variable, and

[TABLE]

The auxiliary variable $\eta$ carries the dimension of a distance and it effectively takes over the role of $dz$. The scintillation operator is recovered by

[TABLE]

The dynamical equation can then be expressed as

[TABLE]

The Wigner functional representation of the dynamical equation in (13) or (17) requires a Wigner functional for the auxiliary function of the scintillation operator $\hat{\mathcal{R}}(z;\eta)$, which represents a linear operation of the form given in (14) and (15). The Wigner functional for such a linear process is given by

[TABLE]

under the assumption that $\mathbf{1}+T=\delta(\mathbf{a}-\mathbf{a}^{\prime})+T(\mathbf{a},\mathbf{a}^{\prime},z)$ is invertable. The derivation of (18) is provided in A.

Note that $(\mathbf{1}+T)|_{\eta=0}=2\delta(\mathbf{a}-\mathbf{a}^{\prime})$ and $(\mathbf{1}-T)|_{\eta=0}=0$. It implies that $(\mathbf{1}+T)^{-1}|_{\eta=0}=\delta(\mathbf{a}-\mathbf{a}^{\prime})/2$. Therefore, if we define $B\equiv 2(\mathbf{1}-T)\diamond(\mathbf{1}+T)^{-1}$, then, by applying derivatives with respect to $\eta$ and setting $\eta=0$, we get $\partial_{\eta}B|_{\eta=0}=-\rmi M(\mathbf{a},\mathbf{a}^{\prime},z)$ and $\partial_{\eta}B^{*}|_{\eta=0}=\rmi M(\mathbf{a}^{\prime},\mathbf{a},z)$, because $M$ is Hermitian.

In terms of Wigner functionals, the dynamical equation in (17) is represented by111The Wigner functional for the product of operators is obtained by expressing the operators via the Weyl transformation in terms of their respective Wigner functionals. The result leads to the functional version of the star-product of the Moyal formalism [35].

[TABLE]

where $W_{\hat{\rho}}(z)[\alpha]$ indicates that the Wigner functional is a function of $z$ and a functional of $\alpha$. Evaluating the derivative with respect to $\eta$ and setting $\eta=0$, we get

[TABLE]

Again, we use an auxiliary function to express the integrand in terms of an exponential function so that we can evaluate the functional integrals

[TABLE]

where $\xi$ is the new auxiliary variable. The argument of the exponential function is linear in $\alpha_{b}$ and $\alpha_{b}^{*}$. Therefore, the functional integration over $\alpha_{b}$ converts the exponential function into a Dirac delta functional, which replaces $\alpha_{a}$ in the argument of the state in terms of $\alpha$ and $M$.

To perform the functional integration over $\alpha_{b}$, we first perform a change of integration variables

[TABLE]

The functional integrations over $\alpha_{b}$ and $\alpha_{a}$ then leads to

[TABLE]

Evaluating the derivative with respect to $\xi$ and setting $\xi=0$, we note that

[TABLE]

Hence, the first-order unitary dynamics for the multi-photon state propagating through a turbulent atmosphere, is obtained

[TABLE]

where we dropped the dependences on $z$ and $\alpha$.

The equation in (25) would not provide the required dynamics when we apply ensemble averaging. The reason is that the refractive index fluctuation has a zero mean, which implies that

[TABLE]

As a result, an ensemble averaging applied to (25), gives

[TABLE]

where $|\mathbf{a}|^{2}\mathbf{1}\equiv|\mathbf{a}|^{2}\delta(\mathbf{a}-\mathbf{a}^{\prime})$. The resulting equation in (27) represents free-space propagation without turbulence.

2.3 Second order

To see the effect of the turbulence after an ensemble averaging, one must go to second order in the scintillation operator. For this purpose, (25) is first integrated over $z$:

[TABLE]

Then we substitute the equation back into itself (making the Fourier domain integrals explicit) to obtain the second-order equation

[TABLE]

In effect, we performed the back-substitution twice, but only kept up to second-order terms, because the third-order terms are at least second order in $dz$. As a result, the $z$-dependences of all Wigner functionals are given as $z_{0}$.

2.4 Ensemble averaging

One can now perform the ensemble averaging process. The zeroth-order terms are unaffected. The first-order terms are given by (26). For the second-order terms, the ensemble averaging involves

[TABLE]

where we substituted in (4) and removed terms that are first order in $N$. One can model the stochastic nature of the refractive index fluctuation by

[TABLE]

where $c$ is the longitudinal spatial frequency component, $\Delta$ is a dimension parameter on the frequency domain and $\chi(\mathbf{a},c)$ is a three-dimensional, normally distributed, complex-valued random function. Since $\tilde{n}$ is a real-valued function, it implies that $\chi^{*}(\mathbf{a},c)=\chi(-\mathbf{a},-c)$. Moreover, it is assumed that this random function is delta-correlated,

[TABLE]

The Kolmogorov power spectral density for the refractive index fluctuations is given by $\Phi_{n}=0.033(2\pi)^{3}C_{n}^{2}|\mathbf{k}|^{-11/3}$, where $C_{n}^{2}$ is the refractive index structure constant and $\mathbf{k}$ is a three-dimensional wave vector — the Fourier domain variable for the refractive index of the medium ($\mathbf{k}=\{2\pi\mathbf{a},2\pi c\}$). The extra power of $(2\pi)^{3}$ comes from a different convention for the definition of the Fourier transform [42].

Next, we evaluate

[TABLE]

where we used (31), (32) and the fact that $\Phi_{n}(\mathbf{a},c)$ is symmetric in all its arguments. Under the Markov approximation, one can replace $\Phi_{n}(\mathbf{a},c_{1})\rightarrow\Phi_{n}(\mathbf{a},0)$, so that

[TABLE]

The expression in (30) then becomes

[TABLE]

where $dz=z-z_{0}$ and where we dropped a term that is second order in $dz$.

2.5 The infinitesimal propagation functional equation

With (26) and (35) substituted into (29), one can apply the limit $dz\rightarrow 0$. The resulting functional differential equation reads

[TABLE]

where

[TABLE]

The equation in (36), which is our main result, represents the infinitesimal propagation functional equation (IPFE) for the evolution of an arbitrary multi-photon state, propagating through atmospheric turbulence under arbitrary conditions. It has the form of a Fokker-Planck equation for the Wigner functional of the state. The first term on the right is a drift term for free-space propagation, the second term is a drift term for the scintillation and the third is the diffusion term for the scintillation.

3 Solving the IPFE

Apart from some trivial solutions, such as the vacuum and thermal states, no general solution has yet been found for the IPFE. Here we make a few observations and suggest approaches that can be used to find solutions for the IPFE.

3.1 Free-space solution

In the case where there is no turbulence, one can use the first-order free-space equation given in (27). A solution for the free-space equation can be obtained by assuming a general form

[TABLE]

where we ignored a possible normalization factor; $A(\mathbf{a},\mathbf{a}^{\prime},z)$, $B(\mathbf{a},\mathbf{a}^{\prime},z)$, and $C(\mathbf{a},\mathbf{a}^{\prime},z)$ are unknown kernel (or covariance) functions; and $\beta(\mathbf{a},z)$ and $\eta(\mathbf{a},z)$ are unknown spectral functions.

Substituting (38) into (27), and separating the resulting equation according to the different combinations of contraction with $\alpha$ and $\alpha^{*}$, we obtain five separate differential equations for the three kernels and the two spectral functions

[TABLE]

Their solutions are

[TABLE]

where $A_{0}(\mathbf{a}_{1},\mathbf{a}_{2})$, $B_{0}(\mathbf{a}_{1},\mathbf{a}_{2})$, and $C_{0}(\mathbf{a}_{1},\mathbf{a}_{2})$ are the initial kernels (or covariance functions) and $\beta_{0}(\mathbf{a})$, and $\eta_{0}^{*}(\mathbf{a})$ are the initial spectral functions. Hence, the free-space solution is

[TABLE]

where

[TABLE]

Another way to view it, is that free-space propagation implies a transformation of the $\alpha$’s, given by

[TABLE]

As such, the result corresponds to simple free-space propagation in a way that one would expect.

3.2 Gaussian solutions

While the solution without turbulence provides an important cross-check to see that the first-order equation is consistent with our expectations, it does not give any indication of the validity of the IPFE. Including the effect of the turbulent medium, one may consider solutions for the IPFE that have the form of a Gaussian Wigner functional. Although, we find only trivial solutions, they do serve as a cross-check for the IPFE.

Consider, for instance the simplest case where the solution has the form

[TABLE]

with $A(\mathbf{a},\mathbf{a}^{\prime},z)$ being an unknown $z$-dependent kernel function. When we substitute this ansatz into (36), evaluate all the functional derivatives and the $z$-derivative, and remove a factor of the original Wigner functional $W_{\hat{\rho}}(z)$ on both sides, we obtain

[TABLE]

We see that the equation does not contain only second-order terms (with single factors of both $\alpha$ and $\alpha^{*}$), but also fourth-order terms (with two factors of $\alpha$ and $\alpha^{*}$ each). As a result, the equation separates into two equations that must be satisfied separately. The second-order equation is

[TABLE]

and the fourth-order equation is

[TABLE]

where we shifted the $\mathbf{a}_{0}$ into the arguments of the kernels so that we can removing all the common factors and the identical contractions with the $\alpha$’s and $\alpha^{*}$’s. Expressing the kernels in terms of sums and differences of their arguments, we have

[TABLE]

It now follows that, if the kernel only depend on the difference in arguments $\mathbf{a}_{d}$, then the sign of $\mathbf{a}_{0}$ in the first term can be changed (since $\Phi_{n}$ is symmetric) and all the terms would cancel. However, this condition would also affect the second-order equation, causing the terms associated with the scintillation to cancel, leaving only the free-space propagation effect. Hence, we end up with

[TABLE]

The only solutions are those where $A(\mathbf{a}-\mathbf{a}^{\prime},z)\propto\delta(\mathbf{a}-\mathbf{a}^{\prime})$. These include the vacuum state and the thermal states.

As a second attempt to solve (36), we introduce some shift terms

[TABLE]

However, after substituting it into (27) and separating the result into separate equations according to the contractions with $\alpha$ and $\alpha^{*}$, one obtains differential equations for the two complex functions given by

[TABLE]

Their solutions are

[TABLE]

However, since $\Lambda$ is a divergent distance, the result is that $\beta(z)=\eta^{*}(z)=0$ for $z>0$. In addition, there are combinations of $\alpha$’s and $\alpha^{*}$’s produced on the right-hand side of the equation that are not represented by any term on the left-hand side. In all these cases, the requirement $\beta(z)=\eta^{*}(z)=0$ is enforced. It thus leads back to the unshifted previous result.

One can proceed to search for a Gaussian solution having an even more complicated expression in the exponent. However, the resulting equation would produce several separate second-order equations and fourth-order equations, which would be difficult to solve, apart from the trivial solution. In the end, the only solutions for the IPFE that can be readily found by assuming a Gaussian form for the Wigner functional are those with isotropic arguments, centered at the origin, such as the vacuum state and the thermal states. Although they do not represent general solutions for the IPFE, they serve as an important cross-check, because it makes sense that these states would not evolve during propagation through a turbulent atmosphere.

3.3 Characteristic functional

Another way to search for a general solution for the IPFE in (36), is to transform the Wigner functional to its corresponding characteristic functional. The hope is that the result would provide insight into the required form for the solution.

The Wigner functional is given in terms of the characteristic functional by

[TABLE]

and the inverse process is given by

[TABLE]

First, we substitute (53) into (36), with $\eta\rightarrow\eta_{0}$, and then we apply the inverse process given in (54). The result has the form

[TABLE]

where

[TABLE]

Note that the positive sign in the second last line comes about due to the different signs of the terms in the exponent in (55).

For the subsequent calculation, we construct $K[\alpha,\eta_{0}]$ with the aid of functional derivatives applied to a source functional of the form

[TABLE]

where $\nu$ and $\mu$ are auxiliary functions. One can define a construction operator that would produce $K[\alpha,\eta_{0}]$ from the source functional by replacing

[TABLE]

in (56). The result is a functional differential operator

[TABLE]

The equation for the characteristic function is then represented by

[TABLE]

The functional integration over $\alpha$ in (57) produces a Dirac delta functional together with an exponential

[TABLE]

When it is substituted into (64), the functional integration over $\eta_{0}$ causes the Dirac delta functional to replace the argument of the characteristic functional with expressions containing the auxiliary functions. The resulting equation has the form

[TABLE]

One can see that the operation of $\mathcal{K}$ on the shifted characteristic functional would lead to a functional differential equation.

We substitute (63) into (66), apply the functional derivatives on the characteristic functional and set $\nu=\mu=0$. The resulting equation reads

[TABLE]

It has the form of a Fokker-Planck equation — the same form as the IPFE for the Wigner functional in (36). As such, the characteristic functional approach does not bring us much nearer to a solution, but it does indicate that the solution of the IPFE and its functional Fourier transform satisfy the same functional differential equation. It thus suggests that some solutions of the IPFE may be eigenstates of the functional Fourier transform. This is indeed true for the Wigner functional of the vacuum state. On the other hand, the Wigner functionals of the thermal states are not equal to their characteristic functionals, but they all have the same form and are therefore (trivial) solutions of these functional differential equations.

3.4 Fixed photon-number solutions

Although the functional Fokker-Planck equation in (36) involves arbitrary multi-photon states, it must also incorporate situations where the state has a fixed number of photons. One can argue that the photon number must remain the same during propagation, because the process from which (36) is derived does not include a loss-mechanism. Therefore, if the states start out with exactly $n$ photons, then they must contain $n$ photons all the time. One can therefore apply the process to a general Fock state, in which the spectrum and its complex conjugate combines to form the density matrix for a state with a fixed number of photons.

For this purpose, we can use a generating functional for the Wigner functionals of fixed-spectrum Fock states. The one-parameter version of this generating functional (88), is derived in B.222The two-parameter version of this generating functional is derived in [30]. Since this generating functional has the form of an isotropic Gaussian, we can use the result obtained in (45), but without removing the factor of the original Wigner functional. One can then evaluate the number of derivatives with respect to the auxiliary variable equal to the number of photons in the state, before setting the auxiliary parameter to zero.

For the single-photon solution, we evaluate one derivative with respect the auxiliary variable, set it to zero and remove the remaining Gaussian factor to obtain

[TABLE]

where $\rho(z)$ is the density matrix in the spatial frequency basis for the single-photon state. The remaining $\alpha$’s and $\alpha^{*}$’s can be removed by functional derivatives. The result would then be appropriately symmetrized. It is equivalent to the single-photon infinitesimal propagation equation (IPE), which has been solved before, using the quadratic structure function approximation [44].

For the bi-photon solution, we evaluate two derivatives with respect the auxiliary variable, set it to zero and remove the remaining Gaussian factor. In this case, one obtains the two-photon equation together with the single-photon equation. After removing the single-photon equation and the $\alpha$’s and $\alpha^{*}$’s with functional derivatives, we obtain the symmetrized equation

[TABLE]

It is the same as the bi-photon IPE, for which a solution has previously been found under the quadratic structure function approximation [45].

These two cases show that the multi-photon Fokker-Planck equation in (36) is consistent with previous results. An equation for any fixed photon-number state can be obtained in this way. However, they become progressively more difficult to solve.

3.5 Polynomial solution

Due to the requirement for normalizability of the state, it is not expected that the solution of the IPFE would be a polynomial of finite order. Such polynomials are unbounded and therefore not normalizable. However, any normalizable functional can be represented as a polynomial of infinite order — a transcendental polynomial.

Here, we express the solution in the form of such a transcendental polynomial and substitute it into the IPFE. The result is an infinite set of uncoupled differential equations. If all these differential equations can somehow be solved, one would have a general solution for the IPFE.

The expression of the transcendental functional polynomial is

[TABLE]

where $\odot$ represents multiple contractions, connecting the kernels $H_{m,n}$ to all the $\alpha^{*}$’s and $\alpha$’s. Note that the kernels are symmetric with respect to any perturbation among $\{\mathbf{a}_{1}...\mathbf{a}_{m}\}$ or among $\{\mathbf{a}_{1}^{\prime}...\mathbf{a}_{n}^{\prime}\}$.

When we substitute the expression for the functional polynomial into the IPFE in (36), we get an expression that can be separated into separate equations for every order. These equations are of the form

[TABLE]

Due to the integer factors, which are produced by the functional derivatives, the equations simplify for lower values of $m$ and $n$. For $m\geq 2$ and $n\geq 2$ the equations are of the general form, given in (71). We note that the solution for $H_{n,m}(z)$ is always the complex conjugate of the solution for $H_{m,n}(z)$.

The first few equations in the sequence are

[TABLE]

The first equation in (72) implies that $H_{0,0}$ is a constant independent of $z$. It can thus act as a normalization constant for the state. The solutions for $H_{1,0}(\mathbf{a}_{1},z)$ and its complex conjugate are

[TABLE]

where $B_{1,0}(\mathbf{a}_{1})$ and $B_{0,1}(\mathbf{a}_{1}^{\prime})$ are the initial spectral functions at $z=0$. These solutions decay to zero on a scale equivalent to $z\sim 1/(k^{2}\Lambda)$. The equation for $H_{1,1}(\mathbf{a}_{1},\mathbf{a}_{1}^{\prime},z)$ is equivalent to the IPE for a single-photon state, which has been solved in [44] under the quadratic structure function approximation. The solution is in the form of an integral

[TABLE]

where $L(\mathbf{a}_{1},\mathbf{a}_{1}^{\prime},\mathbf{u},z)$ is a kernel function given in [44] for the quadratic structure function approximation.

All the remaining equations in the sequence of equations given in (71) still need to be solved. It is not at this time certain that a generic solution can be found for arbitrary orders in (71). If such generic solutions can be found, it is still not certain that such solutions would produce results that can be summed into a closed form to serve as a general solution for the functional differential equation.

4 Conclusion

The evolution of an arbitrary photonic quantum state propagating through a turbulent atmosphere is considered. An infinitesimal propagation approach is used to obtain a Fokker-Planck equation for the Wigner functional of the state — the IPFE. Apart from the $z$-derivatives, all the other derivatives in the equation are functional derivatives.

We do not provide a general solution for the IPFE. However, we show that without turbulence, the solution is in the form as expected for free-space propagation without turbulence. Considering the possibility for a solution in the form of a Gaussian functional, we find that the only solutions of this form are the thermal states, which include the vacuum state.

The use of the characteristic functional to obtain a general solution is also considered. It leads to another functional differential equation for the characteristic functional that has the same form as the IPFE. One can therefore concluded that a general solution of the IPFE would to a set of functional that are mapped back onto itself by the functional Fourier transform.

We show that the functional Fokker-Planck equation is consistent with previous solutions that were obtained for single-photon and bi-photon states. It follows from the observation that the Fokker-Planck equation maintains photon-number, because it does not involve a loss mechanism.

Finally, we propose a general approach to find a solution, expressing it as an infinite order polynomial functional, which can then be decoupled into an infinite set of uncoupled differential equations. It is not sure whether an approach exists with which all such equations can be solved.

Acknowledgement

The research for this work was supported in part by the National Research Foundation (NRF) of South Africa (Grant Number: 118532).

Appendix A Wigner functional of a linear process

Consider a linear operation for a single-photon state, expressed as a generalization of the single-photon projection operator, given by

[TABLE]

where $P(\mathbf{k},\mathbf{k}^{\prime})$ represents the linear transformation to be performed on the angular spectrum (Fourier domain wave function) of the photon and $\left|{\mathbf{k}}\right\rangle$ denotes the elements of a momentum basis. The generalization to arbitrary numbers of particles can then be represented as a sum of tensor products

[TABLE]

where $\exp_{\otimes}(\cdot)$ is defined such that all products in its expansion are tensor products and the first term in the expansion involves only the vacuum state $\left|{\rm vac}\right\rangle\left\langle{\rm vac}\right|$.

To represent the operator as a Wigner functional, we use a coherent state assisted approach [30]. It requires that the spatial degrees of freedom be represented in the Fourier domain, as in (76). Then, we need to compute the overlap of the operator with coherent states of both sides

[TABLE]

where $||\alpha||^{2}\equiv\alpha^{*}\diamond\alpha$.

We now use the result in (77) to compute the Wigner functional for a general linear process. Substituting it into the expression for the Wigner functional (see [30]), and performing the functional integrations over $\alpha_{1}$ and $\alpha_{2}$, we obtain

[TABLE]

where it is assumed that $(\mathbf{1}+P)$ is invertable and

[TABLE]

Appendix B Generating function for fixed-spectrum Fock states

Here, we use the coherent state assisted approach [30] to derive a generating function for the Wigner functionals of the fixed-spectrum Fock states. The latter are defined as

[TABLE]

where $F$ is a normalized complex-valued parameter function, and $\hat{a}^{{\dagger}}$ is the creation operator. The overlap between the density operator for such a Fock state and two fixed-spectrum coherent states is

[TABLE]

where $\langle\alpha_{1},F\rangle$ represents an inner product between the parameter functions of the fixed-spectrum coherent state and the fixed-spectrum Fock state, respectively. We represent the overlap in terms of a generating function

[TABLE]

where $\eta$ is an auxiliary parameter. The overlap is recovered by

[TABLE]

Substituting $\left\langle{\alpha_{1}}\right|\hat{A}\left|{\alpha_{2}}\right\rangle\rightarrow{\cal K}$ into the expression for the Wigner functional (see [30]), we obtain a functional integral expression for the generating function of the Wigner functionals of fixed-spectrum Fock states, given by

[TABLE]

After performing the functional integrations over $\alpha_{1}$ and $\alpha_{2}$, we get

[TABLE]

Thanks to the fact that the parameter functions of the fixed-spectrum Fock states are normalized ($F\diamond F^{*}=1$), one can show that

[TABLE]

and

[TABLE]

It then follows that

[TABLE]

where ${\cal N}_{0}$ is the normalization constant. The Wigner functionals for the individual Fock states are given by

[TABLE]

where $L_{n}(\cdot)$ represents the Laguerre polynomial of order $n$.