Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group
Journal Home About Issues in Progress Current Issue All Issues Feature Issues

Computation of optimal beams in weak turbulence

Qin Li, Anjali Nair, and Samuel N Stechmann

Open Access Open Access

Abstract

When an optical beam propagates through a turbulent medium such as the atmosphere or ocean, the beam will become distorted. It is then natural to seek the best or optimal beam that is distorted least, under some metric such as intensity or scintillation. We seek to maximize the light intensity at the receiver using the paraxial wave equation with weak-fluctuation as the model. In contrast to classical results that typically confine original laser beams to be from a special class, we allow the beam to be general, which leads to an eigenvalue problem of a large-sized matrix with each entry being a multi-dimensional integral. This is an expensive and sometimes infeasible computational task in many practically reasonable settings. To overcome this expense, in a change from past calculations of optimal beams, we transform the calculation from physical space to Fourier space. Since the structure of the turbulence is commonly described in Fourier space, the computational cost is significantly reduced. This also allows us to incorporate some optional turbulence assumptions, such as homogeneous-statistics assumption, small-length-scale cutoff assumption, and Markov assumption, to further reduce the dimension of the numerical integral. The proposed methods provide a computational strategy that is numerically feasible, and results are demonstrated in several numerical examples. These results provide further evidence that special beams can be defined to have beam divergence that is small.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

When optical beams propagate through a random medium, they are subject to distortions that lead to unwanted phenomena like intensity reduction and scintillation [1,2]. It is then often desirable to search for beams that are optimal under certain criteria.

The criterion taken in this paper is to maximize light intensity at the receiver. Mathematically it has been proved to be a very clean problem. It was shown that the optimal beam under this criterion has to be coherent [3,4]. Moreover, it was proved that this optimal, coherent beam is associated with the largest eigenvalue of a matrix or operator:

$$\mathcal{H}=\mathbb{E}[H], \quad\text{with}\quad H(X_1,X_2)=\int\limits_{X'\in\mathcal{R}}h(X_1,X')h^{{\ast}}(X_2,X')\mathrm{d}X'\,.$$
Here $h(X_1,X')$ is the propagator (Green’s function) that propagates light from the origin $X_1$ to $X'$. $\mathcal {R}$ is the receiver region, and the expected value $\mathbb {E}$ takes average over all possible configurations of the random medium. By definition, $H(X_1,X_2)$ is a Hermitian kernel with non-negative eigenvalues. Since the task of finding an optimal beam is equivalent to the task of finding the eigenvector or eigenfunction associated with the largest eigenvalue of $\mathcal {H}=\mathbb {E}[H]$, what remains is to mathematically and computationally formulate $\mathcal {H}$.

This turns out to be a very challenging task in practice. There are two obstacles. First, the expected value $\mathbb {E}$ means all random media configurations need to be taken into account. Numerically, suppose we utilize Monte Carlo to sample many representative configurations; then a different $H$ needs to be evaluated for each of these samples. Second, computing each $H$ is already a major challenge. It amounts to evaluating all Green’s functions $h$ that maps every $X\in \mathcal {A}$, the aperture, to every $X'\in \mathcal {R}$. This is a drastic cost for a problem in three-dimensional space, $\mathbb {R}^{3}$.

These difficulties have made it practically impossible to find the optimal beams. To overcome the difficulties, we consider the following context and approach. As a first aspect, we consider perturbation theory in the weak fluctuation regime. One convenience of studying the weak fluctuation regime is that the nonlinear dependence on the randomness is now linearized around the deterministic setup. Instead of sampling many random configurations of the media, and computing all the associated Green’s functions before taking the expected value $\mathbb {E}$ over the media randomness, now one can find an analytical expression, and the expected value $\mathbb {E}$ gets directly applied to the medium. Eventually it is the moments of the medium’s randomness that enter into the analytical expression, and, if desired, empirical models such as the Kolmogorov spectrum can be employed in a straightforward manner. While analytical expressions for the moments of the solution to the Paraxial Wave Equation (PWE) have been computed before using perturbation theory, they have been for specific classes of sources like Gaussian beams [2]. Here we compute the second moment for general source functions for general weakly turbulent models, and use this to compute optimal beams.

However, such a formulation of the analytical expression for $\mathcal {H}$ is still hard to compute. The fine discretization and the high dimensionality issue encountered in the nonlinear regime are still present here. As will be presented in Section 4.1, this expression results in an $8$-fold integral in a full 3D simulation for each entry of $\mathcal {H}$, and is numerically prohibitive. As some intuition toward reducing the computational cost, one observation is that many empirical models for atmospheric turbulence have special structure on the Fourier domain, and thus it is expected that calculations can be significantly simplified if conducted on the Fourier side. The Fourier-space perspective, however, has not been taken previously in the literature on optimal beams. In our calculation, we fully take advantage of the encoded turbulence structure in the Fourier domain, and convert our calculation to that of $\mathcal {H}_\mathrm {f}$, the Fourier-space counterpart of $\mathcal {H}$; see Section 4.2. This allows us to incorporate some well-known assumptions on the medium’s structure with ease. We use homogeneous assumption, small-length-scale cutoff assumption, and Markov assumption sequentially, and observe the reduction from an $8$-fold integral, to a $6$-fold, to a $2$-fold and finally to a $1$-fold integral respectively, making it gradually more and more numerically feasible, see Section 6. Some properties of $\mathcal {H}$ and $\mathcal {H}_\mathrm {f}$ are shown in Section 5, and numerical evidences are provided in Section 7.

The methods here provide an additional contribution to the literature on designing optimal beams in the presence of turbulence. As mentioned above, a mathematical formulation for optimal beams was given in [3,4]. Few numerical examples of optimal beams have been presented. In the present paper, we find optimal beams for the paraxial wave equation in the weak turbulence setting where the refractive index is allowed to be a general random function as a representation of turbulent fluctuations. Past work has investigated optimal beams under other frameworks that introduce additional assumptions, such as the phase screen model and the extended Huygens–Fresnel (eHF) principle [5–9], or searching within a special class of beams such as the Gaussian Schell-model beams [10].

More broadly, the problem of finding optimal beams is motivated by much past work on increasing beam intensity or reducing scintillation. From past work on these topics, it is generally known that coherent beams maximize the received intensity, although coherent beams do not perform well with respect to other metrics like scintillation; indeed, turbulence has a degrading effect on many types of beams [4,11–15]. To reduce the effects of scintillation, past work has identified the importance of partially coherent beams [12,13,16–21]. The problem of minimizing scintillation is more challenging from a computational perspective, and we therefore focus here on the more tractable problem of maximizing intensity.

2. Problem setup and notation

Consider a laser beam source with optical field $\phi$ that is located in a transmitter region $\mathcal {A}$, and let $U$ denote the field received at the receiver region $\mathcal {R}$. For waves that are governed by a linear equation, $U$ can be represented as:

$$U(X')=\int\limits_{X\in\mathcal{A}}h_\omega(X,X')\phi(X)\mathrm{d}X\,.$$
This $h_\omega (X,X')$ is termed the propagator that sends light information from the aperture location $X$ to the receiver location $X'$, so it naturally includes the medium information. This propagator depends on the specific configuration of the media, as indicated by the subindex $\omega$. We surpress this subindex for the conciseness of the notation when the context is clear.

There are two sources of randomness. The source $\phi$ is generated by the laser and contains random fluctuation, and the medium also presents turbulence [22]. When a specific source and medium configuration is fixed, and the light intensity can be calculated as

$$I_i(X')=|U(X')|^{2}=\int\nolimits_{X_{1},X_2\in\mathcal{A}}h(X_1,X')\phi(X_1)h^{{\ast}}(X_2,X')\phi^{{\ast}}(X_2)\mathrm{d}X_1\mathrm{d}X_2,\quad X'\in\mathcal{R}\, ,$$
where we use superscript $^{*}$ to denote the complex conjugate. Taking the average with respect to both the source fluctuations and media turbulence, we have the averaged intensity:
$$I(X')=\mathbb{E}\Big[\int\limits_{X_{1}\in\mathcal{A}}\int\limits_{X_2\in\mathcal{A}}h(X_1,X')J(X_1,X_2)h^{{\ast}}(X_2,X')\mathrm{d}X_1\mathrm{d}X_2\Big],\quad X'\in\mathcal{R},$$
where $J(X_1,X_2)=\langle \phi (X_1)\phi ^{\ast }(X_2)\rangle$ is the mutual intensity function of the source. Throughout the paper we use $\langle \cdot \rangle$ to denote the averaging with respect to the source randomness, and $\mathbb {E}[\cdot ]$ to denote the averaging with respect to the randomness in the medium.

Designing optimal laser beam amounts to tuning $J$ that achieves the highest intensity. To do so, the intensity maximization problem is formulated as

$$\max_{J} I_\mathcal{R}=\max_{J}\int\limits_{X'\in\mathcal{R}}I(X')\mathrm{d}X', \quad\text{subject to}\quad I_0 = 1,$$
where $I_0=\int\nolimits _{X\in \mathcal {A}}J(X,X)\mathrm {d}X=\int\nolimits _{X\in \mathcal {A}}\langle \phi (X)\phi ^{\ast }(X)\rangle \mathrm {d}X$ is the light intensity at the initial state.

A useful reformulation of (5) follows from noting that, by definition, $J$ is an Hermitian matrix. Hence, following earlier work [23], we can write $J$ in terms of its coherent mode expansion as

$$J(X_1,X_2)=\sum_k\alpha_k\psi_k(X_1)\psi_k^{{\ast}}(X_2), \quad X_1,X_2\in\mathcal{A},$$
where $\{\alpha _k\}$ are non-negative weights and $\{\psi _k\}$ is a set of functions that are orthonormal over the transmitter region. Considering all candidates for $J$ is then equivalent to considering all $\{\alpha _k,\psi _k\}$ pairs [24–26]. This formulation allows us to rewrite the constraint as
$$I_0=\sum_k\alpha_k=1,$$
and the objective function $I_\mathcal {R}$ becomes:
$$I_\mathcal{R}=\sum_k\alpha_k\int\limits_{X_1\in\mathcal{A}}\int\limits_{X_2\in\mathcal{A}}\psi_k(X_1)\mathbb{E}[H(X_1,X_2)]\psi_k^{{\ast}}(X_2)\mathrm{d}X_1\mathrm{d}X_2,$$
where for all $X_1,X_2\in \mathcal {A}$ we use the same notation as in (1) and set $\mathcal {H} = \mathbb {E}[H]$ with $H(X_1,X_2)=\int\limits _{X'\in \mathcal {R}}h(X_1,X')h^{\ast }(X_2,X')\mathrm {d}X'$. By definition, $H(X_1,X_2)$ is a Hermitian kernel with non-negative eigenvalues.

With this reformulation, one has an analytical solution for maximizing (8) when constrained on (7). As summarized in [4], when $\mathbb {E}[H]$ is given, the optimal solution to (5) has the form of a coherent mode, namely:

$$J(X_1, X_2) = \alpha_1\psi_1(X_1)\psi_1^{{\ast}}(X_2), \quad\text{with}\quad \alpha_1 = 1,$$
where $\psi _1$ is the eigenfunction of $\mathbb {E}[H]$ associated with the largest eigenvalue. Since finding eigenvalues and eigenfunctions of a given matrix is mathematically straightforward, the only obstacle of identifying the optimal beam configuration lies in formulating $\mathbb {E}[H]$.

3. Calculation of transfer function

While many of the main ideas here are applicable to general types of waves, in what follows we consider here the setting of optics with the paraxial wave equation (PWE), for concreteness. The PWE is also known as the parabolic approximation [27–29].

As a starting point, we recall the form of the PWE:

$$\nabla_\perp^{2} A +2ik\frac{\partial A}{\partial z} + k^{2}\left(\frac{n^{2}}{n_r^{2}}-1\right)A =0,$$
where $A$ is the complex signal amplitude, $k$ is the reference wave number of the source, and $n=n(x,y,z)$ is the index of refraction, which is a function of turbulent anomalies of air temperature. Here $X=(x,y)$ denotes the perpendicular direction, with $\nabla _\perp ^{2}=\partial _x^{2}+\partial _y^{2}$, and $z$ denotes the transverse direction of propagation. The source is placed at $z=0$ and the receiver is located at $z=Z$. As a consequence:
$$A(X,z=0)=\phi(X), \quad\text{and}\quad U(X') = A(X',z=Z)\,.$$

In the weak fluctuation regime, the refraction index of the medium centers around a constant $n_r$:

$$n^{2}=n_r^{2} + \epsilon n_1^{2}$$
where $\epsilon$ encodes the small amplitude of the fluctuation. (To arrive at (10), start from the ansatz $n = n_r +\tilde {n}$, so that $n^{2}$ expands as $n^{2} = n_r^{2} + 2\tilde {n}n_r + \tilde {n}^{2}$. Under weak fluctuations, $\tilde {n}$ is substantially smaller than $n_r$. So we exclude $\tilde {n}^{2}$ from the calculations, and change notation from $2\tilde {n}n_r$ to $\epsilon n_1^{2}$.) In this setting, we can use perturbation theory [2,30], and follow the classical asymptotic expansion technique [31,32] to expand $A$ to be:
$$A=A_0+\epsilon A_1+\epsilon^{2} A_2+\cdots\,.$$
Inserting this ansatz into (9) and balancing each order of $\epsilon$, we have:
$$\epsilon^{0}: \qquad\nabla_\perp^{2} A_0 +2ik\frac{\partial A_0}{\partial z} =0,$$
$$\epsilon^{1}: \qquad\nabla_\perp^{2} A_1 +2ik\frac{\partial A_1}{\partial z} ={-}k^{2}\frac{n_1^{2}}{n_r^{2}}A_0,$$
$$\epsilon^{2}: \qquad\nabla_\perp^{2} A_2 +2ik\frac{\partial A_2}{\partial z} ={-}k^{2}\frac{n_1^{2}}{n_r^{2}}A_1\,.$$
Denoting $V(x,y,z)=\frac {n_1^{2}}{n_r^{2}}$, we can express the solutions to $A_i$ explicitly, as summarized in the following subsection.

3.1 Hierarchical solution to PWE in the weak fluctuation regime

Suppose $X\in \mathbb {R}^{d}$. Recall the equation for $A_0$ in (12) with $\phi$ as the initial source term:

$$\nabla_\perp^{2} A_0 +2ik\frac{\partial A_0}{\partial z}=0, \quad A(X,0)=\phi(X)\,.$$
For $K\in \mathbb {R}^{d}$, the Fourier transform is $\hat {f}(K,z)=\int\limits _{X\in \mathbb {R}^{d}}f(X,z)e^{-iK\cdot X}\mathrm {d}X$ and thus the inverse is $f(X,z)=\frac {1}{(2\pi )^{d}}\int\limits _{K\in \mathbb {R}^{d}}\hat {f}(K,z)e^{iK\cdot X}\mathrm {d}K$. We perform the Fourier transform of Eq. (15) to have: $-\|K\|^{2}\hat {A}_0+2ik\frac {\partial \hat {A}_0}{\partial z}=0$. This ODE (ordinary differential equation) has an explicit solution:
$$\hat{A}_0(K,z)=\hat{\phi}(K)\hat{G}(K,z)\,\quad \text{with}\quad \hat{G}(K,z)=\exp(-\frac{iz}{2k}\|K\|^{2})\,.$$
$\hat {G}$ is the Fourier transform of the uniform-medium Green’s function given by
$$G(X,z)=\Big(\frac{1}{2\pi}\Big)^{d/2}\Big(\frac{k}{iz}\Big)^{d/2}\exp\Big(\frac{ik}{2z}\|X\|^{2}\Big)\,.$$
Taking the inverse Fourier transform we obtain:
$$A_0(X,z)=\phi*G(X,z),$$
where $\ast$ is the convolution sign. Similarly, by applying the Fourier transform to the equation for $A_1$ in (13), and define $S(K,z)=\widehat {VA_0}(K,z)$, we have: $-\|K\|^{2}\hat {A}_1+2ik\frac {\partial \hat {A}_1}{\partial z}=-k^{2}S$. Its solution is
$$\hat{A}_1(K,z)=\frac{ik}{2}\int_{z'=0}^{z}S(K,z')\hat{G}(K,z-z')\mathrm{d}z',$$
and consequently:
$$A_1(X,z)=\frac{ik}{2}\int_{z'=0}^{z}\left[(VA_0)({\cdot} ,z')*G({\cdot} ,z-z')\right](X)\mathrm{d}z',$$
where $G(\cdot,z-z')$ denotes the Green’s function evaluated at the shifted $z$ coordinate $z-z'$, and
$$\left[VA_0({\cdot},z')\ast G({\cdot},z-z')\right](X) = \int_{X'\in\mathbb{R}^{d}}VA_0(X',z')G(X-X',z-z'))\mathrm{d}X'\,.$$
The solution to $A_2$ is similarly obtained:
$$A_2(X,z)=\frac{ik}{2}\int_{z'=0}^{z}dz'\left[(VA_1)({\cdot} ,z')*G({\cdot} ,z-z')\right](X)\,.$$

3.2 Calculation of the propagator

From the definition in (2) and the expansion in (11), we express

$$h = h_0 + \epsilon h_1+\epsilon^{2} h_2+\text{high orders}$$
where $h_i(X,X')$ is the propagator that collects the contribution of $\phi (X)$ in $A_i(X',Z)$. According to (18), we have
$$U(X')=A_0(X',Z)=\int\limits_{X\in\mathcal{A}}\phi(X)G(X'-X,Z)\mathrm{d}X,$$
which naturally makes
$$h_0(X,X')=G(X'-X,Z)\,.$$
For the next order, the rearrangement of (20) gives
$$h_1(X,X')=\frac{ik}{2}\int_{z_1=0}^{Z}\mathrm{d}z_1\int\limits_{X_1'\in \mathbb{R}^{d}}V(X_1',z_1)G(X'-X_1',Z-z_1)G(X_1'-X,z_1)\mathrm{d}X_1'\,.$$
One way to interpret this formula is to set a screen at $z=z_1$ and collect all plane wave contribution from $(z=0,X)$ to $(z_1,X'_1)$, and view them as wave fronts by continuing propagating the plane waves to $(z=Z,X')$, coinciding with the Huygens–Fresnel principle. This interpretation suggests that $h_1$ collects the information of waves that gets scattered once.

Similarly, we have the expression for $h_2$ as

$$\begin{aligned}h_2(X,X')&={-}\frac{k^{2}}{4}\int\limits_{z_1=0}^{Z}\int\limits_{X_1'\in\mathbb{R}^{d}}V(X_1',z_1)G(X'-X_1',Z-z_1)\mathrm{d}X_1'\mathrm{d}z_1\\ &\int\limits_{z_2=0}^{z_1}\int\limits_{X_2'\in\mathbb{R}^{d}}V(X_2',z_2)G(X_1'-X_2',z_1-z_2)G(X_2'-X,z_2)\mathrm{d}X_2'\mathrm{d}z_2 \end{aligned}.$$
Similar to the interpretation above, one can view it as a collection of waves that scatter twice, once at $(z_1,X_1')$ and once at $(z_2,X_2')$. The media information is included only through the scattering points.

3.3 Calculation of the propagator on the Fourier space

The calculation can be repeated on the Fourier domain as well. Since the presentation of Eqs. (12)–14 on Fourier domain are all linear, there exists a propagator, termed $g$, so that the solution:

$$\hat{U}(K)=\hat{A}(K,Z)=\int\limits_{K'\in\mathbb{R}^{d}}\hat{\phi}(K')g(K',K)\mathrm{d}K'\,.$$
As done in the physical domain, we expand $g$ in terms of powers of $\epsilon$ to have
$$g = g_0 + \epsilon g_1 + \epsilon^{2} g_2 + \text{high orders},$$
with $g_i$ presenting the contribution from $\hat {A}_i$. Comparing to (16) it is straightforward to have:
$$g_0(K',K)=\hat{G}(K,Z)\delta(K-K')\,.$$

To compute $g_{1}$, we recall (19) and plug in (16):

$$\begin{aligned}\hat{U}_1(K)&=\frac{ik}{2}\int\limits_{z_1=0}^{Z}S(K,z_1)\hat{G}(K,Z-z_1)\mathrm{d}z_1\\ &=\frac{ik}{2(2\pi)^{d}}\int\limits_{z_1=0}^{Z}\hat{G}(K,Z-z_1)\mathrm{d}z_1\int\limits_{K'\in\mathbb{R}^{d}}\hat{V}(K-K',z_1)\hat{G}(K',z_1)\hat{\phi}(K')\mathrm{d}K', \end{aligned}$$
which suggests:
$$g_1(K',K)=\frac{ik}{2(2\pi)^{d}}\int\limits_{z_1=0}^{Z}\hat{G}(K,Z-z_1)\hat{G}(K',z_1)\hat{V}(K-K',z_1)\mathrm{d}z_1\,.$$

Similarly, we have

$$\begin{aligned}g_2(K',K)=&-\Big(\frac{1}{2\pi}\Big)^{2d}\frac{k^{2}}{4}\int\limits_{z_1=0}^{Z}\hat{G}(K,Z-z_1)\mathrm{d}z_1\int\limits_{K{\prime\prime}\in\mathbb{R}^{d}}\hat{V}(K-K{\prime\prime},z_1)\mathrm{d}K{\prime\prime}\\ &\times\int\limits_{z_2=0}^{z_1}\hat{G}(K{\prime\prime},z_1-z_2)\hat{G}(K',z_2)\hat{V}(K{\prime\prime}-K',z_2)\mathrm{d}z_2\,. \end{aligned}$$
The formula (28), (30) and (31) are suggesting that the scattering takes place at location $z_1$ for once-scattered wave, and $z_1$, $z_2$ for twice-scattered wave respectively, and only through scattering, the wave picks up the frequency information from the media.

4. Calculation of $\mathcal {H}$

The explicit formulation of the propagators, $h$ on the physical space and $g$ on the Fourier domain, allows us to compute $\mathcal {H}$. We derive the formula in this section.

4.1 Representation in physical space

Noting the definition (1) for the expansion in (22), we rewrite:

$$\mathcal{H}=\mathcal{H}_{00}+\epsilon^{2}\big(\mathcal{H}_{02}+\mathcal{H}_{11}+\mathcal{H}_{20}\big),$$
where $\mathcal {H}_{ij}(X_1,X_2)=\mathbb {E}\Big [\int\limits _{X'\in \mathcal {R}}h_i(X_1,X')h^{\ast }_j(X_2,X')\mathrm {d}X'\Big ]$. We should note that $\mathcal {H}_{01}=0$ and $\mathcal {H}_{10}=0$ and thus are dropped out of the expansion. We now calculate each term. First, recalling (24), we have:
$$\begin{aligned}\mathcal{H}_{00}(X_1,X_2)&=\int\limits_{X'\in\mathcal{R}}G(X_1-X',Z)G^{{\ast}}(X_2-X',Z)\mathrm{d}X'\\ &=\Big(\frac{1}{2\pi}\Big)^{d}\Big(\frac{k}{Z}\Big)^{d}\int\limits_{X'\in\mathcal{R}}\exp\Big(\frac{ik}{2Z}(\|X_1-X'\|^{2}-\|X_2-X'\|^{2})\Big)\mathrm{d}X'\\ &=G(X_1,Z)G^{{\ast}}(X_2,Z)\int\limits_{X'\in\mathcal{R}}\exp\Big(-\frac{ik}{Z}X'\cdot(X_1-X_2)\Big)\mathrm{d}X'\,. \end{aligned}$$
Similarly, defining
$$\mathbb{E}[V(X_1,z_1)V^{{\ast}}(X_2,z_2)]=\Gamma_V(X_1,z_1,X_2,z_2),$$
and plugging (25) in $\mathcal {H}_{11}$, we have:
$$\begin{aligned}\mathcal{H}_{11}(X_1,X_2) &=\frac{k^{2}}{4}\int\limits_{z_1,z_2=0}^{Z}\mathrm{d}z_1\mathrm{d}z_2\int\limits_{\mathbb{R}^{d}\times \mathbb{R}^{d}}\Gamma_V(X_1',z_1,X_2',z_2)G(X_1-X_1',z_1)G^{{\ast}}(X_2-X_2',z_2)\mathrm{d}X_1'\mathrm{d}X_2'\\ &\times\int\limits_{X'\in\mathcal{R}}G(X'-X_1',Z-z_1)G^{{\ast}}(X'-X_2',Z-z_2)\mathrm{d}X'. \end{aligned}$$
The formulations of $\mathcal {H}_{02}$ and $\mathcal {H}_{20}$ naturally follow as
$$\begin{aligned}\mathcal{H}_{02}(X_1,X_2)&={-}\frac{k^{2}}{4}\int\limits_{z_1=0}^{Z}\int\limits_{z_2=0}^{z_1}\mathrm{d}z_2\mathrm{d}z_1\int\limits_{X_1'\in\mathbb{R}^{d}}\int\limits_{X_2'\in\mathbb{R}^{d}}\Gamma_V(X_1,'z_1,X_2',z_2)G^{{\ast}}(X_1'-X_2',z_1-z_2)G^{{\ast}}(X_2'-X_2,z_2)\mathrm{d}X_1'\mathrm{d}X_2'\\ &\times\int\limits_{X'\in\mathcal{R}}G(X'-X_1,Z)G^{{\ast}}(X'-X_1',Z-z_1)\mathrm{d}X', \end{aligned}$$
and $\mathcal {H}_{20}=\mathcal {H}^{\ast }_{02}$. We note that this is a very complicated formulation and can hardly be of practical use in reality. More specifically, for each fixed $(X_1,X_2)$, the computation of $\mathcal {H}_{02}$, for example, amounts to $8$ dimensional integral when $d=2$.

4.2 Representation in the Fourier domain

The computational cost is prohibitive if the calculation is conducted on the physical domain. However, in reality, assumptions on the atmospheric turbulence are typically represented on the Fourier domain of $\Gamma _V$. Naturally, if one can repeat the process on the Fourier space and incorporate the assumptions, computational difficulty can potentially be reduced. We explore such possibility in this section.

Recall that the intensity to be maximized is:

$$I_\mathcal{R}=\int\limits_{X'\in\mathcal{R}}\mathbb{E}\left[\langle U(X')U^{{\ast}}(X')\rangle\right]\mathrm{d}X'=\int\limits_{X'\in\mathbb{R}^{d}}\mathbb{E}\left[\langle U(X')w(X')U^{{\ast}}(X')w^{{\ast}}(X')\rangle\right]\mathrm{d}X',$$
where $w$ is the indicator function that takes value $1$ within the window $x\in \mathcal {R}$, and $0$ outside. Using Parseval’s equality, this translates to:
$$I_\mathcal{R}=\frac{1}{(2\pi)^{d}}\int\limits_{K\in\mathbb{R}^{d}}\mathbb{E}\left[\langle W(K)W^{{\ast}}(K)\rangle\right]\mathrm{d}K,$$
where $W=\widehat {Uw}$ denotes the Fourier transform:
$$W(K)=\widehat{Uw}(K) = \frac{1}{(2\pi)^{d}}\hat{U}\ast\hat{w}\, .$$
Expanding this convolution in (37):
$$I_\mathcal{R} =\frac{1}{(2\pi)^{3d}}\int\limits_{K_1,K_2\in\mathbb{R}^{d}}\mathbb{E}\big[\langle\hat{U}(K_1)\hat{U}^{{\ast}}(K_2)\rangle\big]\left( \;\int\limits_{K\in\mathbb{R}^{d}}\hat{w}(K-K_1)\hat{w}^{{\ast}}(K-K_2)\mathrm{d}K \right) \mathrm{d}K_1\mathrm{d}K_2\,,$$
we note that the window function term can be simplified:
$$\begin{aligned}\frac{1}{(2\pi)^{d}} \int\limits_{K\in\mathbb{R}^{d}}\hat{w}(K-K_1)\hat{w}^{{\ast}}(K-K_2)\mathrm{d}K&=\int\limits_{X'\in\mathcal{R}}\int\limits_{X{\prime\prime}\in\mathcal{R}}\mathrm{d}X'\mathrm{d}X{\prime\prime}e^{i(K_1\cdot X'-K_2\cdot X{\prime\prime})}\delta(X'-X{\prime\prime})\\ &=\int\limits_{X'\in\mathcal{R}}\mathrm{d}X'e^{i(K_1-K_2)\cdot X'}=\hat{w}(K_1-K_2)\, \end{aligned}$$
and that $\hat {U}(K)=\sum _j\epsilon ^{j}\hat {U}_j(K)$, we adopt the representation (27) to further rewrite $I_\mathcal {R}$ to:
$$\begin{aligned}I_\mathcal{R} &=\sum_{ij}\frac{\epsilon^{i+j}}{(2\pi)^{2d}}\int\limits_{K_1,K_2\in\mathbb{R}^{d}}\mathbb{E}\big[\langle\hat{U}_i(K_1)\hat{U}_j^{{\ast}}(K_2)\rangle\big]\hat{w}(K_1-K_2)\mathrm{d}K_1\mathrm{d}K_2\\ &=\sum_{ij}\frac{\epsilon^{i+j}}{(2\pi)^{2d}}\int\limits_{K_1,K_2\in\mathbb{R}^{d}}\langle\hat{\phi}(K_1)\hat{\phi}^{{\ast}}(K_2)\rangle\mathrm{d}K_1\mathrm{d}K_2\\ &\times\int\limits_{K_1',K_2'\in\mathbb{R}^{d}}\mathbb{E}\big[g_i(K_1,K_1')g_j^{{\ast}}(K_2,K_2')\big]\hat{w}(K_1'-K_2')\mathrm{d}K_1'\mathrm{d}K_2'\,. \end{aligned}$$
Recall that $J(X_1,X_2)=\langle \phi (X_1)\phi ^{*}(X_2)\rangle$ and defining
$$\hat{J}(K_1,K_2)=\int\limits_{X_1\in\mathbb{R}^{d}}\int\limits_{X_2\in\mathbb{R}^{d}}J(X_1,X_2)e^{{-}i(K_1\cdot X_1-K_2\cdot X_2)}\mathrm{d}X_2\mathrm{d}X_1,$$
Eq. (40) becomes
$$\begin{aligned}I_\mathcal{R}=&\sum_{ij}\frac{\epsilon^{i+j}}{(2\pi)^{2d}}\int\limits_{K_1,K_2\in\mathbb{R}^{d}}\hat{J}(K_1,K_2)\mathrm{d}K_1\mathrm{d}K_2\\ &\times\int\limits_{K_1',K_2'\in\mathbb{R}^{d}}\mathbb{E}\big[g_i(K_1,K_1')g_j^{{\ast}}(K_2,K_2')\big]\hat{w}(K_1'-K_2')\mathrm{d}K_1'\mathrm{d}K_2'\\ =&\frac{1}{(2\pi)^{d}}\int\limits_{K_1,K_2\in\mathcal{R}^{d}}\hat{J}(K_1,K_2)\mathcal{H}_{\mathrm{f}}(K_1,K_2)\mathrm{d}K_1\mathrm{d}K_2, \end{aligned}$$
where we define
$$\begin{aligned}\mathcal{H}_{\mathrm{f},ij}&=\frac{1}{(2\pi)^{d}}\int\limits_{K_1'\in\mathbb{R}^{d}}\int\limits_{K_2'\in\mathbb{R}^{d}}\mathbb{E}\big[g_i(K_1,K_1')g_j^{{\ast}}(K_2,K_2')\big]\hat{w}(K_1'-K_2')\mathrm{d}K_1'\mathrm{d}K_2', \end{aligned}$$
and set
$$\mathcal{H}_{\mathrm{f}}=\mathcal{H}_{\mathrm{f},00}+\epsilon^{2}\Big(\mathcal{H}_{\mathrm{f},02}+\mathcal{H}_{\mathrm{f},11}+\mathcal{H}_{\mathrm{f},20}\Big)\,.$$

With the same argument as provided in [4], to produce the optimal beam that provides the maximum light intensity, $\hat {J}$ should be composed of coherent beam, in the sense that

$$\hat{J}(K_1,K_2)=\psi(K_1)\psi^{{\ast}}(K_2),$$
where $\psi (K)$ is the eigenfunction of $\mathcal {H}_{\mathrm {f}}$ associated with the largest eigenvalue. The statement translates the problem to computing $\mathcal {H}_{\mathrm {f}}$ and its eigenfunctions. This can be readily done. Recall (28) we have:
$$\begin{aligned}\mathcal{H}_{\mathrm{f},00}(K_1,K_2)=\frac{1}{(2\pi)^{d}}\hat{G}(K_1,Z)\hat{G}^{{\ast}}(K_2,Z)\hat{w}(K_1-K_2)\,. \end{aligned}$$
To find $\mathcal {H}_{\mathrm {f},11}$, we use (30):
$$\begin{aligned}&\mathcal{H}_{\mathrm{f},11}(K_1,K_2)=\Big(\frac{1}{2\pi}\Big)^{3d}\frac{k^{2}}{4}\int_{z_1=0}^{Z}\int_{z_2=0}^{Z}\hat{G}(K_1,z_1)\hat{G}^{{\ast}}(K_2,z_2)\mathrm{d}z_2\mathrm{d}z_1\\ &\times\int\limits_{K_1',K_2'\in\mathbb{R}^{d}}\hat{G}(K_1',Z-z_1)\hat{G}^{{\ast}}(K_2',Z-z_2)\hat{\Gamma}_V(K_1'-K_1,z_1,K_2'-K_2,z_2)\hat{w}(K_1'-K_2')\mathrm{d}K_1'\mathrm{d}K_2' \end{aligned}.$$
Here $\hat {\Gamma }_V$ is the covariance of the random medium in Fourier space, given by
$$\hat{\Gamma}_V(K_1,z_1,K_2,z_2)=\mathbb{E}\big[\hat{V}(K_1,z_1)\hat{V}^{{\ast}}(K_2,z_2)\big]\,.$$

Similarly, using (31) we have

$$\begin{aligned}&\mathcal{H}_{\mathrm{f},02}(K_1,K_2)={-}\Big(\frac{1}{2\pi}\Big)^{3d}\frac{k^{2}}{8}\hat{G}(K_1,Z)\int\limits_{z_1=0}^{Z}\int\limits_{z_2=0}^{Z}\hat{G}^{{\ast}}(K_2,z_2)\mathrm{d}z_2\mathrm{d}z_1\\ &\times\int\limits_{K_1',K_2'\in\mathbb{R}^{d}}\hat{G}^{{\ast}}(K_2',z_1-z_2)\hat{G}^{{\ast}}(K_1',Z-z_1)\hat{\Gamma}_V(K_2-K_2',z_2,K_1'-K_2',z_1)\hat{w}(K_1-K_1')\mathrm{d}K_1'\mathrm{d}K_2' \end{aligned}.$$
Similarly:
$$\mathcal{H}_{\mathrm{f},20}(K_1,K_2)=\mathcal{H}_{\mathrm{f},02}^{{\ast}}(K_1,K_2).$$

We end our discussion by pointing out that even before any assumptions on the turbulence gets incorporated, the Fourier domain quantity $\mathcal {H}_{\mathrm {f},ij}$ is already easier than that of the physical-space quantity $\mathcal {H}_{ij}$. Indeed, as presented in (45–48), the $\mathcal {H}_{\mathrm {f},ij}$ quantities are four-folded integrals once $\hat {\Gamma }$ is given. The extra integral in the $\mathcal {H}_{ij}$ formula in (36) in $X\in \mathcal {R}$ is absorbed in $\hat {w}$ and has been completed analytically. More specifically, if $d=2$, for every fixed $(K_1,K_2)$, the calculation of $\mathcal {H}_{\mathrm {f},ij}$ is a $6$-dimensional integration, as compared to $8$-dimensional as shown in (36).

5. Properties of $\mathcal {H}$

5.1 Relation between $\mathcal {H}$ and $\mathcal {H}_{\mathrm {f}}$

While $\mathcal {H}$ and $\mathcal {H}_{\mathrm {f}}$ were computed above from different perspectives (in physical space and Fourier space), they are related. To see the connection between them, note that the total average intensity can be written in terms of the mutual intensity function as

$$\begin{aligned}I&=\frac{1}{(2\pi)^{d}}\int\limits_{K_1\in\mathbb{R}^{d}}\int\limits_{K_2\in\mathbb{R}^{d}}\hat{J}(K_1,K_2)\mathcal{H}_{\mathrm{f}}(K_1,K_2)\mathrm{d}K_1\mathrm{d}K_2\\ &=\frac{1}{(2\pi)^{d}}\int\limits_{K_1\in\mathbb{R}^{d}}\int\limits_{K_2\in\mathbb{R}^{d}}\mathcal{H}_{\mathrm{f}}(K_1,K_2)\mathrm{d}K_1\mathrm{d}K_2\int\limits_{X_1\in\mathcal{A}}\int\limits_{X_2\in\mathcal{A}}J(X_1,X_2)e^{{-}i(K_1\cdot X_1-K_2\cdot X_2)}\mathrm{d}X_1\mathrm{d}X_2\\ &=\frac{1}{(2\pi)^{d}}\int\limits_{X_1\in\mathcal{A}}\int\limits_{X_2\in\mathcal{A}}J(X_1,X_2)\mathrm{d}X_1\mathrm{d}X_2\int\limits_{K_1\in\mathbb{R}^{d}}\int\limits_{K_2\in\mathbb{R}^{d}}\mathcal{H}_{\mathrm{f}}(K_1,K_2)e^{{-}i(K_1\cdot X_1-K_2\cdot X_2)}\mathrm{d}K_1\mathrm{d}K_2 \end{aligned}$$
which shows that
$$\mathcal{H}(X_1,X_2)=\frac{1}{(2\pi)^{d}}\int\limits_{K_1\in\mathbb{R}^{d}}\int\limits_{K_2\in\mathbb{R}^{d}}\mathcal{H}_{\mathrm{f}}(K_1,K_2)e^{{-}i(K_1\cdot X_1-K_2\cdot X_2)}\mathrm{d}K_1\mathrm{d}K_2\, ,$$
so that $\mathcal {H}$ and $\mathcal {H}_{\mathrm {f}}$ are themselves related by a Fourier transformation.

5.2 Energy conservation

We now show that, by keeping terms up to $O(\epsilon ^{2})$ in the expansion, energy is conserved. This is an important reason for keeping the $O(\epsilon ^{2})$ terms in addition to the $O(\epsilon )$ terms. Indeed, one criticism over employing perturbation theory in wave propagation is that it often loses preservation of some crucial physical properties, such as energy conservation [33]. However, energy conservation is retained here.

To see so, let the receiver occupy the entire space (i.e., $\mathcal {R}=\mathbb {R}^{d}$), in which case the window function becomes:

$$\hat{w}(K_1-K_2)=\int\limits_{X'\in\mathcal{R}^{d}}\mathrm{d}X'e^{i(K_1-K_2)\cdot X'}=(2\pi)^{d}\delta(K_1-K_2)\,.$$
Inserting this into (44) and (45–48), we have
$$\mathcal{H}_{\mathrm{f},00}(K_1,K_2)=\delta(K_1-K_2), \quad \text{and}\quad \mathcal{H}_{\mathrm{f},11}+\mathcal{H}_{\mathrm{f},02}+\mathcal{H}_{\mathrm{f},20}=0,$$
so that
$$\begin{aligned}I&=\frac{1}{(2\pi)^{d}}\int\limits_{K_1\in\mathbb{R}^{d}}\int\limits_{K_1\in\mathbb{R}^{d}}\hat{J}(K_1,K_2)\delta(K_1-K_2)\mathrm{d}K_1\mathrm{d}K_2\\ &=\frac{1}{(2\pi)^{d}}\int\limits_{K_1\in\mathbb{R}^{d}}\hat{J}(K_1,K_1)\mathrm{d}K_1=\int\limits_{X\in\mathcal{A}}J(X,X)\mathrm{d}X=I_0\,. \end{aligned}$$
We note that the energy conservation holds independent of $\epsilon$. This means that, even when $\epsilon$ and the associated approximation errors are relatively large, energy conservation still holds true and provides a certain level of physical realism.

6. Some optional simplifications

One advantage of utilizing the Fourier-space presentation is that it makes it much easier to incorporate the classical assumptions on atmosphere turbulence since these assumptions are typically specified in the Fourier space. This would allow us to further reduce the computation.

In this subsection, we describe several common assumptions and the resulting simplifications to the formulas for the $\mathcal {H}_{\mathrm {f},ij}$ quantities. Then, in a later section, a selection of these cases will be investigated via numerical calculations.

6.1 Homogeneous-statistics assumption

One classical assumption is the homogeneous (or stationary) property. This is to assume the covariance $\Gamma _V$ of the medium has the structure of

$$\Gamma_V(X_1,z_1,X_2,z_2)=f(X_1-X_2,z_1-z_2)\, .$$
In this case, we also have
$$\hat{\Gamma}_V(K_1,z_1,K_2,z_2)=(2\pi)^{d}F(K_1,z_1-z_2)\delta(K_1-K_2),$$
where $F(K,z)=\hat {f}(X,z)$ is the Fourier transform. This newly induced $\delta$ function in the $K$ domain helps to eliminate one-fold of integration. For example, when inserted into (45), we have:
$$\begin{aligned}\mathcal{H}_{\mathrm{f},11}(K_1,K_2)&=\frac{k^{2}}{4(2\pi)^{d}}\mathcal{H}_{\mathrm{f},00}(K_1,K_2)\int_{z_1=0}^{Z}\int_{z_2=0}^{Z}\mathrm{d}z_2\mathrm{d}z_1\\ &\times\int\limits_{K'\in\mathbb{R}^{d}}F(K',z_1-z_2)\hat{G}(K',z_1-z_2)\exp\Big(-\frac{i}{k}K'\cdot(K_1z_1-K_2z_2)\Big)\mathrm{d}K'\\ &=\frac{k^{2}}{4(2\pi)^{d}}\mathcal{H}_{\mathrm{f},00}(K_1,K_2)\int_{z_1=0}^{Z}\int_{z_2=0}^{Z}\widehat{F\hat{G}}\big((K_1z_1-K_2z_2)/k,z_1-z_2\big)\mathrm{d}z_2\mathrm{d}z_1\,. \end{aligned}$$
We rewrote the innermost integral in terms of a Fourier transform in the last equation.

Similarly,

$$\mathcal{H}_{\mathrm{f},02}(K_1,K_2)={-}\frac{k^{2}}{8(2\pi)^{d}}\mathcal{H}_{\mathrm{f},00}(K_1,K_2)\int_{z_1=0}^{Z}\int_{z_2=0}^{Z}\widehat{F\hat{G}}\big(K_2(z_2-z_1)/k,z_1-z_2\big)\mathrm{d}z_2\mathrm{d}z_1\,.$$
and $\mathcal {H}_{\mathrm {f},20}(K_1,K_2)=\mathcal {H}_{\mathrm {f},02}^{\ast }(K_1,K_2)$. We note that this has further simplified the computation to a $4$-dimensional integral, with two dimensions absorbed into the Fourier transform when $d=2$. Moreover, the Fourier transform component, $\widehat {F\hat {G}}$, though being a $2$-dimensional integral, can be calculated through FFT, which further reduces an $N^{2}$ computational cost to $N\log {N}$.

6.2 Decorrelation-in-$z$ assumption

Another assumption that one might choose to make, in addition to the earlier assumption of homogeneity, is regarding the decay rate of the covariance in $z$. This assumption gets widely used, for instance, in [30]. It states that for some characteristic length $\delta _z\ll Z$,

$$F(K,z)\sim 0, \quad \mbox{for}\quad |z|>\delta_z\,.$$
Note that, for comparison, this is a slightly relaxed assumption compared to the Markov approximation. If the Markovian approximation is assumed in the $z$ direction, then the random medium at every $z$ point is statistically independent. The assumption in (58), in contrast, allows a non-trivial correlation length in the $z$ direction, up to the length scale $\delta _z$.

To utilize this additional assumption, it is helpful to define the transformation

$$\eta=\frac{z_1+z_2}{2}, \quad \mu=z_1-z_2.$$
We then perform the change of variable to have
$$\begin{aligned}\mathcal{H}_{\mathrm{f},11}(K_1,K_2)=&\frac{k^{2}}{4(2\pi)^{d}}\mathcal{H}_{\mathrm{f},00}(K_1,K_2)\int\limits_{K'\in\mathbb{R}^{d}}\mathrm{d}K'\int_{\eta=0}^{Z/2} \exp\Big(-\frac{i}{k}K'\cdot(K_1-K_2)\eta\Big)\mathrm{d}\eta\\ &\times\int_{\mu={-}2\eta}^{2\eta}F(K',\mu)\hat{G}(K',\mu)\exp\Big(-\frac{i}{2k}K'\cdot(K_1+K_2)\mu\Big)\mathrm{d}\mu\\ &+\frac{k^{2}}{4(2\pi)^{d}}\mathcal{H}_{\mathrm{f},00}(K_1,K_2)\int\limits_{K'\in\mathbb{R}^{d}}\mathrm{d}K'\int_{\eta=Z/2}^{Z} \exp\Big(-\frac{i}{k}K'\cdot(K_1-K_2)\eta\Big)\mathrm{d}\eta\\ &\times\int_{\mu={-}2(Z-\eta)}^{2(Z-\eta)}F(K',\mu)\hat{G}(K',\mu)\exp\Big(-\frac{i}{2k}K'\cdot(K_1+K_2)\mu\Big)\mathrm{d}\mu\\ \approx&\frac{k^{2}}{4(2\pi)^{d}}\mathcal{H}_{\mathrm{f},00}(K_1,K_2)\int\limits_{K'\in\mathbb{R}^{d}}\mathrm{d}K'\int_{\eta=0}^{Z} \exp\Big(-\frac{i}{k}K'\cdot(K_1-K_2)\eta\Big)\mathrm{d}\eta\\ &\times\int_{\mu\in\mathbb{R}}F(K',\mu)\hat{G}(K',\mu)\exp\Big(-\frac{i}{2k}K'\cdot(K_1+K_2)\mu\Big)\mathrm{d}\mu \end{aligned}$$
where the second estimate comes from (58). The estimate includes the extra integration area $\int ^{\delta _z/2}_{0}\int _{2\eta }^{\delta _z}\mathrm {d}\mu \mathrm {d}\eta \sim \mathcal {O}(\delta _z^{2})\to 0$ for small $\delta _z$.

6.3 Small-length-scale cutoff assumption

In addition to the assumption from (58), we now suppose the turbulence is characterized by the Kolmogorov model and with $l_0$ and $L_0$ being the inner and outer scales of turbulence, respectively. This means, if we write the power spectral density, the complete Fourier transform of the covariance function, given by:

$$\begin{aligned}\Phi(K,K_z)&=\epsilon^{2} \int\limits_{z\in\mathbb{R}}\int\limits\limits_{X\in\mathbb{R}^{d}}f(X,z)e^{{-}i(K\cdot X+K_zz)}\mathrm{d}X\mathrm{d}z\\ &=\epsilon^{2} \int\limits_{z\in\mathbb{R}}F(K,z)e^{{-}K_zz}\mathrm{d}z \end{aligned}.$$
There are a few assumptions in place:

Firstly, we assume $\Phi$ is a power law in the inertial subrange $1/L_0\ll |(K,K_z)|=\sqrt {\|K\|^{2}+{K_z}^{2}}\ll 1/l_0$. Letting $L_0=\infty$, we run the von Karman spectrum multiplied by a Gaussian factor approximation outside the inertial subrange to ensure that the power spectrum decays rapidly for wavenumbers larger than $1/l_0$ [2]. These assumptions makes $F(K,z)$ negligible for $\|K\|\gg 1/l_0$.

Secondly, as noted in [30], we also assume $F(K, z)$ is negligible also for $\|K\|z>1$. These assumptions together suggest that $F(K,z)$ is not negligible only when:

$$\frac{\|K\|^{2}z}{2}<\frac{\|K\|}{2}<\frac{1}{2l_0}\,.$$

However, in this region, we recall the definition of $\hat {G}$ in (16), we will see that $\hat {G}\sim 1$. This can be seen evaluating

$$\hat{G} =\exp\{-\frac{i z}{2k}\|K\|^{2}\}\sim \exp(-\frac{i}{2kl_0})\sim\exp\{0\}=1,$$
where we used the assumption that $1/(2kl_0)\ll 1$. This assumption is realistic in the atmosphere for waves in the optical/IR regime according to [2]. Plug this approximation back into Eq. (60), we have:
$$\begin{aligned}&\mathcal{H}_{\mathrm{f},11}(K_1,K_2)\\ \approx&\frac{Zk^{2}}{4(2\pi)^{d}}\mathcal{H}_{\mathrm{f},00}(K_1,K_2)\int\limits_{K'\in\mathbb{R}^{d}}\exp\Big(-\frac{i}{2k}K'\cdot(K_1-K_2)Z\Big)sinc\Big(\frac{K'}{2k}\cdot(K_1-K_2)Z\Big)\mathrm{d}K'\\ &\times\int_{\mu\in\mathbb{R}}F(K',\mu)\exp\Big(-\frac{i}{2k}K'\cdot(K_1+K_2)\mu\Big)\mathrm{d}\mu\\ =&\frac{Zk^{2}}{4(2\pi)^{d}\epsilon^{2}}\mathcal{H}_{\mathrm{f},00}(K_1,K_2)\int\limits_{K'\in\mathbb{R}^{d}}\exp\Big(-\frac{i}{2k}K'\cdot(K_1-K_2)Z\Big)sinc\Big(\frac{K'}{2k}\cdot(K_1-K_2)Z\Big)\\ &\times\Phi(K',K'\cdot(K_1+K_2)/2k)\mathrm{d}K'\,. \end{aligned}$$
Similarly, we have
$$\mathcal{H}_{\mathrm{f},02}(K_1,K_2)\approx{-}\frac{Zk^{2}}{8(2\pi)^{d}\epsilon^{2}}\mathcal{H}_{\mathrm{f},00}(K_1,K_2)\int\limits_{K'\in\mathbb{R}^{d}}\Phi(K',-K'\cdot K_2/k)\mathrm{d}K'$$
and $\mathcal {H}_{\mathrm {f},20}=\mathcal {H}^{\ast }_{\mathrm {f},02}$.

6.4 Markov approximation

One further simplifies the calculation when Markov approximation is imposed.

Let the covariance $f$ take the form of

$$f(x,z)=f_1(x)f_2(z),$$
then Eqs. (56) and (57) get re-written as
$$\mathcal{H}_{\mathrm{f},11}=\frac{k^{2}}{4(2\pi)^{d}}\mathcal{H}_{\mathrm{f},00}(K_1,K_2)\int_{z_1=0}^{Z}\int_{z_2=0}^{Z}f_2(z_1-z_2)\widehat{F_1\hat{G}({\cdot}}, z_1-z_2)\big((K_1z_1-K_2z_2)/k\big)\mathrm{d}z_2\mathrm{d}z_1,$$
and
$$\mathcal{H}_{\mathrm{f},02}(K_1,K_2)={-}\frac{k^{2}}{8(2\pi)^{d}}\mathcal{H}_{\mathrm{f},00}(K_1,K_2)\int_{z_1=0}^{Z}\int_{z_2=0}^{Z}f_2(z_1-z_2)\widehat{F_1\hat{G}({\cdot}}, z_1-z_2)\big(K_2(z_2-z_1)/k\big)\mathrm{d}z_2\mathrm{d}z_1,$$
where $F_1=\hat {f}_1$. Markov approximation means that $f_2=\delta (z)$, meaning the turbulence for every $z$ point is completely independent, then:
$$f(x,z)=f_1(x)\delta(z),$$
which further reduces the computation of $\mathcal {H}_{\mathrm {f},11}$ and $\mathcal {H}_{\mathrm {f},02}$ to
$$\mathcal{H}_{\mathrm{f},11}=\frac{Zk^{2}}{4}\mathcal{H}_{\mathrm{f},00}(K_1,K_2)\int_{y=0}^{1}f_1\big((K_2-K_1)yZ/k\big)\mathrm{d}y,$$
and
$$\mathcal{H}_{\mathrm{f},02}(K_1,K_2)={-}\frac{Zk^{2}}{8}f_1(0)\mathcal{H}_{\mathrm{f},00}(K_1,K_2),$$
with $\mathcal {H}_{\mathrm {f},20}(K_1,K_2)=\mathcal {H}_{\mathrm {f},02}^{\ast }(K_1,K_2)$. Note that (68) is a single integral and is computationally cheap.

We should note, however, when the assumption is this strong, the computation on the physical domain is similarly simple. Indeed, define $\eta =\frac {X_1'+X_2'}{2}$ and $\mu =X_1'-X_2'$, (35) becomes:

$$\begin{aligned}\mathcal{H}_{11}&=\frac{k^{2}}{4}\int\limits_{z_1=0}^{Z}\Big(\frac{k}{2\pi(Z-z_1)}\Big)^{d}G(X_1,z_1)G^{{\ast}}(X_2,z_1)\mathrm{d}z_1\int\limits_{\mu\in\mathbb{R}^{d}}f(\mu)\exp\big(-\frac{ik\mu(X_1+X_2)}{2z_1}\big)\mathrm{d}\mu\\ &\times\int\limits_{X'\in\mathcal{R}}\exp\big(-\frac{ikX'\cdot\mu}{Z-z_1}\big)\mathrm{d}X'\int\limits_{\eta\in\mathbb{R}^{d}}\exp\Big(ik\eta\big(\frac{Z\mu-(X_1-X_2)(Z-z_1)}{z_1(Z-z_1)}\big)\Big)\mathrm{d}\eta\\ &=\frac{Zk^{2}}{4}\mathcal{H}_{00}(X_1,X_2)\int\limits_{y=0}^{1}f_1\big((X_1-X_2)y\big)\mathrm{d}y\,. \end{aligned}$$
Similarly, (36) becomes:
$$\mathcal{H}_{02}={-}\frac{Zk^{2}}{8}f_1(0)\mathcal{H}_{00}(X_1,X_2)\,.$$
This gives a one-fold integral and is numerically easy as well. This means when Markov approximation holds true, the computation on the physical domain is similarly easy with its counterpart from Fourier domain.

7. Numerical examples

In this section, we present numerical examples to demonstrate the proposed methods for computing the optimal beam and its associated mutual intensity function.

To set up the computation, we assume a 2D domain, with the transmitter region being $X\in [-r,r]$ with $r=0.05$ m, and the receiver being located at $Z=3000$ m. The wave frequency is set to be $k=2\pi \times 10^{6}$ rad/m. For discretization we use a mesh size of $\Delta x=0.001$m. This gives 100 grid points in the $x$ direction within the transmitter region. The computational domain is $(x,z)\in [-L,L]\times [0,Z]$, with $L=1$ m. The domain half-width $L$ is set to be much larger than $r$ to ensure that the waves at the computational boundary are negligible. We use periodic boundary conditions, meaning $A(-L,z)=A(L,z)$, and $\partial _xA(-L,z)=\partial _xA(L,z)$ in the $x$ direction. From the perspective of Fourier space, the wave numbers range from $[-\pi /\Delta x,\pi /\Delta x)$ with $\Delta K=\pi /L$. This leads to 2000 grid points in Fourier space.

In what follows, in subsection 7.1, we compare $\mathcal {H}$ computed in physical space, and $\mathcal {H}_{\mathrm {f}}$ computed in Fourier space. It will be seen that the numerical results of $\mathcal {H}$ and $\mathcal {H}_{\mathrm {f}}$ are on top of each other, which suggests a fine enough resolution is being used. In subsection 7.2, we present the shapes of optimal beams when different assumptions are incorporated. The results suggest that the small-length-scale cutoff assumption, which significantly reduces the numerical cost, provides accurate approximations in the calculation of $\mathcal {H}$. In subsection 7.3, we show the behavior of the optimal beams under different levels of the strength of the turbulence.

7.1 Validation of numerical resolution

To get started, we first numerically verify the equivalence between the computation provided from the physical space (33) and that from the Fourier space (44). We assume $\epsilon =0$ so there is no turbulence in the medium. For a receiver region of $X\in \mathcal {R}=[-R,R]$ with $R$ being the radius of the receiver region, we rewrite (33) to be:

$$\mathcal{H}_{00}(X_1,X_2)=\Big(\frac{Rk}{\pi Z}\Big)\exp\Big(\frac{ik}{2Z}(X_1^{2}-X_2^{2})\Big)sinc\Big(\frac{Rk}{Z}(X_1-X_2)\Big),\quad (X_1,X_2)\in[{-}r,r]\,.$$
Similarly, for this particular setup, we rewrite (44) to be:
$$\mathcal{H}_{\mathrm{f},00}(K_1,K_2)=\frac{R}{\pi}\exp\Big(-\frac{iZ}{2k}(K_1^{2}-K_2^{2})\Big)sinc\big(R(K_1-K_2)\big)\,.$$
In Fig. 1, we plot $\mathcal {H}$ with different $\mathcal {R}$, the size of the receiver. As shown in the plots, as $\mathcal {R}$ increases, $\mathcal {H}$ becomes closer and closer to the identity matrix, with more and more eigenvalues closer to $1$. Physically this means that all modes from the transmitter arrive at the receiver with the intensity preserved.

 figure: Fig. 1.

Fig. 1. $\mathcal {H}$ computed using (72) when the receiver size $R$ is $0.01$, $0.03$, $0.06$, and $0.09$ m. The transmitter size is fixed at $r=0.05$ m.

Download Full Size | PDF

Computing the optimal beams using (72) and (73) should agree. This is shown in Fig. 2 and Fig. 3. In particular, in Fig. 2 we demonstrate the intensity of the optimal beam for different $R$, and the agreement of the first nine eigenvalues. We should note that the optimal beams give higher intensity than the focused beam (using initial data $A(X,z=0)=\phi (X)=e^{-\frac {ik}{2Z}X^{2}}$ as a complex Gaussian) for all $R$, suggesting the focused beams are not optimal. One interesting phenomenon to be observed here is that the optimal beam achieves the full intensity when the receiver size is only about $0.05$m. This is the same size as the transmitter, indicating that the beam divergence is small. This observation resonates with the calculation shown in Fig. 4 of [3] where it suggests the full intensity can be captured when $R\sim \frac {2Z}{kr}$, which agrees with our computation. In Fig. 3, we plot the profile of the beams. Once again, the calculation given on the physical space and that given on the Fourier space agree with each other. These agreements suggest the numerical resolution is fine enough for the numerical experiments to be trusted.

 figure: Fig. 2.

Fig. 2. In the panel on the left, we show the comparison of total intensity over various receiver sizes, and the focused beam solution. Receiver size is in meters. In the panel on the right, we show the agreement of the first nine eigenvalues given by $\mathcal {H}$ and $\mathcal {H}_\mathrm {f}$, taking $R=0.09$m.

Download Full Size | PDF

 figure: Fig. 3.

Fig. 3. Eigenfunctions provided by computing $\mathcal {H}$ (blue line) and $\mathcal {H}_{\mathrm {f}}$ (red line) are on top of each other. The left and the middle panel show the real and imaginary parts of the first six engenvectors, and the panel on the right shows the light intensity received at the receiver with the $R = 0.09$m for each of the six eigenfunctions.

Download Full Size | PDF

 figure: Fig. 4.

Fig. 4. The plot in blue is the relative error in Fourier space vs physical space calculations for $\mathcal {H}$ in a uniform medium (with respect to the Frobenious norm). The red plot shows the relative difference in the presence and absence of turbulence from Fourier space calculations. $\epsilon$ is fixed at $5\times 10^{-8}$ ($R$ is in meters)

Download Full Size | PDF

In Fig. 4, we plot the difference of $\mathcal {H}$ computed using the brute-force calculation (72), and that computed using the Fourier transform (50) from a simplified $\mathcal {H}_\mathrm {f}$ (73). For all $R$ this error is significantly smaller than the difference between $\mathcal {H}$ computed for the different $\epsilon$ ($\epsilon =0$ vs. $\epsilon = 5\times 10^{-8}$, where further details of the $\epsilon >0$ case are described below). This is further evidence that the numerical errors are small compared to the changes in the solutions brought about by turbulent fluctuations.

7.2 Shapes of optimal beams, and cost comparisons

In the second examples, we investigate cases with turbulent fluctuations in the medium, and we utilize the homogeneous assumption discussed in Section 6.1. We use the same setup as above, and assume $\epsilon =5\times 10^{-8}$ (slightly larger than $0.3/k$). As a start, we assume $\Gamma$ to be a Gaussian function:

$$\Gamma(X_1,z_1,X_2,z_2)=\exp\Big(-\frac{(X_1-X_2)^{2}+(z_1-z_2)^{2}}{l_0^{2}}\Big)$$
with the correlation length $l_0=0.1$ m [2]. Under this assumption, we can compute $F$ to be:
$$F(K,z)=\sqrt{\pi l_0^{2}}\exp\Big(-\frac{l_0^{2}K^{2}}{4}\Big)\exp\Big(-\frac{z^{2}}{l_0^{2}}\Big)\,.$$
When the homogeneous assumption holds true, we take the formulas in (56) and (57) and insert
$$\widehat{F\hat{G}}(K',z_1-z_2)=\sqrt{\pi l_0^{2}}\sqrt{\frac{\pi}{s(z_1-z_2)}} \exp\Big(-\frac{(z_1-z_2)^{2}}{l_0^{2}}\Big)\exp\Big(-\frac{K'^{2}}{4s(z_1-z_2)}\Big),$$
where $s(z)=\frac {l_0^{2}}{4}+\frac {iz}{2k}$. Note that $\exp [-(z_1-z_2)^{2}/l_0^{2}]$ is negligible when $|z_1-z_2|>3l_0$, so one can perform integration in $z_1$ and $z_2$ in the domain of $(z_1,z_2)\in (0,Z]\times (z_1-3l_0,z_1)$ in the region of $z_2<z_1$. The same simplification can be used for $z_2>z_1$. With turbulence added in this way, we discover that the optimal beam is slightly wider than it is in the uniform-medium case, as illustrated in Fig. 5 (left panel) computed using $R=0.05$m.

 figure: Fig. 5.

Fig. 5. In the panel on the left, we show the initial intensity profile of optimal beam using $R=0.05$m. In the panel on the right, we show the comparison of total intensity at the receiver using optimal beams and focused beams under the homogeneous random medium assumption and small-length-scale cutoff assumption ($R$ is in meters).

Download Full Size | PDF

We furthermore test the accuracy of the small-length-scale cutoff assumption (see Section 6.3). For this purpose, we use (62), (63) with

$$\Phi(K,K_z)=\epsilon^{2}\pi l_0^{2}\exp\Big(-\frac{l_0^{2}(K^{2}+K_z^{2})}{4}\Big)\,.$$
Noting that $\Phi (K,K_z)$ is negligible when $\|K\|>6/l_0$, we set $K\in [-K_{\max },K_{\max }]$ with $\Delta K=1/4l_0$ and $K_{\max }=6/l_0$. In Fig. 5 (right panel), we test the light intensity given by the optimal beam using different assumptions in the calculation. The computation generated by using small-length-scale cutoff assumption agrees very well with that generated using the homogeneous assumption only, for all values of $R$. One should note, however, the computation using the homogeneous random medium assumption, as shown in Eq. (56), uses $2$-dimensional integral (for $d=1$ setup here), takes around $16$ minutes, while the same computation took around $3$ seconds using the small-length-scale cutoff assumption, where the computation is $1$-folded integral, suggested in (61). This holds true for every entry of $\mathcal {H}_\mathrm {f}$, and thus brings a significant savings in computation.

Also shown in Fig. 5 is a comparison of the optimal beam and the focused beam. In all cases, the focused beam gives weaker light intensity at the receiver in comparison to the optimal beam. Also, while one might think that the optimization result is a trivial result because it has an approximately Gaussian profile of intensity, it is important to note that the optimal beam is complex-valued and the phase information is crucial. The complex-valued beam is not itself a Gaussian profile. See Fig. 3 for an example of a non-Gaussian, complex-valued beam profile which has a Gaussian profile of intensity.

Figure 5 is evidence that the optimal beam has a beam divergence that is small. More specifically, note that the optimal beam has an intensity of approximately 1 after traveling a distance of $Z=3000$ m, from a transmitter of size $r=0.05$ m to a receiver of the same size ($R=0.05$ m). Hence, beam divergence must be small in order to allow the full intensity of the beam to reach such a small receiver. Such a phenomenon has been seen in previous calculations of optimal beams, in the case of a uniform non-turbulent medium and a phase screen model of turbulence [3]. In Fig. 5, this phenomenon is also seen in the case of the turbulence setup of the present paper.

7.3 Sensitivity studies

In this subsection, we study the relation between the optimal beam and the intensity of the turbulence. We would like to see the changes in the optimal beam that arise as greater turbulent fluctuations are introduced. As seen in Fig. 6, the total intensity of optimal beam drops as $\epsilon$ increases, for all choices of $R$. For small $R$, this drop is significant. The profile of the optimal beam also changes as $\epsilon$ changes, and the change is most prominent for small $R$, as shown in Fig. 7.

 figure: Fig. 6.

Fig. 6. Optimal intensity at the receiver as $\epsilon$ increases ($R$ is in meters)

Download Full Size | PDF

 figure: Fig. 7.

Fig. 7. Initial intensity profile of optimal beams for different receiver sizes for $\epsilon =5\times 10^{-8}$ ($R$ is in meters)

Download Full Size | PDF

8. Concluding discussion

In this paper, we computed profiles of optimal beams that achieve the highest intensity at the receiver. Mathematically this amounts to solving for eigenfunctions of $\mathcal {H}$. In most realistic settings, $\mathcal {H}$ is numerically infeasible to compute, with each entry calling for an $8$-fold integral. We proposed to convert the calculation to the Fourier domain where assumptions on the turbulent medium can be naturally incorporated. By introducing assumptions of spatially homogeneous statistics of the random medium, small-length-scale cutoff assumption, and Markov assumption, the $8$-fold integral is replaced by $6$-fold, $2$-fold and $1$-fold integrals respectively, with the numerical cost significantly reduced. This research generalizes the existing literature that concerns mostly the Markov approximation, or special classes of mutual intensity functions, and now allows for a general beam profile in a much more general turbulent medium structure. The methods proposed here point toward the possibility of computing general optimal beams.

The numerical examples here suggest that optimal beams can have nearly the full intensity transmitted, and small beam divergence. Similar results had also been seen in past studies under different numerical experimental setups [3]. These types of results show the influence of initial beam properties on the downstream characteristics of the beam.

Funding

Office of Naval Research (N00014-21-1-2119, N00014-21-1-2140).

Acknowledgments

The authors thank Svetlana Avramov-Zamurovic for helpful comments and discussion. The research of Q.L. is partially supported by Office of Naval Research (ONR) grant N00014-21-1-2140, and the research of A.N. and S.N.S. is partially supported by ONR grant N00014-21-1-2119.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data and code underlying the results presented in this paper are available in [34].

References

1. J. W. Strohbehn, Laser beam propagation in the atmosphere (Springer-Verlag, New York, 1978).

2. L. C. Andrews and R. L. Phillips, Laser Beam Propagation through Random Media (SPIE Press, 2005).

3. T. J. Schulz, “Iterative transform algorithm for the computation of optimal beams,” J. Opt. Soc. Am. A 21(10), 1970–1974 (2004). [CrossRef]  

4. T. J. Schulz, “Optimal beams for propagation through random media,” Opt. Lett. 30(10), 1093–1095 (2005). [CrossRef]  

5. B. Liu, “Optimal beam forming for laser beam propagation through random media,” Ph. D. Thesis, Michigan Technological University (2006).

6. J. Zhou, J. Wu, and Q. Hu, “Optimal transmission modes under atmosphere turbulence with transmitter/receiver aperture size constraint,” Opt. Express 26(25), 33333–33348 (2018). [CrossRef]  

7. D. Slepian, “Prolate spheroidal wave functions, fourier analysis and uncertainty—iv: extensions to many dimensions; generalized prolate spheroidal functions,” Bell Syst. Tech. J. 43(6), 3009–3057 (1964). [CrossRef]  

8. A. Belmonte and J. M. Kahn, “Approaching fundamental limits to free-space communication through atmospheric turbulence,” in Broadband Access Communication Technologies XII, vol. 10559 (SPIE, 2018), pp. 70–76.

9. J. Shapiro, S. Guha, and B. Erkmen, “Ultimate channel capacity of free-space optical communications,” J. Opt. Netw. 4(8), 501–516 (2005). [CrossRef]  

10. D. G. Voelz and X. Xiao, “Metric for optimizing spatially partially coherent beams for propagation through turbulence,” Opt. Eng. 48(3), 036001 (2009). [CrossRef]  

11. Y. Cai and S. He, “Average intensity and spreading of an elliptical gaussian beam propagating in a turbulent atmosphere,” Opt. Lett. 31(5), 568–570 (2006). [CrossRef]  

12. P. Polynkin, A. Peleg, L. Klein, T. Rhoadarmer, and J. Moloney, “Optimized multiemitter beams for free-space optical communications through turbulent atmosphere,” Opt. Lett. 32(8), 885–887 (2007). [CrossRef]  

13. X. Qian, W. Zhu, and R. Rao, “Numerical investigation on propagation effects of pseudo-partially coherent gaussian schell-model beams in atmospheric turbulence,” Opt. Express 17(5), 3782–3791 (2009). [CrossRef]  

14. O. Korotkova, S. Avramov-Zamurovic, C. Nelson, R. Malek-Madani, Y. Gu, and G. Gbur, “Scintillation reduction in multi-gaussian schell-model beams propagating in atmospheric turbulence,” in Laser Communication and Propagation through the Atmosphere and Oceans III, vol. 9224 (SPIE, 2014), pp. 190–196.

15. L. Borcea, J. Garnier, and K. Sølna, “Multimode communication through the turbulent atmosphere,” J. Opt. Soc. Am. A 37(5), 720–730 (2020). [CrossRef]  

16. A. Dogariu and S. Amarande, “Propagation of partially coherent beams: turbulence-induced degradation,” Opt. Lett. 28(1), 10–12 (2003). [CrossRef]  

17. B. Chen, Z. Chen, and J. Pu, “Propagation of partially coherent bessel-gaussian beams in turbulent atmosphere,” Opt. Laser Technol. 40(6), 820–827 (2008). [CrossRef]  

18. X. Ji, X. Chen, and B. Lü, “Spreading and directionality of partially coherent hermite-gaussian beams propagating through atmospheric turbulence,” J. Opt. Soc. Am. A 25(1), 21–28 (2008). [CrossRef]  

19. G. Gbur and E. Wolf, “Spreading of partially coherent beams in random media,” J. Opt. Soc. Am. A 19(8), 1592–1598 (2002). [CrossRef]  

20. O. Korotkova, L. C. Andrews, and R. L. Phillips, “Model for a partially coherent gaussian beam in atmospheric turbulence with application in lasercom,” Opt. Eng. 43(2), 330–341 (2004). [CrossRef]  

21. J. C. Ricklin and F. M. Davidson, “Atmospheric optical communication with a gaussian schell beam,” J. Opt. Soc. Am. A 20(5), 856–866 (2003). [CrossRef]  

22. G. Gbur, “Partially coherent beam propagation in atmospheric turbulence,” J. Opt. Soc. Am. A 31(9), 2038–2045 (2014). [CrossRef]  

23. E. Wolf, “New theory of partial coherence in the space–frequency domain. part i: spectra and cross spectra of steady-state sources,” J. Opt. Soc. Am. 72(3), 343–351 (1982). [CrossRef]  

24. F. Gori, “Mode propagation of the field generated by collett-wolf schell-model sources,” Opt. Commun. 46(3-4), 149–154 (1983). [CrossRef]  

25. T. Shirai, A. Dogariu, and E. Wolf, “Mode analysis of spreading of partially coherent beams propagating through atmospheric turbulence,” J. Opt. Soc. Am. A 20(6), 1094–1102 (2003). [CrossRef]  

26. T. Habashy, A. T. Friberg, and E. Wolf, “Application of the coherent-mode representation to a class of inverse source problems,” Inverse problems 13(1), 47–61 (1997). [CrossRef]  

27. F. D. Tappert, “The parabolic approximation method,” in Wave Propagation and Underwater Acoustics, (Springer, 1977), Lecture Notes in Physics, vol 70., pp. 224–287.

28. A. C. Radder, “On the parabolic equation method for water-wave propagation,” J. Fluid Mech. 95(1), 159–176 (1979). [CrossRef]  

29. M. D. White, “High-order parabolic beam approximation for aero-optics,” J. Comput. Phys. 229(15), 5465–5485 (2010). [CrossRef]  

30. S. Clifford, “The classical theory of wave propagation in a turbulent medium,” in Laser beam propagation in the atmosphere, (Springer, 1978), pp. 9–43.

31. S. Orszag and C. M. Bender, Advanced mathematical methods for scientists and engineers (McGraw-Hill New York, 1978).

32. J. K. Kevorkian and J. D. Cole, Multiple scale and singular perturbation methods, vol. 114 (Springer Science & Business Media, 2012).

33. M. Charnotskii, “Extended huygens–fresnel principle and optical waves propagation in turbulence: discussion,” J. Opt. Soc. Am. A 32(7), 1357–1365 (2015). [CrossRef]  

34. A. Nair, Q. Li, and S. Stechmann, “Matlab code for optimal beams in weak turbulence,” figshare(2022), https://doi.org/10.6084/m9.figshare.20439354.v1.

Data availability

Data and code underlying the results presented in this paper are available in [34].

34. A. Nair, Q. Li, and S. Stechmann, “Matlab code for optimal beams in weak turbulence,” figshare(2022), https://doi.org/10.6084/m9.figshare.20439354.v1.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (7)
Fig. 1.
Fig. 1. $\mathcal {H}$ computed using (72) when the receiver size $R$ is $0.01$, $0.03$, $0.06$, and $0.09$ m. The transmitter size is fixed at $r=0.05$ m.

View in Article | Download Full Size | PDF

Fig. 2.
Fig. 2. In the panel on the left, we show the comparison of total intensity over various receiver sizes, and the focused beam solution. Receiver size is in meters. In the panel on the right, we show the agreement of the first nine eigenvalues given by $\mathcal {H}$ and $\mathcal {H}_\mathrm {f}$, taking $R=0.09$m.

View in Article | Download Full Size | PDF

Fig. 3.
Fig. 3. Eigenfunctions provided by computing $\mathcal {H}$ (blue line) and $\mathcal {H}_{\mathrm {f}}$ (red line) are on top of each other. The left and the middle panel show the real and imaginary parts of the first six engenvectors, and the panel on the right shows the light intensity received at the receiver with the $R = 0.09$m for each of the six eigenfunctions.

View in Article | Download Full Size | PDF

Fig. 4.
Fig. 4. The plot in blue is the relative error in Fourier space vs physical space calculations for $\mathcal {H}$ in a uniform medium (with respect to the Frobenious norm). The red plot shows the relative difference in the presence and absence of turbulence from Fourier space calculations. $\epsilon$ is fixed at $5\times 10^{-8}$ ($R$ is in meters)

View in Article | Download Full Size | PDF

Fig. 5.
Fig. 5. In the panel on the left, we show the initial intensity profile of optimal beam using $R=0.05$m. In the panel on the right, we show the comparison of total intensity at the receiver using optimal beams and focused beams under the homogeneous random medium assumption and small-length-scale cutoff assumption ($R$ is in meters).

View in Article | Download Full Size | PDF

Fig. 6.
Fig. 6. Optimal intensity at the receiver as $\epsilon$ increases ($R$ is in meters)

View in Article | Download Full Size | PDF

Fig. 7.
Fig. 7. Initial intensity profile of optimal beams for different receiver sizes for $\epsilon =5\times 10^{-8}$ ($R$ is in meters)

View in Article | Download Full Size | PDF

Equations (87)

Equations on this page are rendered with MathJax. Learn more.

H = E [ H ] , with H ( X 1 , X 2 ) = ∫ X ′ ∈ R h ( X 1 , X ′ ) h ∗ ( X 2 , X ′ ) d X ′ .
U ( X ′ ) = ∫ X ∈ A h ω ( X , X ′ ) ϕ ( X ) d X .
I i ( X ′ ) = | U ( X ′ ) | 2 = ∫ X 1 , X 2 ∈ A h ( X 1 , X ′ ) ϕ ( X 1 ) h ∗ ( X 2 , X ′ ) ϕ ∗ ( X 2 ) d X 1 d X 2 , X ′ ∈ R ,
I ( X ′ ) = E [ ∫ X 1 ∈ A ∫ X 2 ∈ A h ( X 1 , X ′ ) J ( X 1 , X 2 ) h ∗ ( X 2 , X ′ ) d X 1 d X 2 ] , X ′ ∈ R ,
max J I R = max J ∫ X ′ ∈ R I ( X ′ ) d X ′ , subject to I 0 = 1 ,
J ( X 1 , X 2 ) = ∑ k α k ψ k ( X 1 ) ψ k ∗ ( X 2 ) , X 1 , X 2 ∈ A ,
I 0 = ∑ k α k = 1 ,
I R = ∑ k α k ∫ X 1 ∈ A ∫ X 2 ∈ A ψ k ( X 1 ) E [ H ( X 1 , X 2 ) ] ψ k ∗ ( X 2 ) d X 1 d X 2 ,
J ( X 1 , X 2 ) = α 1 ψ 1 ( X 1 ) ψ 1 ∗ ( X 2 ) , with α 1 = 1 ,
∇ ⊥ 2 A + 2 i k ∂ A ∂ z + k 2 ( n 2 n r 2 − 1 ) A = 0 ,
A ( X , z = 0 ) = ϕ ( X ) , and U ( X ′ ) = A ( X ′ , z = Z ) .
n 2 = n r 2 + ϵ n 1 2
A = A 0 + ϵ A 1 + ϵ 2 A 2 + ⋯ .
ϵ 0 : ∇ ⊥ 2 A 0 + 2 i k ∂ A 0 ∂ z = 0 ,
ϵ 1 : ∇ ⊥ 2 A 1 + 2 i k ∂ A 1 ∂ z = − k 2 n 1 2 n r 2 A 0 ,
ϵ 2 : ∇ ⊥ 2 A 2 + 2 i k ∂ A 2 ∂ z = − k 2 n 1 2 n r 2 A 1 .
∇ ⊥ 2 A 0 + 2 i k ∂ A 0 ∂ z = 0 , A ( X , 0 ) = ϕ ( X ) .
A ^ 0 ( K , z ) = ϕ ^ ( K ) G ^ ( K , z ) with G ^ ( K , z ) = exp ⁡ ( − i z 2 k ‖ K ‖ 2 ) .
G ( X , z ) = ( 1 2 π ) d / 2 ( k i z ) d / 2 exp ⁡ ( i k 2 z ‖ X ‖ 2 ) .
A 0 ( X , z ) = ϕ ∗ G ( X , z ) ,
A ^ 1 ( K , z ) = i k 2 ∫ z ′ = 0 z S ( K , z ′ ) G ^ ( K , z − z ′ ) d z ′ ,
A 1 ( X , z ) = i k 2 ∫ z ′ = 0 z [ ( V A 0 ) ( ⋅ , z ′ ) ∗ G ( ⋅ , z − z ′ ) ] ( X ) d z ′ ,
[ V A 0 ( ⋅ , z ′ ) ∗ G ( ⋅ , z − z ′ ) ] ( X ) = ∫ X ′ ∈ R d V A 0 ( X ′ , z ′ ) G ( X − X ′ , z − z ′ ) ) d X ′ .
A 2 ( X , z ) = i k 2 ∫ z ′ = 0 z d z ′ [ ( V A 1 ) ( ⋅ , z ′ ) ∗ G ( ⋅ , z − z ′ ) ] ( X ) .
h = h 0 + ϵ h 1 + ϵ 2 h 2 + high orders
U ( X ′ ) = A 0 ( X ′ , Z ) = ∫ X ∈ A ϕ ( X ) G ( X ′ − X , Z ) d X ,
h 0 ( X , X ′ ) = G ( X ′ − X , Z ) .
h 1 ( X , X ′ ) = i k 2 ∫ z 1 = 0 Z d z 1 ∫ X 1 ′ ∈ R d V ( X 1 ′ , z 1 ) G ( X ′ − X 1 ′ , Z − z 1 ) G ( X 1 ′ − X , z 1 ) d X 1 ′ .
h 2 ( X , X ′ ) = − k 2 4 ∫ z 1 = 0 Z ∫ X 1 ′ ∈ R d V ( X 1 ′ , z 1 ) G ( X ′ − X 1 ′ , Z − z 1 ) d X 1 ′ d z 1 ∫ z 2 = 0 z 1 ∫ X 2 ′ ∈ R d V ( X 2 ′ , z 2 ) G ( X 1 ′ − X 2 ′ , z 1 − z 2 ) G ( X 2 ′ − X , z 2 ) d X 2 ′ d z 2 .
U ^ ( K ) = A ^ ( K , Z ) = ∫ K ′ ∈ R d ϕ ^ ( K ′ ) g ( K ′ , K ) d K ′ .
g = g 0 + ϵ g 1 + ϵ 2 g 2 + high orders ,
g 0 ( K ′ , K ) = G ^ ( K , Z ) δ ( K − K ′ ) .
U ^ 1 ( K ) = i k 2 ∫ z 1 = 0 Z S ( K , z 1 ) G ^ ( K , Z − z 1 ) d z 1 = i k 2 ( 2 π ) d ∫ z 1 = 0 Z G ^ ( K , Z − z 1 ) d z 1 ∫ K ′ ∈ R d V ^ ( K − K ′ , z 1 ) G ^ ( K ′ , z 1 ) ϕ ^ ( K ′ ) d K ′ ,
g 1 ( K ′ , K ) = i k 2 ( 2 π ) d ∫ z 1 = 0 Z G ^ ( K , Z − z 1 ) G ^ ( K ′ , z 1 ) V ^ ( K − K ′ , z 1 ) d z 1 .
g 2 ( K ′ , K ) = − ( 1 2 π ) 2 d k 2 4 ∫ z 1 = 0 Z G ^ ( K , Z − z 1 ) d z 1 ∫ K ′ ′ ∈ R d V ^ ( K − K ′ ′ , z 1 ) d K ′ ′ × ∫ z 2 = 0 z 1 G ^ ( K ′ ′ , z 1 − z 2 ) G ^ ( K ′ , z 2 ) V ^ ( K ′ ′ − K ′ , z 2 ) d z 2 .
H = H 00 + ϵ 2 ( H 02 + H 11 + H 20 ) ,
H 00 ( X 1 , X 2 ) = ∫ X ′ ∈ R G ( X 1 − X ′ , Z ) G ∗ ( X 2 − X ′ , Z ) d X ′ = ( 1 2 π ) d ( k Z ) d ∫ X ′ ∈ R exp ⁡ ( i k 2 Z ( ‖ X 1 − X ′ ‖ 2 − ‖ X 2 − X ′ ‖ 2 ) ) d X ′ = G ( X 1 , Z ) G ∗ ( X 2 , Z ) ∫ X ′ ∈ R exp ⁡ ( − i k Z X ′ ⋅ ( X 1 − X 2 ) ) d X ′ .
E [ V ( X 1 , z 1 ) V ∗ ( X 2 , z 2 ) ] = Γ V ( X 1 , z 1 , X 2 , z 2 ) ,
H 11 ( X 1 , X 2 ) = k 2 4 ∫ z 1 , z 2 = 0 Z d z 1 d z 2 ∫ R d × R d Γ V ( X 1 ′ , z 1 , X 2 ′ , z 2 ) G ( X 1 − X 1 ′ , z 1 ) G ∗ ( X 2 − X 2 ′ , z 2 ) d X 1 ′ d X 2 ′ × ∫ X ′ ∈ R G ( X ′ − X 1 ′ , Z − z 1 ) G ∗ ( X ′ − X 2 ′ , Z − z 2 ) d X ′ .
H 02 ( X 1 , X 2 ) = − k 2 4 ∫ z 1 = 0 Z ∫ z 2 = 0 z 1 d z 2 d z 1 ∫ X 1 ′ ∈ R d ∫ X 2 ′ ∈ R d Γ V ( X 1 , ′ z 1 , X 2 ′ , z 2 ) G ∗ ( X 1 ′ − X 2 ′ , z 1 − z 2 ) G ∗ ( X 2 ′ − X 2 , z 2 ) d X 1 ′ d X 2 ′ × ∫ X ′ ∈ R G ( X ′ − X 1 , Z ) G ∗ ( X ′ − X 1 ′ , Z − z 1 ) d X ′ ,
I R = ∫ X ′ ∈ R E [ ⟨ U ( X ′ ) U ∗ ( X ′ ) ⟩ ] d X ′ = ∫ X ′ ∈ R d E [ ⟨ U ( X ′ ) w ( X ′ ) U ∗ ( X ′ ) w ∗ ( X ′ ) ⟩ ] d X ′ ,
I R = 1 ( 2 π ) d ∫ K ∈ R d E [ ⟨ W ( K ) W ∗ ( K ) ⟩ ] d K ,
W ( K ) = U w ^ ( K ) = 1 ( 2 π ) d U ^ ∗ w ^ .
I R = 1 ( 2 π ) 3 d ∫ K 1 , K 2 ∈ R d E [ ⟨ U ^ ( K 1 ) U ^ ∗ ( K 2 ) ⟩ ] ( ∫ K ∈ R d w ^ ( K − K 1 ) w ^ ∗ ( K − K 2 ) d K ) d K 1 d K 2 ,
1 ( 2 π ) d ∫ K ∈ R d w ^ ( K − K 1 ) w ^ ∗ ( K − K 2 ) d K = ∫ X ′ ∈ R ∫ X ′ ′ ∈ R d X ′ d X ′ ′ e i ( K 1 ⋅ X ′ − K 2 ⋅ X ′ ′ ) δ ( X ′ − X ′ ′ ) = ∫ X ′ ∈ R d X ′ e i ( K 1 − K 2 ) ⋅ X ′ = w ^ ( K 1 − K 2 )
I R = ∑ i j ϵ i + j ( 2 π ) 2 d ∫ K 1 , K 2 ∈ R d E [ ⟨ U ^ i ( K 1 ) U ^ j ∗ ( K 2 ) ⟩ ] w ^ ( K 1 − K 2 ) d K 1 d K 2 = ∑ i j ϵ i + j ( 2 π ) 2 d ∫ K 1 , K 2 ∈ R d ⟨ ϕ ^ ( K 1 ) ϕ ^ ∗ ( K 2 ) ⟩ d K 1 d K 2 × ∫ K 1 ′ , K 2 ′ ∈ R d E [ g i ( K 1 , K 1 ′ ) g j ∗ ( K 2 , K 2 ′ ) ] w ^ ( K 1 ′ − K 2 ′ ) d K 1 ′ d K 2 ′ .
J ^ ( K 1 , K 2 ) = ∫ X 1 ∈ R d ∫ X 2 ∈ R d J ( X 1 , X 2 ) e − i ( K 1 ⋅ X 1 − K 2 ⋅ X 2 ) d X 2 d X 1 ,
I R = ∑ i j ϵ i + j ( 2 π ) 2 d ∫ K 1 , K 2 ∈ R d J ^ ( K 1 , K 2 ) d K 1 d K 2 × ∫ K 1 ′ , K 2 ′ ∈ R d E [ g i ( K 1 , K 1 ′ ) g j ∗ ( K 2 , K 2 ′ ) ] w ^ ( K 1 ′ − K 2 ′ ) d K 1 ′ d K 2 ′ = 1 ( 2 π ) d ∫ K 1 , K 2 ∈ R d J ^ ( K 1 , K 2 ) H f ( K 1 , K 2 ) d K 1 d K 2 ,
H f , i j = 1 ( 2 π ) d ∫ K 1 ′ ∈ R d ∫ K 2 ′ ∈ R d E [ g i ( K 1 , K 1 ′ ) g j ∗ ( K 2 , K 2 ′ ) ] w ^ ( K 1 ′ − K 2 ′ ) d K 1 ′ d K 2 ′ ,
H f = H f , 00 + ϵ 2 ( H f , 02 + H f , 11 + H f , 20 ) .
J ^ ( K 1 , K 2 ) = ψ ( K 1 ) ψ ∗ ( K 2 ) ,
H f , 00 ( K 1 , K 2 ) = 1 ( 2 π ) d G ^ ( K 1 , Z ) G ^ ∗ ( K 2 , Z ) w ^ ( K 1 − K 2 ) .
H f , 11 ( K 1 , K 2 ) = ( 1 2 π ) 3 d k 2 4 ∫ z 1 = 0 Z ∫ z 2 = 0 Z G ^ ( K 1 , z 1 ) G ^ ∗ ( K 2 , z 2 ) d z 2 d z 1 × ∫ K 1 ′ , K 2 ′ ∈ R d G ^ ( K 1 ′ , Z − z 1 ) G ^ ∗ ( K 2 ′ , Z − z 2 ) Γ ^ V ( K 1 ′ − K 1 , z 1 , K 2 ′ − K 2 , z 2 ) w ^ ( K 1 ′ − K 2 ′ ) d K 1 ′ d K 2 ′ .
Γ ^ V ( K 1 , z 1 , K 2 , z 2 ) = E [ V ^ ( K 1 , z 1 ) V ^ ∗ ( K 2 , z 2 ) ] .
H f , 02 ( K 1 , K 2 ) = − ( 1 2 π ) 3 d k 2 8 G ^ ( K 1 , Z ) ∫ z 1 = 0 Z ∫ z 2 = 0 Z G ^ ∗ ( K 2 , z 2 ) d z 2 d z 1 × ∫ K 1 ′ , K 2 ′ ∈ R d G ^ ∗ ( K 2 ′ , z 1 − z 2 ) G ^ ∗ ( K 1 ′ , Z − z 1 ) Γ ^ V ( K 2 − K 2 ′ , z 2 , K 1 ′ − K 2 ′ , z 1 ) w ^ ( K 1 − K 1 ′ ) d K 1 ′ d K 2 ′ .
H f , 20 ( K 1 , K 2 ) = H f , 02 ∗ ( K 1 , K 2 ) .
I = 1 ( 2 π ) d ∫ K 1 ∈ R d ∫ K 2 ∈ R d J ^ ( K 1 , K 2 ) H f ( K 1 , K 2 ) d K 1 d K 2 = 1 ( 2 π ) d ∫ K 1 ∈ R d ∫ K 2 ∈ R d H f ( K 1 , K 2 ) d K 1 d K 2 ∫ X 1 ∈ A ∫ X 2 ∈ A J ( X 1 , X 2 ) e − i ( K 1 ⋅ X 1 − K 2 ⋅ X 2 ) d X 1 d X 2 = 1 ( 2 π ) d ∫ X 1 ∈ A ∫ X 2 ∈ A J ( X 1 , X 2 ) d X 1 d X 2 ∫ K 1 ∈ R d ∫ K 2 ∈ R d H f ( K 1 , K 2 ) e − i ( K 1 ⋅ X 1 − K 2 ⋅ X 2 ) d K 1 d K 2
H ( X 1 , X 2 ) = 1 ( 2 π ) d ∫ K 1 ∈ R d ∫ K 2 ∈ R d H f ( K 1 , K 2 ) e − i ( K 1 ⋅ X 1 − K 2 ⋅ X 2 ) d K 1 d K 2 ,
w ^ ( K 1 − K 2 ) = ∫ X ′ ∈ R d d X ′ e i ( K 1 − K 2 ) ⋅ X ′ = ( 2 π ) d δ ( K 1 − K 2 ) .
H f , 00 ( K 1 , K 2 ) = δ ( K 1 − K 2 ) , and H f , 11 + H f , 02 + H f , 20 = 0 ,
I = 1 ( 2 π ) d ∫ K 1 ∈ R d ∫ K 1 ∈ R d J ^ ( K 1 , K 2 ) δ ( K 1 − K 2 ) d K 1 d K 2 = 1 ( 2 π ) d ∫ K 1 ∈ R d J ^ ( K 1 , K 1 ) d K 1 = ∫ X ∈ A J ( X , X ) d X = I 0 .
Γ V ( X 1 , z 1 , X 2 , z 2 ) = f ( X 1 − X 2 , z 1 − z 2 ) .
Γ ^ V ( K 1 , z 1 , K 2 , z 2 ) = ( 2 π ) d F ( K 1 , z 1 − z 2 ) δ ( K 1 − K 2 ) ,
H f , 11 ( K 1 , K 2 ) = k 2 4 ( 2 π ) d H f , 00 ( K 1 , K 2 ) ∫ z 1 = 0 Z ∫ z 2 = 0 Z d z 2 d z 1 × ∫ K ′ ∈ R d F ( K ′ , z 1 − z 2 ) G ^ ( K ′ , z 1 − z 2 ) exp ⁡ ( − i k K ′ ⋅ ( K 1 z 1 − K 2 z 2 ) ) d K ′ = k 2 4 ( 2 π ) d H f , 00 ( K 1 , K 2 ) ∫ z 1 = 0 Z ∫ z 2 = 0 Z F G ^ ^ ( ( K 1 z 1 − K 2 z 2 ) / k , z 1 − z 2 ) d z 2 d z 1 .
H f , 02 ( K 1 , K 2 ) = − k 2 8 ( 2 π ) d H f , 00 ( K 1 , K 2 ) ∫ z 1 = 0 Z ∫ z 2 = 0 Z F G ^ ^ ( K 2 ( z 2 − z 1 ) / k , z 1 − z 2 ) d z 2 d z 1 .
F ( K , z ) ∼ 0 , for | z | > δ z .
η = z 1 + z 2 2 , μ = z 1 − z 2 .
H f , 11 ( K 1 , K 2 ) = k 2 4 ( 2 π ) d H f , 00 ( K 1 , K 2 ) ∫ K ′ ∈ R d d K ′ ∫ η = 0 Z / 2 exp ⁡ ( − i k K ′ ⋅ ( K 1 − K 2 ) η ) d η × ∫ μ = − 2 η 2 η F ( K ′ , μ ) G ^ ( K ′ , μ ) exp ⁡ ( − i 2 k K ′ ⋅ ( K 1 + K 2 ) μ ) d μ + k 2 4 ( 2 π ) d H f , 00 ( K 1 , K 2 ) ∫ K ′ ∈ R d d K ′ ∫ η = Z / 2 Z exp ⁡ ( − i k K ′ ⋅ ( K 1 − K 2 ) η ) d η × ∫ μ = − 2 ( Z − η ) 2 ( Z − η ) F ( K ′ , μ ) G ^ ( K ′ , μ ) exp ⁡ ( − i 2 k K ′ ⋅ ( K 1 + K 2 ) μ ) d μ ≈ k 2 4 ( 2 π ) d H f , 00 ( K 1 , K 2 ) ∫ K ′ ∈ R d d K ′ ∫ η = 0 Z exp ⁡ ( − i k K ′ ⋅ ( K 1 − K 2 ) η ) d η × ∫ μ ∈ R F ( K ′ , μ ) G ^ ( K ′ , μ ) exp ⁡ ( − i 2 k K ′ ⋅ ( K 1 + K 2 ) μ ) d μ
Φ ( K , K z ) = ϵ 2 ∫ z ∈ R ∫ X ∈ R d f ( X , z ) e − i ( K ⋅ X + K z z ) d X d z = ϵ 2 ∫ z ∈ R F ( K , z ) e − K z z d z .
‖ K ‖ 2 z 2 < ‖ K ‖ 2 < 1 2 l 0 .
G ^ = exp ⁡ { − i z 2 k ‖ K ‖ 2 } ∼ exp ⁡ ( − i 2 k l 0 ) ∼ exp ⁡ { 0 } = 1 ,
H f , 11 ( K 1 , K 2 ) ≈ Z k 2 4 ( 2 π ) d H f , 00 ( K 1 , K 2 ) ∫ K ′ ∈ R d exp ⁡ ( − i 2 k K ′ ⋅ ( K 1 − K 2 ) Z ) s i n c ( K ′ 2 k ⋅ ( K 1 − K 2 ) Z ) d K ′ × ∫ μ ∈ R F ( K ′ , μ ) exp ⁡ ( − i 2 k K ′ ⋅ ( K 1 + K 2 ) μ ) d μ = Z k 2 4 ( 2 π ) d ϵ 2 H f , 00 ( K 1 , K 2 ) ∫ K ′ ∈ R d exp ⁡ ( − i 2 k K ′ ⋅ ( K 1 − K 2 ) Z ) s i n c ( K ′ 2 k ⋅ ( K 1 − K 2 ) Z ) × Φ ( K ′ , K ′ ⋅ ( K 1 + K 2 ) / 2 k ) d K ′ .
H f , 02 ( K 1 , K 2 ) ≈ − Z k 2 8 ( 2 π ) d ϵ 2 H f , 00 ( K 1 , K 2 ) ∫ K ′ ∈ R d Φ ( K ′ , − K ′ ⋅ K 2 / k ) d K ′
f ( x , z ) = f 1 ( x ) f 2 ( z ) ,
H f , 11 = k 2 4 ( 2 π ) d H f , 00 ( K 1 , K 2 ) ∫ z 1 = 0 Z ∫ z 2 = 0 Z f 2 ( z 1 − z 2 ) F 1 G ^ ( ⋅ ^ , z 1 − z 2 ) ( ( K 1 z 1 − K 2 z 2 ) / k ) d z 2 d z 1 ,
H f , 02 ( K 1 , K 2 ) = − k 2 8 ( 2 π ) d H f , 00 ( K 1 , K 2 ) ∫ z 1 = 0 Z ∫ z 2 = 0 Z f 2 ( z 1 − z 2 ) F 1 G ^ ( ⋅ ^ , z 1 − z 2 ) ( K 2 ( z 2 − z 1 ) / k ) d z 2 d z 1 ,
f ( x , z ) = f 1 ( x ) δ ( z ) ,
H f , 11 = Z k 2 4 H f , 00 ( K 1 , K 2 ) ∫ y = 0 1 f 1 ( ( K 2 − K 1 ) y Z / k ) d y ,
H f , 02 ( K 1 , K 2 ) = − Z k 2 8 f 1 ( 0 ) H f , 00 ( K 1 , K 2 ) ,
H 11 = k 2 4 ∫ z 1 = 0 Z ( k 2 π ( Z − z 1 ) ) d G ( X 1 , z 1 ) G ∗ ( X 2 , z 1 ) d z 1 ∫ μ ∈ R d f ( μ ) exp ⁡ ( − i k μ ( X 1 + X 2 ) 2 z 1 ) d μ × ∫ X ′ ∈ R exp ⁡ ( − i k X ′ ⋅ μ Z − z 1 ) d X ′ ∫ η ∈ R d exp ⁡ ( i k η ( Z μ − ( X 1 − X 2 ) ( Z − z 1 ) z 1 ( Z − z 1 ) ) ) d η = Z k 2 4 H 00 ( X 1 , X 2 ) ∫ y = 0 1 f 1 ( ( X 1 − X 2 ) y ) d y .
H 02 = − Z k 2 8 f 1 ( 0 ) H 00 ( X 1 , X 2 ) .
H 00 ( X 1 , X 2 ) = ( R k π Z ) exp ⁡ ( i k 2 Z ( X 1 2 − X 2 2 ) ) s i n c ( R k Z ( X 1 − X 2 ) ) , ( X 1 , X 2 ) ∈ [ − r , r ] .
H f , 00 ( K 1 , K 2 ) = R π exp ⁡ ( − i Z 2 k ( K 1 2 − K 2 2 ) ) s i n c ( R ( K 1 − K 2 ) ) .
Γ ( X 1 , z 1 , X 2 , z 2 ) = exp ⁡ ( − ( X 1 − X 2 ) 2 + ( z 1 − z 2 ) 2 l 0 2 )
F ( K , z ) = π l 0 2 exp ⁡ ( − l 0 2 K 2 4 ) exp ⁡ ( − z 2 l 0 2 ) .
F G ^ ^ ( K ′ , z 1 − z 2 ) = π l 0 2 π s ( z 1 − z 2 ) exp ⁡ ( − ( z 1 − z 2 ) 2 l 0 2 ) exp ⁡ ( − K ′ 2 4 s ( z 1 − z 2 ) ) ,
Φ ( K , K z ) = ϵ 2 π l 0 2 exp ⁡ ( − l 0 2 ( K 2 + K z 2 ) 4 ) .
Select as filters
Select Topics Cancel
Top
       
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.
Privacy | Terms of Use