# Ermakov-Pinney equation for time-varying mass systems

Bijan Bagchi  Physics Department, Shiv Nadar University, Gautam Buddha Nagar,
Uttar Pradesh 201314, India  Achal Vinod  Physics Department, Shiv Nadar
University, Gautam Buddha Nagar,
Uttar Pradesh 201314, India

We extend Fring-Tenney approach of constructing invariants of constant mass
time-dependent system to the case of a time-dependent mass particle. From a
coupled set of equations described in terms of guiding parameter functions, we
track down a modified Ermakov-Pinney equation involving a time-dependent mass
function. As a concrete example we focus on an exponential choice of the mass
function.

††E-mails<EMAIL_ADDRESS><EMAIL_ADDRESS>
Keywords: Time-dependent Schrödinger equation, Lewis-Riesenfeld method,
Ermakov-Pinney equation

##  1 Introduction

For some time now, there have been several efforts to study time-dependent
quantum systems from different perspectives. One such consists in the use of
optical traps (see for some representative articles [1, 2, 3, 4, 5]). In the
non-relativistic context, the governing equation of concern is the time-
dependent Schrödinger equation with the explicit presence of an optical
potential. It is usually taken to be in a separable form which is made up of a
time-dependent modulation term together with some trapping shape described by
a spatial profile [6, 7].

Interest in exploring the dynamics of time-dependent quantum systems started
with the method of seeking Hermitian invariants [8, 9, 10] for them. Recently,
Fring and Tenney [11, 12] made a systematic study of time-independent
approximations for a class of time-dependent optical potentials. Their general
strategy was to write down an approximate Lewis-Riesenfeld scheme [9] and
explore, among other things, time-dependent potentials with a Stark term.
While constructing time-independent invariants, even for non-Hermitian systems
[12, 13, 14], one of the central issues was to derive the highly nonlinear
Ermakov-Pinney $(EP)$ equation [15, 16]. Interestingly, they could obtain a
regular solution to it, determine the electric field and finally construct the
invariants. Note that they concentrated only on a constant mass system. For
other approaches towards calculating dynamical invariants we refer to [17].

Historically, derivation of the $EP$ equation (see [18] in which some symmetry
aspects of this equation were explored) was first undertaken in an early work
of Ermakov [15] while studying integrablity of certain classes of ordinary
second-order differential equation

$\ddot{x}+\omega^{2}x=\frac{\lambda^{2}}{x^{3}},\quad\lambda>0$ (1.1)

We see that the $EP$ equation describes an oscillator in the presence of an
inverse-cube force term.

Much later, Pinney [16] provided the solution and gave a complete
representation of it in the form

$x(t)=\left(Ar^{2}+2Brs+Cs^{2}\right)^{\frac{1}{2}}$ (1.2)

where $r,s$ are both functions of $t$ and stand as the independent solutions
of the standard harmonic oscillator $\ddot{x}+\omega^{2}x=0$ while $A,B,C$
satisfy the constraint $AC-B^{2}=\frac{h^{2}}{W}$, $W$ being the Wronskian.

The purpose of this article is to pursue an analog procedure as that of [11]
but consider a time-dependent mass instead of a constant one. The variation of
the mass with time is a well studied topic in problems of unstable particles
[19, 20] in the realm of particle physics. A modification to the form of the
Dirac equation in seeking invariants has also been considered for time-
dependent masses in Dirac motivated relativistic systems [21]. Cases of both
time-dependent and position-dependent mass have also been explored in the
literature (see, for instance, [22, 23]).

In the following, we basically utilize invariants to solve the time-dependent
Schrödinger equation when the guiding time-dependent potential is
factorizable111For a Hamiltonian which is time-dependent such a situation does
not in general occur in terms of a Gaussian profile apart from a time-
dependent piece. We see that a modified EP-equation emerges as an auxiliary
equation whose form differs from what is obtained in the constant mass case.
However, its highly nonlinear nature precludes us from finding any exact
analytical solution. We have therefore gone for a specific case-study of the
mass profile as given by an exponential representation. It reduces the basic
equations to a tractable form thus facilitating a simple numerical assessment.

This work is organized as follows:

In section 2, we highlight the consistency conditions that emerge when the
basic equation for an invariant Hermitian operator $\hat{I}$ is subjected to
the Hamiltonian obeying the time-dependent Schrödinger equation. In section 3
we derive a modified $EP$ equation by installing time-dependence in the mass.
Such an equation constitutes our subsequent point of analysis which is taken
up numerically. In Section 4 we present a particular class of solution
corresponding to the exponential choice of the mass. Finally, in section 5, we
present a summary.

##  2 Constraining equations for time-dependent coefficients

Let us begin with the time-dependent non-relativistic Schrödinger equation
which in one-dimension reads

$i\frac{\partial}{\partial
t}\psi(x,t)=\left[-\frac{1}{2m}\frac{\partial^{2}}{\partial
x^{2}}+V(x,t)\right]\psi(x,t)$ (2.1)

where $m=m(t)$. We work with natural units $\hbar=c=1$.

For a general quantum mechanical Hamiltonian $\hat{H}$, the time evolution
operator $U(t,t_{0})$ obeys [24]

$i\dot{U}(t,t_{0})=\hat{H}U(t,t_{0}),\quad U(t_{0},t_{0})=I,\quad t>t_{0}$
(2.2)

If the Hamiltonian does not depend on time then the solution of the above
equation is simply

$U(t,t_{0})=e^{-i\hat{H}(t-t_{0})}$ (2.3)

where Hamiltonian is considered to be Hermitian.

Against the background of the evolution equation obeying $(2.2)$, we note that
the time-dependent invariant operator $\hat{I}(t)$ satisfies

$\frac{d\hat{I}}{dt}=\dot{\hat{I}}(t)+\frac{1}{i}[\hat{I},H]=0$ (2.4)

where $H$ is time-dependent $\hat{H}=\hat{H}(x,t)$. As a concrete model of our
analysis we consider the following representative Hamiltonian

$\hat{H}(x,t)=\frac{p^{2}}{2m(t)}+V(x,t)+xE(t)$ (2.5)

where the mass and the potential have been taken to be time-dependent. The
latter is assumed to be factorizable222For a wider class of time-dependent
potentials that are factorizable and also exactly solvable see [25, 26] in
that $V(x,t)=\kappa(t)V(x)$, where $\kappa$ is a function of $t$.
Additionally, a Stark like term is included in $(2.5)$ which points to the
presence of an electric field $E(t)$. The angular frequency $\omega$ is taken
to be constant.

The choice of the representation for the invariant is flexible. In the
following we adopt the form

$\hat{I}(t)=\frac{1}{2}\left[\alpha(t)p^{2}+\gamma(t)x+\delta(t)\\{x,p\\}+\varepsilon(t)x^{2}\right]$
(2.6)

where $\\{x,p\\}=xp+px$ and the coefficient functions $\alpha,\gamma,\delta$
and $\varepsilon$ are time-dependent functions. Of course, these are distinct
from $m(t)$ or $\kappa(t)$ and have to be solved for through the consistency
equations.

We now address the simplest class of time-dependent Hamiltonian as induced by
a varying mass harmonic oscillator potential. In such a case $\kappa$
corresponds to $m(t)$ while $V(x)=\frac{1}{2}\omega^{2}x^{2}$. Substituting in
$(2.4)$ we obtain through the use of $(2.5)$ the commutator

$[\hat{I},\hat{H}]=\frac{1}{2}\bigg{(}\alpha
m(-\hbar^{2}\omega^{2}-2i\hbar\omega^{2}xp)+\alpha
E(-2i\hbar)p+\frac{\gamma}{2m}(2i\hbar)p+\frac{\delta}{2m}(4i\hbar)p^{2}\\\
+\delta m(-2i\hbar)\omega^{2}x^{2}+\delta
E(-2i\hbar)x+\frac{\varepsilon}{2m}(4i\hbar xp+2\hbar^{2})\bigg{)}$ (2.7)

When substituted in $(2.4)$ the following set of conditions is easily obtained

$\displaystyle\dot{\alpha}=-\frac{2\delta}{m}$ (2.8)
$\displaystyle\dot{\delta}=\alpha m\omega^{2}-\frac{\varepsilon}{m}$ (2.9)
$\displaystyle\dot{\varepsilon}=2\delta m\omega^{2}$ (2.10)
$\displaystyle\dot{\gamma}=2\delta E$ (2.11)

We also have an additional relation

$\displaystyle\gamma=2m\alpha E$ (2.12)

along with a closed form for $E(t)$ given by

$E(t)=\frac{1}{m}e^{\int Q(t)dt},\quad Q(t)=\frac{3\delta}{\alpha m}$ (2.13)

We should mention that in the constant mass case, the set of equations
$(2.8)-(2.12)$ reduce to those of the corresponding ones in [11].

##  3 Ermakov-Pinney equation and its modification

We aim at writing down a modified EP equation in the presence of $m(t)$.
Redefining $\alpha=\sigma^{2}$ the remaining parameters read

$\displaystyle\delta=-m\dot{\sigma}\sigma$ (3.1)
$\displaystyle\gamma=2\sigma^{2}mE(t)$ (3.2)
$\displaystyle\varepsilon=m^{2}\dot{\sigma}^{2}+m^{2}\sigma\ddot{\sigma}+m\dot{m}\dot{\sigma}\sigma+m^{2}\sigma^{2}\omega^{2}$
(3.3)

Introducing $\tau\equiv\tau(\sigma)$ as an arbitrary function of $\sigma$
which is defined through

$\displaystyle
m^{2}\sigma\ddot{\sigma}+m\dot{m}\dot{\sigma}\sigma+m^{2}\sigma^{2}\omega^{2}=\frac{\tau
m^{2}}{\sigma^{2}}$ (3.4)

it is clear that $(3.4)$ stands as the signature of a dissipative EP equation
for the time-varying mass case

$\ddot{\sigma}+\omega^{2}\sigma+\frac{\dot{m}}{m}\dot{\sigma}=\frac{\tau}{\sigma^{3}}$
(3.5)

We immediately see that when $m$ is stationary and $\tau$ is a constant equal
to $\lambda^{2}$ the standard equation $(1.1)$ is recovered. $(3.3)$ now reads

$\varepsilon=m^{2}\dot{\sigma}^{2}+\frac{\tau m^{2}}{\sigma^{2}}$ (3.6)

Clearly, no unique set of solutions of the parameter functions exists
satisfying $\sigma$-dependent equations $(3.1)$, $(3.2)$ and $(3.6)$, for a
given functional choice of the mass function $m(t)$ and the electric field
$E(t)$.

Identifying the EP equation in the extended form $(3.5)$ constitutes the
central result of this note. Although not exactly coincident with $(3.5)$, a
dissipative $EP$ equation was also derived by Fring and Tenney [12] but in a
different context

$\ddot{\chi}+f^{2}\chi-\frac{\dot{f}}{f}\dot{\chi}=\frac{f^{2}}{\chi^{3}},\quad\chi=\chi(t)$
(3.7)

where $f$ is a time-dependent parametrizing function. $f$ seems to play the
analog of the mass function but notice that while the coefficient of $\sigma$
in $(3.5)$ is a constant, the same in the corresponding $(3.7)$ is a function
of time.

Our next endeavour would be to analyse $(3.6)$ by taking its time-derivative
and comparing with $(2.10)$. Towards this end we obtain the pair of equations

$\displaystyle\dot{\varepsilon}=2m\dot{m}\dot{\sigma}^{2}+2\dot{\sigma}\ddot{\sigma}m^{2}+\frac{1}{\sigma^{2}}(2m\dot{m}\tau+m^{2}\tau^{\prime}\dot{\sigma})-\frac{2m^{2}\tau\dot{\sigma}}{\sigma^{3}}$
(3.8) $\displaystyle\dot{\varepsilon}=-2m^{2}\omega^{2}\sigma\dot{\sigma}$
(3.9)

where $(3.1)$ is used. Equating them we straightforwardly obtain a
differential equation for $\tau$

$\tau^{\prime}+2\frac{\dot{m}}{m\dot{\sigma}}\tau=0$ (3.10)

where the prime denotes a derivative with respect to $\sigma$.

##  4 Solving the modified $EP$ equation

Since $\tau=\tau(\sigma)$ let us assume a monomial form for it

$\tau=\tau_{0}\sigma^{4}$ (4.1)

where $\tau_{0}$ is a constant assumed to be $>0$. The EP equation, in
principle, admits of complex solutions given by the general form

$\sigma(t)=\xi(t)e^{i\eta(t)},\quad\xi,\eta\in\Re$ (4.2)

where $\xi$ and $\eta$ are linked by a set of coupled equations

$\displaystyle\ddot{\xi}+\frac{\dot{m}}{m}\dot{\xi}+(\Omega^{2}-\dot{\eta}^{2})=0$
(4.3) $\displaystyle
2\dot{\xi}\dot{\eta}+\xi\left(\ddot{\eta}+\frac{\dot{m}}{m}\dot{\eta}\right)=0$
(4.4)

where $\Omega^{2}=(\omega^{2}-\tau_{0})$. These are very complicated equations
to tackle. So we concentrate on real solutions only and in this regard
consider a mass function $m(t)$ of the exponential type

$m(t)=m_{0}e^{-qt},\quad q>0$ (4.5)

where its damping character is in keeping with [21, 22] but with respect to
the t-variable. For $\sigma$ we find the solution

$\sigma(t)=e^{qt/2}$ (4.6)

subject to the the following constraint on $\tau_{0}$

$\tau_{0}=\omega^{2}-\frac{q^{2}}{4}$ (4.7)

$\tau$ is insensitive to the sign of $q$. For completeness we note that
$\delta$ turns out to be a constant while the coefficient function $\gamma$
and the electric field $E$ show a decaying behaviour $\backsim e^{-qt/2}$.

We display the behavior of the mass variation in Figure 1 corresponding to
different values of $\omega$. We took the input value of $\tau_{0}=0.01$. The
results reflect a typical exponential trajectory of $m$.

Figure 1: $m(t)$ variation for $\omega=0.1,0.3,0.5,0.7,0.9$

##  5 Summary

In this article we considered the possibility of introducing a time-dependent
mass while constructing an invariant for a time-dependent quantum Hamiltonian.
Assuming a time-dependent harmonic oscillator with a Stark-like Hamiltonian we
were able to derive the invariant conditions in their entirety. We then
obtained the time-dependent mass version of the $EP$ equation that included an
additional term as compared to the constant-mass case. Because of its strongly
nonlinear character we inquired into simple parameter variations to assess the
mass function. For the specific case of a damping mass with a certain input
value of $\tau_{0}$, its graphical behavior was plotted over a wide range of
$t$.

##  6 Acknowledgment

We are indebted to Andreas Fring and Rebecca Tenney for pointing out errors in
the previous version of the draft.

## References

  * [1] Ashkin A 1970 Acceleration and trapping of particles by radiation pressure Phys. Rev. Lett. 24 156

  * [2] Neuman K C and Block S M 2004 Optical trapping Review of scientific instruments 75 2787

  * [3] Ashkin A, Dziedzic J M, Bjorkholm J E and Chu S 1986 Observation of a single-beam gradient force optical trap for dielectric particles Optics Letters 11 288

  * [4] Metcalf H J and Van der Straten P 2007 Laser cooling and trapping of neutral atoms, The Optics Encyclopedia: Basic Foundations and Practical Applications.

  * [5] Ashkin A, Dziedzic J M and Yamane T 1987 Optical trapping and manipulation of single cells using infrared laser beams Nature 330 769\.

  * [6] Ohara K M, Gehm M E, Granade S R and Thomas J E 2001 Scaling laws for evaporative cooling in time-dependent optical traps Phys. Rev. A64 051403

  * [7] Fulton R, Bishop A I, Shneider M N and Barker P F 2006 Controlling the motion of cold molecules with deep periodic optical potentials Nature Physics 2 465

  * [8] Lewis H R, Jr 1968 Class of exact invariants for classical and quantum time‐dependent harmonic oscillators J. Math. Phys. 9 1976

  * [9] Lewis H R, Jr and Riesenfeld W 1969 An exact quantum theory of the time dependent harmonic oscillator and of a charged particle time dependent electromagnetic field J. Math. Phys. 10 1458

  * [10] Lai Y -Z, Liang J -Q, Müller-Kirsten H J W and Zhou J -G 1996 Time-dependent quantum systems and the invariant Hermitian operator Phys. Rev. A53 3691

  * [11] Fring A and Tenney R 2020 Time-independent approximations for time-dependent optical potentials Eur. Phys. J. Plus 135 163

  * [12] Fring A and Tenney R 2021 Perturbative approach for strong and weakly coupled time-dependent non-Hermitian quantum systems Phys Scrip 96 045211

  * [13] Fring A and Frith T 2017 Mending the broken PT-regime via an explicit time-dependent Dyson map Phys. Lett. A381 2318

  * [14] Koussa W and Maamache M 2020 Pseudo-Invariant Approach for a Particle in a Complex Time-Dependent Linear Potential Int. J. Theor. Phys. 59 1490

  * [15] Ermakov V 1880 Second order differential equations. Conditions of complete integrability Universita Izvestia Kiev Series III 9 1

  * [16] Pinney E 1950 The nonlinear differential equation $y^{{}^{\prime\prime}}+p(x)y+\frac{c}{y^{3}}=0$ Proc. Amer. Math. Soc. 1 681

  * [17] Bertin M C, Pimentek B M and Ramirez J A 2012 Construction of time-dependent dynamical invariants: A new approach J. Math. Phys. 53 042104

  * [18] Morris R M and Leach P G L 2015 The Ermakov-Pinney Equation: its varied origins and the effects of the introduction of symmetry-breaking functions arXiv:1510.08992v1 (math.CA)

  * [19] Bekkar H, Benamira F and Maamache M 2003 Comment on “Solution of the Schrödinger equation for the time-dependent linear potential” Phys. Rev. A 68 016101

  * [20] Landim R R and Guedes I 2000 Wave functions for a Dirac particle in a time-dependent potential Phys. Rev. A61 054101

  * [21] Zhang Z-G 2007 Complete solutions for a class of time-dependent Dirac equations Phys. Scr. 76 349

  * [22] Pedrosa I A 1997 Exact wave functions of a harmonic oscillator with time-dependent mass and frequency Phys. Rev. A 55 3219

  * [23] Bagchi B, Gorain P, Quesne C and Roychoudhury R 2005 New approach to (quasi)-exactly solvable Schrödinger equations with a position-dependent effective mass Europhys. Lett. 72 155

  * [24] Sakurai J J 1994 Modern quantum mechanics Addison-Wesley Publishing Company, Inc.

  * [25] Rogers J and Spector D 1992 Quantum mechanics with explicit time-dependence Phys. Lett. A170 344

  * [26] Efthimion C J and Spector D 1994 Separation of variables and exactly soluble time-dependent potentials in quantum mechanics Phys. Rev. A49 2301


>     m(rf) pipe(sparse) opt(max_features(5))                              || ///
    >     m(gradboost) pipe(sparse) opt(n_estimators(250) learning_rate(0.01)) ,  ///
    >     njobs(4)
    Learner Z1_pystacked added successfully.

#### Step 3: Perform cross-fitting.

  *

        . ddml crossfit
    Cross-fitting E[y|X] equation: net_tfa
    Cross-fitting fold 1 2 3 4 5 ...completed cross-fitting
    Cross-fitting E[D|X] equation: p401
    Cross-fitting fold 1 2 3 4 5 ...completed cross-fitting
    Cross-fitting E[Z|X]: e401
    Cross-fitting fold 1 2 3 4 5 ...completed cross-fitting

    . ddml extract, show(pystacked)

    mean pystacked weights across folds/resamples for Z1_pystacked (e401)
                   learner  mean_weight
          ols            1    .02435374
      lassocv            2    5.963e-19
      ridgecv            3    .48160664
           rf            4    .03685229
    gradboost            5    .45718733

    mean pystacked MSEs across folds/resamples for Z1_pystacked (e401)
                 learner   mean_MSE
          ols          1  .20042018
      lassocv          2  .19661349
      ridgecv          3  .19661916
           rf          4  .21495839
    gradboost          5  .19662694

    mean pystacked weights across folds/resamples for Y1_pystacked (net_tfa)
                   learner  mean_weight
          ols            1     .0948191
      lassocv            2    .13620045
      ridgecv            3    .68323138
           rf            4    .07307122
    gradboost            5    .00888029

    mean pystacked MSEs across folds/resamples for Y1_pystacked (net_tfa)
                 learner   mean_MSE
          ols          1  3.256e+09
      lassocv          2  2.991e+09
      ridgecv          3  2.988e+09
           rf          4  3.418e+09
    gradboost          5  3.280e+09

    mean pystacked weights across folds/resamples for D1_pystacked (p401)
                   learner  mean_weight
          ols            1    .06282466
      lassocv            2    .13679104
      ridgecv            3     .4428182
           rf            4    .06142842
    gradboost            5    .29613768

    mean pystacked MSEs across folds/resamples for D1_pystacked (p401)
                 learner   mean_MSE
          ols          1  .17196874
      lassocv          2  .17048203
      ridgecv          3  .17061646
           rf          4  .18560235
    gradboost          5  .17094582

#### Step 4: Estimate causal effects.

  *

        . ddml estimate

    DDML estimation results:
    spec  r     Y learner     D learner         b        SE     Z learner
     opt  1  Y1_pystacked  D1_pystacked 13528.643(1726.023)
    opt = minimum MSE specification for that resample.

    Min MSE DDML model
    y-E[y|X]  = Y1_pystacked_1                         Number of obs   =      9915
    D-E[D|X,Z]= D1_pystacked_1
    Z-E[Z|X]  = Z1_pystacked_1

         net_tfa  Coefficient  Std. err.      z    P>|z|     [95% conf. interval]

            p401    13528.64   1726.023     7.84   0.000      10145.7    16911.59
           _cons   -42.68924    531.132    -0.08   0.936    -1083.689    998.3104

#### One-line syntax.

  *

        . qui qddml $Y ($X) ($D = $Z), model(iv) cmd(pystacked)

###  J.3 Interactive IV Model interactiveiv

  *

        . use "sipp1991.dta", clear

    . global Y net_tfa

    . global X age inc educ fsize marr twoearn db pira hown

    . global Z e401

    . global D p401

##### Step 1: Initialize ddml model.

  *

        . set seed 123

    . ddml init interactiveiv

##### Step 2: Add supervised machine learners for estimating conditional
expectations.

  *

        . * add learners for E[Y|X,Z=0] and E[Y|X,Z=0]
    . ddml E[Y|X,Z]: pystacked $Y $X                                           || ///
    >     method(ols)                                                          || ///
    >     m(lassocv) xvars(c.($X)##c.($X))                                     || ///
    >     m(ridgecv) xvars(c.($X)##c.($X))                                     || ///
    >     m(rf) pipe(sparse) opt(max_features(5))                              || ///
    >     m(gradboost) pipe(sparse) opt(n_estimators(250) learning_rate(0.01)) ,  ///
    >     njobs(4)
    Learner Y1_pystacked added successfully.

    . * add learners for E[D|X,Z=0] and E[D|X,Z=0]
    . ddml E[D|X,Z]: pystacked $D $X                                           || ///
    >     method(ols)                                                          || ///
    >     m(lassocv) xvars(c.($X)##c.($X))                                     || ///
    >     m(ridgecv) xvars(c.($X)##c.($X))                                     || ///
    >     m(rf) pipe(sparse) opt(max_features(5))                              || ///
    >     m(gradboost) pipe(sparse) opt(n_estimators(250) learning_rate(0.01)) ,  ///
    >     njobs(4)
    Learner D1_pystacked added successfully.

    . * add learners for E[Z|X,]
    . ddml E[Z|X]: pystacked $Z $X                                             || ///
    >     method(ols)                                                          || ///
    >     m(lassocv) xvars(c.($X)##c.($X))                                     || ///
    >     m(ridgecv) xvars(c.($X)##c.($X))                                     || ///
    >     m(rf) pipe(sparse) opt(max_features(5))                              || ///
    >     m(gradboost) pipe(sparse) opt(n_estimators(250) learning_rate(0.01)) ,  ///
    >     njobs(4)
    Learner Z1_pystacked added successfully.

##### Step 3: Perform cross-fitting.

  *

        . ddml crossfit
    Cross-fitting E[y|X,Z] equation: net_tfa
    Cross-fitting fold 1 2 3 4 5 ...completed cross-fitting
    Cross-fitting E[D|X,Z] equation: p401
    Cross-fitting fold 1 2 3 4 5 ...completed cross-fitting
    Cross-fitting E[Z|X]: e401
    Cross-fitting fold 1 2 3 4 5 ...completed cross-fitting

##### Step 4: Estimate causal effects.

  *

        . ddml estimate

    DDML estimation results (LATE):
    spec  r    Y0 learner    Y1 learner    D0 learner    D1 learner         b        SE     Z learner
     opt  1  Y1_pystacked  Y1_pystacked  D1_pystacked  D1_pystacked 11568.523(1613.594)  Z1_pystacked
    opt = minimum MSE specification for that resample.

    E[y|X,D=0]   = Y1_pystacked0_1                     Number of obs   =      9915
    E[y|X,D=1]   = Y1_pystacked1_1
    E[D|X,Z=0]   = D1_pystacked0_1
    E[D|X,Z=1]   = D1_pystacked1_1
    E[Z|X]       = Z1_pystacked_1

                                Robust
         net_tfa  Coefficient  std. err.      z    P>|z|     [95% conf. interval]

            p401    11568.52   1613.594     7.17   0.000     8405.938    14731.11

    Warning: 13 propensity scores trimmed to lower limit .01.

##### One-line syntax.

  *

        . qui qddml $Y ($X) ($D=$Z), model(late) cmd(pystacked)


Date of publication xxxx 00, 0000, date of current version xxxx 00, 0000.
10.1109/ACCESS.2017.DOI

This work was supported in part by the Optical Satellite Communications
Consortium Canada (OSC).

Corresponding author: Eylem Erdogan (e-mail: eylem.erdogan@medeniyet.edu.tr).

# Site Diversity in Downlink Optical Satellite Networks Through Ground Station
Selection

EYLEM ERDOGAN1  IBRAHIM ALTUNBAS2  GUNES KARABULUT KURT2,4  MICHEL BELLEMARE3
GUILLAUME LAMONTAGNE3  HALIM YANIKOMEROGLU4  Department of Electrical and
Electronics Engineering, Istanbul Medeniyet University, 34700, Uskudar,
Istanbul, Turkey  Electronics and Communication Engineering Department,
Istanbul Technical University, 34469, Maslak, Istanbul, Turkey  MDA, Satellite
Systems, Sainte-Anne-de-Bellevue, H9X 3R2, QC, Canada  Department of Systems
and Computer Engineering, Carleton University, Ottawa, K1S 5B6, ON, Canada

###### Abstract

Recent advances have shown that satellite communication (SatCom) will be an
important enabler for next generation terrestrial networks as it can provide
numerous advantages, including global coverage, high speed connectivity,
reliability, and instant deployment. An ideal alternative for radio frequency
(RF) satellites is its free-space optical (FSO) counterpart. FSO or laser
SatCom can mitigate the problems occurring in RF SatCom, while providing
important advantages, including reduced mass, lower consumption, better
throughput, and lower costs. Furthermore, laser SatCom is inherently resistant
to jamming, interception, and interference. Owing to these benefits, this
paper focuses on downlink laser SatCom, where the best ground station (GS) is
selected among numerous candidates to provide reliable connectivity and site
diversity. To quantify the performance of the proposed scheme, we derive
closed-form outage probability and ergodic capacity expressions for two
different practical GS deployment scenarios. Thereafter, asymptotic analysis
is conducted to obtain the overall site diversity order, and aperture
averaging is studied to illustrate the impact of aperture diameter on the
overall performance. Furthermore, we investigate the site diversity order for
a constellation of satellites that are communicating with the best GS by using
opportunistic scheduling. Finally, important design guidelines that can be
useful in the design of practical laser SatComs are outlined.

###### Index Terms:

Laser satellite communication, site diversity, free-space optical
communication, atmospheric turbulence and attenuation.

=-15pt

##  I Introduction

Satellite communication (SatCom) has become an important part of aerial
networks in recent years due to its capabilities, which include flawless
wireless connectivity, wide service coverage, and high-fidelity services for
all the users around the world. An important feature of SatCom is to
simultaneously transfer the signal rapidly around the Earth by providing
distance-insensitive point-to-multipoint communications [1]. So far,
satellites have been used for television coverage, data communication,
navigation, weather forecasts, climate and environmental monitoring, space
science, and so on [2]. Satellites can be divided into three main categories
depending on their altitudes and orbit types, low Earth orbit (LEO), medium
Earth orbit (MEO) and geostationary Earth orbit (GEO). GEO satellites are
located at about $36000$ km above the Earth’s equator to provide wide
coverage, whereas LEO satellites, which circle around the Earth at lower
altitudes, have faster rotations, require less power, and are cheaper than GEO
and MEO satellites [3]. The technology behind SatCom is mainly based on radio-
frequency (RF) systems, where $100$ MHz to $50$ GHz frequencies are used
depending on the types and applications of satellites. In recent years, SatCom
has emerged to provide high-speed, seamless broadband Internet connectivity
around the globe as many different companies have started to launch
constellations of satellites. For instance, SpaceX’s Starlink recently began
launching LEO satellites operating at high frequencies [4]. Similarly, OneWeb
has launched about $104$ LEO satellites in $2020$, which operate in the
$K_{a}$-$K_{u}$ band of the frequency spectrum [5]. These recent developments
suggest that RF-SatCom will become a key enabler in the integration of aerial
and terrestrial networks in future wireless communication systems. In RF-
SatCom, the most important drawbacks are cost, regulatory restrictions and
limited available bandwidth, as it requires high data rates and broader
bandwidths to connect anyone at anytime. Furthermore, RF-SatCom is prone to
interference, jamming, and interception, which pose security risks, especially
for military communications.

As a solution to these problems, free-space optical (FSO) or laser SatCom has
attracted considerable interest both in recent academic and industry
publications [6]. In laser SatCom, $20$ to $375$ THz spectrum111In laser
SatCom, only a small portion of the frequency spectrum can be used due to huge
atmospheric losses. can be used to provide very high throughput in satellite-
to-ground (downlink) and ground-to-satellite (uplink) communications [7].
Laser SatCom can also provide significant advantages compared to its RF
counterpart, including smaller antennas, reduced mass, lower consumption,
better throughput, and lower costs. Furthermore, the narrow beam used in
optical systems can provide secure communication, and it is immune to jamming,
interception, interference. In addition, laser SatCom does not require a
spectrum operating license for frequency use due to its inherent nature [8].
Due to these advantages, laser SatCom is expected to become a key enabler for
future optical satellite systems, particularly for satellite-to-ground
(downlink) communications, where line-of-sight (LOS) connectivity can be
established perfectly [9]. In downlink laser SatCom, the major adverse effects
are atmospheric turbulence, atmospheric attenuation, and angle of arrival
(AoA) fluctuations [7]. The latter can be attenuated by using variable-focus
lenses, which can adjust the beam size [10]. Also, aperture averaging, where
the scintillation is spatially averaged over the aperture, can be used to
reduce the adverse effects of atmospheric turbulence [11]. Finally,
atmospheric attenuation due to adverse weather conditions can be resolved by
achieving site diversity, in which the number of optical links can be enhanced
by using multiple ground stations [12, 13, 14, 15]. It can also be resolved by
attaining spatial diversity with the aid of multiple apertures on a single GS,
using appropriate combining techniques, such as selection combining, maximum
ratio combining or equal gain combining [16, 17].

Due to the advantages outlined above, optical satellites have become an
important topic in the recent literature. In [18, 19, 20, 21], the issue of
optical ground-to-satellite (uplink) communication was considered, and
important performance metrics, including outage probability, error
probability, and link capacity, were obtained in the presence beam wandering,
climatic effects, and atmospheric attenuation. By contrast, [16] and [22]
focused on downlink SatCom, where the former investigated the impact of
spatial diversity and aperture averaging, and the latter examined the
throughput. Furthermore, [23] and [24] proposed using FSO communication as a
feeder link where the ground station (GS) feeds the satellite through a high
capacity link.

In this paper, we pursue a different line of inquiry by providing a detailed
analysis for the downlink optical SatCom systems. More precisely, to deal with
the adverse weather conditions, we try to maximize the site diversity by
selecting the best GS that provides the best channel conditions. It is
important to note that site diversity, which can be crucial for creating
seamless connectivity between satellite-to-ground links, has already been
discussed in the literature e.g., [12, 13, 14, 15] and the references therein.
However, these previous works focused on the optimum GS selection problem from
the network layer point-of-view. To the best of our knowledge, none of these
works have considered the physical layer performance of the optical SatCom
with GS selection. Given the importance of site diversity, and to fill the gap
in the literature, we focus on the physical layer performance of the downlink
optical SatCom with multiple GSs222In the proposed setup, we assume that all
GSs are connected through fibre-optic wires, and the GS with the highest
communication reliability can share the information to the other GSs.. To
quantify the performance of the proposed setup, we consider an aggregate
channel model consisting of atmospheric turbulence and atmospheric
attenuation, where we obtained two performance indicators, outage probability
and ergodic capacity. We also consider the scenario of a constellation of
satellites in which the best satellite communicates with the best GS by using
opportunistic scheduling to maximize the site diversity. More precisely, the
paper makes the following specific contributions:

  * •

We focus on the physical layer performance of the downlink optical SatCom,
where the best GS is selected among a set of $\mathcal{K}$ sites that are
available for communication, to minimize the outage probability and to
maximize the ergodic capacity. For this network, we consider an aggregate
channel model consisting of turbulence induced fading, and atmospheric
attenuation due to Mie scattering and geometrical scattering, and we obtain
the instantaneous signal-to-noise ratio (SNR).

  * •

After obtaining the instantaneous SNR, we derive new closed-form outage
probability and ergodic capacity expressions to characterize the overall
performance of the proposed scheme. We further elaborate the system at high
SNR to obtain the overall diversity order.

  * •

We provide two new GS deployment scenarios: ground level deployment and high
ground windy weather deployment. We also provide some interesting system
design guidelines and consider aperture averaging to mitigate the adverse
effects of turbulence induced fading.

  * •

We also obtain the diversity order when a set of satellites form a
constellation and communicate with the best GS by using opportunistic
scheduling.

The remainder of the paper is organized as follows: In Section II, the system
model and problem formulation are outlined. In Section III, the outage
probability and ergodic capacity analyses are presented, and the overall
diversity order of the considered multi-site system is obtained. In Section
IV, numerical results are provided, and Section V concludes the paper.

Figure 1: Optical downlink SatCom with multiple ground stations.

##  II System Model and Problem Formulation

We consider here a downlink optical SatCom where a low Earth orbit (LEO)
satellite deployed on a circular orbit at $500$ km altitude seeks to
communicate with the best GS selected among a set of $\mathcal{K}$ sites that
are available for communication, as depicted in Fig. 1. In the proposed model,
we assume that all sites are in fixed locations providing perfect LOS
connectivity, and the best GS with the highest signal-to-noise (SNR) ratio is
selected to enhance the system performance by creating site diversity [25].
Prior to the data transmission, the satellite, which is moving at a speed of
$8$ km/s, takes aim at the best GS, and the GS is aligned with the incoming
beam to compensate for the bore-sight pointing errors. In this setup, the
aggregated channel model is taken into consideration which consists of
atmospheric attenuation ($I_{j}^{(a)}\big{|}_{j=1}^{\mathcal{K}}$) and
turbulence induced fading ($I_{j}^{(t)}\big{|}_{j=1}^{\mathcal{K}}$).
Mathematically speaking, the aggregated channel of the $j$-th GS can be
expressed as [26, 27]

$\displaystyle I_{j}=I_{j}^{(a)}I_{j}^{(t)},$ (1)

and the instantaneous SNR ($\gamma_{j}$) can be expressed as

$\displaystyle\gamma_{j}=\frac{P_{S}}{N_{0}}I_{j}^{2},$ (2)

where $P_{S}$ is the transmitted satellite signal power, and $N_{0}$ is the
one sided noise power spectral density. The following subsections present the
atmospheric attenuation and turbulence induced fading models for the proposed
model.

###  II-A Atmospheric Attenuation Model

In downlink optical SatCom, two different scattering effects can be observed
as the beam propagates to the GS. The first is Mie scattering, which redirects
the transmitted signal from its intended direction when the signal wavelength
is equal to the diameter of the particles in the medium. The second is
geometrical scattering, which causes reflection, refraction, and scattering
when the size of the particles in the medium is much greater than the signal
wavelength333It is important to note that Rayleigh scattering may also
adversely affect the optical SatCom. However, it can be negligible for systems
operating below $375$ THz, as per the ITU-R report [7]..

####  II-A1 Atmospheric Attenuation Due to Mie Scattering

Mie scattering is specified as the primary source of losses in downlink
optical SatCom operating between $150$ and $375$ THz frequencies
($\lambda=0.8-2$ $\mu$m) and it is largely caused by microscopic particles of
water [7, Sec. (3.1)]. The following expression, which can precisely model the
Mie scattering effects, is appropriate for GSs located at altitudes between
$0$ and $5$ km above the mean sea level [7]:

$\displaystyle\rho{{}^{\prime}}=ah_{E}^{3}+bh_{E}^{2}+ch_{E}+d,$ (3)

where $\rho{{}^{\prime}}$ denotes the extinction ratio, $h_{E}$ stands for the
height of the GS above the mean sea level (km), and $a$, $b$, $c$ and $d$ are
the wavelength $\lambda$ ($\mu$m)-dependent empirical coefficients, which can
be expressed as

$\displaystyle a=-0.000545\lambda^{2}+0.002\lambda-0.0038$ $\displaystyle
b=0.00628\lambda^{2}-0.0232\lambda+0.0439$ $\displaystyle
c=-0.028\lambda^{2}+0.101\lambda-0.18$ $\displaystyle
d=-0.228\lambda^{3}+0.922\lambda^{2}-1.26\lambda+0.719,$ (4)

and the atmospheric attenuation due to Mie scattering
$\Big{(}I_{j}^{(m)}\Big{)}$ can be expressed as444 Throughout the paper, we
assume that the height of each GS above the mean sea level is the same. We
also assume that all GSs are propagating at the same wavelength. However, the
elevation angles and the propagation distances may vary depending on the
location of each GS.

$\displaystyle
I_{j}^{(m)}=\exp\Big{(}-\frac{\rho{{}^{\prime}}}{\sin(\theta_{j})}\Big{)},$
(5)

where $\theta_{j}$ is the elevation angle of the $j$-th GS.

TABLE I: Geometrical scattering parameters for various types of clouds at $1550$ nm.  Cloud type |  $N$ (cm-3) |  $\mathcal{L}_{W}$ ($g/m^{-3}$) |  $V$ (km)
---|---|---|---
Cumulus | $250$ | $1.0$ | $0.0280$
Stratus | $250$ | $0.29$ | $0.0626$
Stratocumulus | $250$ | $0.15$ | $0.0959$
Altostratus | $400$ | $0.41$ | $0.0369$
Nimbostratus | $200$ | $0.65$ | $0.0429$
Cirrus | $0.025$ | $0.06405$ | $64.66$
Thin cirrus | $0.5$ | $3.128\times 10^{-4}$ | $290.69$

####  II-A2 Atmospheric Attenuation due to Geometrical Scattering

Geometrical scattering is used to model the attenuation that is close to the
surface of the Earth and is caused by fog or dense clouds. In this model,
visibility ($V_{j}$), which is an important factor for determining geometrical
scattering, can be modeled in terms of liquid water content
($\mathcal{L}_{W_{j}}$) and cloud number concentration ($N_{j}$) as [26]

$\displaystyle V_{j}=\frac{1.002}{(\mathcal{L}_{W_{j}}N_{j})^{0.6473}}.$ (6)

These parameters are summarized in Table I for various cloud formations. In
geometrical scattering, the attenuation can be expressed by using the Beer-
Lambert law,

$\displaystyle I_{j}^{(g)}=\exp(-\Theta_{j}L_{j}),$ (7)

where $L_{j}$ is the propagation distance, and $\Theta_{j}$ is the attenuation
coefficient, which can be expressed as [28, Sect. (3)]

$\displaystyle\Theta_{j}=\Bigg{(}\frac{3.91}{V_{j}}\Bigg{)}\Bigg{(}\frac{\lambda}{550}\Bigg{)}^{-\psi_{j}},$
(8)

where $\psi$ is the particle size related coefficient given according to Kim’s
model as

$\displaystyle\psi_{j}=\begin{cases}1.6,&V_{j}>50\\\ 1.3&6<V_{j}<50\\\
0.16V_{j}+0.34&1<V_{j}<6\\\ V-0.5&0.5<V_{j}<1\\\ 0&V_{j}<0.5,\end{cases}$ (9)

and the atmospheric attenuation $I_{j}^{(a)}$ for the $j$-th GS can be
expressed as [18]

$\displaystyle
I_{j}^{(a)}=I_{j}^{(g)}I_{j}^{(m)}=\exp(-\sigma_{j}L_{j})\exp\Big{(}-\frac{\rho^{\prime}}{\sin(\theta_{j})}\Big{)}.$
(10)

TABLE II: List of Notations and Parameters Parameter | Definition
---|---
$D_{G}$ | Hard receiver aperture diameter (m) [t]
$\mathcal{K}$ | Set of GSs that are available for communication [t]
$L$ | Propagation distance [t]
$\zeta$ | Zenith angle [t]
$\theta$ | Elevation angle [t]
$v_{r}$ | Ground r.m.s. wind speed (m/s) [t]
$v_{g}$ | Ground wind speed (m/s) [t]
$h_{E}$ | Height of the GS above mean sea level (km) [t]
$h$ | Altitude [t]
$H$ | Altitude of the satellite (m [t] )
$h_{0}$ | Height of the GS above ground level (m) [t]
$\rho_{c}$ | Atmospheric correlation width [t]
$\lambda$ | Wavelength [t]
$\alpha,\beta$ | Shape parameters of the EW fading [t]
$\eta$ | Fading severity parameter of the EW fading [t]
$k$ | Optical wave number [t]
$V$ | Visibility (km) [t]
$\psi$ | Particle size related coefficient [t]
$N$ | Cloud number concentration [t]
$\Theta$ | Attenuation coefficient [t]
$\mathcal{L}_{W}$ | Liquid water content [t]
$\sigma_{R}^{2}$ | Rytov variance [t]
$C_{n}^{2}$ | Refractive index constant [t]
$\sigma_{I}^{2}$ | Scintillation index [t]
$\gamma_{th}$ | Predefined threshold for acceptable communication quality [t]

###  II-B Turbulence Induced Fading Model

In this paper, the turbulence induced downlink channel experiences an
exponentiated Weibull fading channel, where the corresponding probability
density function (PDF) and cumulative distribution function (CDF) can be
expressed as [29]

$\displaystyle f_{I_{j}}^{(t)}(I)$
$\displaystyle=\frac{\alpha_{j}\beta_{j}}{\eta_{j}}\Bigg{(}\frac{I}{\eta_{j}}\Bigg{)}^{\beta_{j}-1}\exp\Bigg{[}-\Bigg{(}\frac{I}{\eta_{j}}\Bigg{)}^{\beta_{j}}\Bigg{]}$
$\displaystyle\times\Bigg{(}1-\exp\Bigg{[}-\Bigg{(}\frac{I}{\eta_{j}}\Bigg{)}^{\beta_{j}}\Bigg{]}\Bigg{)}^{\alpha_{j}-1}$
(11)

and

$\displaystyle{F}_{I_{j}}^{(t)}(I)=\Bigg{(}1-\exp\Bigg{[}-\Bigg{(}\frac{I}{\eta_{j}}\Bigg{)}^{\beta_{j}}\Bigg{]}\Bigg{)}^{\alpha_{j}},$
(12)

where $\alpha_{j}$, $\beta_{j}$ are the shape parameters and $\eta_{j}$ is the
scale parameter of the $j$-th GS. The expressions for $\alpha_{j}$,
$\beta_{j}$, and $\eta_{j}$ can be expressed as [30]

$\displaystyle\alpha_{j}=\frac{7.220\times\sigma_{I_{j}}^{2/3}}{\Gamma\Big{(}2.487\sigma_{I_{j}}^{2/6}-0.104\Big{)}},$
$\displaystyle\beta_{j}=1.012\Big{(}\alpha\sigma_{I_{j}}^{2}\Big{)}^{-13/25}+0.142$
$\displaystyle\eta_{j}=\frac{1}{\alpha\Gamma\Big{(}1+1/\beta_{j}\Big{)}g_{1}(\alpha_{j},\beta_{j})},$
(13)

where $g_{1}(\alpha_{j},\beta_{j})$ is the $\alpha$ and $\beta$ dependent
constant variable, which can be written as [30]

$\displaystyle
g_{1}(\alpha_{j},\beta_{j})=\sum_{k=0}^{\infty}\frac{(-1)^{k}\Gamma(\alpha_{j})}{k!(k+1)^{1+1/\beta_{j}}\Gamma(\alpha_{j}-k)},$
(14)

and $\sigma_{I_{j}}^{2}$ denotes the scintillation index of the $j$-th GS,
which can be given by [31, Sect. (12)]

$\displaystyle\sigma_{I_{j}}^{2}=\exp\Bigg{[}\frac{0.49\sigma_{R}^{2}}{\big{(}1+1.11\sigma_{R_{j}}^{12/5}\big{)}^{7/6}}+\frac{0.51\sigma_{R_{j}}^{2}}{\big{(}1+0.69\sigma_{R_{j}}^{12/5}\big{)}^{5/6}}\Bigg{]}-1,$
(15)

and the Rytov variance $\sigma_{R_{j}}^{2}$ can be expressed as [31, Sect.
(12)]

$\displaystyle\sigma_{R_{j}}^{2}=2.25k^{7/6}\sec^{11/6}(\zeta_{j})\int_{h_{0}}^{H}C_{n_{j}}^{2}(h)(h-h_{0})^{5/6}dh,$
(16)

where $k=\frac{2\pi}{\lambda}$ is the optical wave number, $\zeta_{j}$ is the
zenith angle of the $j$-th GS, $h_{0}$ stands for the height of the GS above
ground level, $H$ is the altitude of the satellite, and $C_{n_{j}}^{2}(h)$ is
the altitude ($h$) dependent refractive index constant, which can be written
as [32]

$\displaystyle C_{n_{j}}^{2}(h)$ $\displaystyle=8.148\times
10^{-56}v_{r_{j}}^{2}h^{10}e^{-h/1000}+2.7\times 10^{-16}e^{-h/1500}$
$\displaystyle+C_{0}e^{-h/100}\hskip 8.5359ptm^{-{2}/{3}},$ (17)

where $v_{r_{j}}=\sqrt{v_{g_{j}}^{2}+30.69v_{g_{j}}+348.91}$ is the r.m.s
ground wind speed in $m/s$ , $v_{g_{j}}$ is the ground wind speed in $m/s$ for
the $j$-th GS, and $C_{0}=1.7\times 10^{-14}$ is the nominal value of the
refractive index constant at ground level. Table II summarizes all notations
and parameters.

###  II-C Problem Formulation

This section formulates the GS selection strategy for the proposed setup. In
what follows, the GS with the highest instantaneous SNR is selected to
maximize the site diversity. Mathematically speaking, it can be formulated as

$\displaystyle j^{*}=\arg\max_{1\leq
j\leq\mathcal{K}}\Big{[}\gamma_{j}\Big{]},$ (18)

where $j^{*}$ is the selected GS index. By doing so, outage probability
($P_{\text{out}}$) can be minimized as

$\displaystyle
P_{\text{out}}=\Pr[\gamma\leq\gamma_{th}]=\Pr\bigg{[}\max_{1\leq
j\leq\mathcal{K}}(\gamma_{j})\leq\gamma_{th}\bigg{]},$ (19)

where $\gamma_{th}$ is the predefined threshold for acceptable communication
quality. Furthermore, ergodic capacity, which can assess the ergodic channel
capacity, can be formulated as [33]

$\displaystyle{\mathcal{C}}_{\text{erg}}={\mathbb{E}}\left[\log_{2}\Big{(}1+{\max_{1\leq
j\leq\mathcal{K}}(\gamma_{j})}\Big{)}\right],$ (20)

where ${\mathbb{E}}\left[\cdot\right]$ is the expectation operation. It is
important to note that it is almost impossible to find an exact ergodic
capacity expression for the proposed scenario by using (20). Therefore, with
the aid of Jensen’s inequality, we propose two approximate ergodic capacity
bounds as given here;

$\displaystyle{\mathcal{C}}_{\text{erg}}^{\text{B}_{1}}\approx\max_{1\leq
j\leq\mathcal{K}}{\mathbb{E}}\left[\log_{2}(1+{\gamma_{j}})\right],$ (21a)
$\displaystyle{\mathcal{C}}_{\text{erg}}^{\text{B}_{2}}\approx\log_{2}\Big{(}{\mathbb{E}}\left[1+\max_{1\leq
j\leq\mathcal{K}}({\gamma_{j}})\right]\Big{)},$ (21b)

##  III Performance Analysis

This section derives new closed-form outage probability and ergodic capacity
expressions for the proposed system.

###  III-A Outage Probability Analysis

Outage probability (OP) can be defined as the probability that the SNR will
fall below a predefined threshold, $\gamma_{th}$, for acceptable communication
quality. By substituting (12) into (19), with the aid of (18), the OP can be
expressed as

$\displaystyle
P_{\text{out}}=\prod_{j=1}^{\mathcal{K}}\Bigg{(}1-\exp\Bigg{[}-\Bigg{(}\frac{\gamma_{th}}{\Big{(}\eta_{j}I_{j}^{(a)}\Big{)}^{2}\bar{\gamma}_{j}^{(t)}}\Bigg{)}^{\beta_{j}/2}\Bigg{]}\Bigg{)}^{\alpha_{j}},$
(22)

where
$\bar{\gamma}_{j}^{(t)}=\frac{P_{S}}{N_{0}}\mathbb{E}\Big{[}\Big{(}I_{j}^{(t)}\Big{)}^{2}\Big{]}$
is the average SNR. By applying the Binomial theorem, and after a few
manipulations, a tractable OP expression can be found as

$\displaystyle
P_{\text{out}}=\prod_{j=1}^{\mathcal{K}}\sum_{\rho=0}^{\infty}\binom{\alpha_{j}}{\rho}(-1)^{\rho}\exp\Bigg{[}-\rho\Bigg{(}\frac{\gamma_{th}}{\Big{(}I_{j}^{(a)}\eta_{j}\Big{)}^{2}\bar{\gamma}_{j}^{(t)}}\Bigg{)}^{\frac{\beta_{j}}{2}}\Bigg{]}.$
(23)

To gain further insights about the system behavior, the OP can be analyzed at
high SNR. To do so, we first invoke the high SNR assumption of
$\exp(-x/a)\approx 1-x/a$ into (22) as

$\displaystyle
P_{\text{out}}^{\infty}=\prod_{j=1}^{\mathcal{K}}\Bigg{[}\Bigg{(}\frac{\gamma_{th}}{\Big{(}\eta_{j}I_{j}^{(a)}\Big{)}^{2}\bar{\gamma}_{j}^{(t)}}\Bigg{)}^{\alpha_{j}\beta_{j}/2}\Bigg{]}.$
(24)

Then, if we express $\bar{\gamma}_{j}^{(t)}=\kappa_{j}\bar{\gamma}$, where
$\kappa_{j}\big{|}_{j=1}^{\mathcal{K}}$ is constant, after a few
manipulations, the above expression can be written as

$\displaystyle
P_{\text{out}}^{\infty}=\prod_{j=1}^{\mathcal{K}}\Bigg{[}\Bigg{(}\frac{1}{\Big{(}\eta_{j}I_{j}^{(a)}\Big{)}^{2}\kappa_{j}}\Bigg{)}^{\alpha_{j}\beta_{j}/2}\Bigg{]}\Bigg{(}\frac{\gamma_{th}}{\bar{\gamma}}\Bigg{)}^{\sum_{j=1}^{\mathcal{K}}\alpha_{j}\beta_{j}/2}.$
(25)

At high SNR, the OP can be expressed as [34]

$\displaystyle
P_{\text{out}}^{\infty}=\mathcal{G}_{c}(\bar{\gamma})^{-\mathcal{G}_{d}},$
(26)

where $\mathcal{G}_{c}$ determines the shift of the curve in SNR. The
diversity order $\mathcal{G}_{d}$ can be defined as the slope of the OP curve.
For the considered multi-site system, $\mathcal{G}_{d}$ can be obtained as

$\displaystyle\mathcal{G}_{d}={\sum_{j=1}^{\mathcal{K}}\alpha_{j}\beta_{j}/2}$
(27)

If we consider a practical system setup in which a set of $\mathcal{Z}$
satellites form a constellation and communicate with the best GS, which
provides the highest SNR, by using opportunistic scheduling in which the best
satellite providing best channel characteristics scheduled to communicate with
the selected GS depending on the weather conditions, turbulence induced fading
and LOS characteristics. In that case, the diversity order of the multi-
satellite system can be obtained as

$\displaystyle\mathcal{G}_{d}^{\text{C}}=\sum_{z=1}^{\mathcal{Z}}\sum_{j=1}^{\mathcal{K}}\alpha_{j}\beta_{j}/2\approx\mathcal{Z}\mathcal{K}\alpha\beta/2.$
(28)

Proof: Please see Appendix A.

Furthermore, in this section, we propose two different practical deployment
scenarios and investigate the OP performance of the proposed setup when there
are $\mathcal{K}=20$ set of GSs that are available for communication.

####  III-A1 Case 1 - Ground level deployment scenario

In the first setup, we assume that all GSs are deployed at the ground level
($h_{0}=h_{E}=0$ m) above the mean sea and ground level, affecting from the
nominal ground r.m.s. wind speed, which is $v_{g}=2.8$ m/s. Assuming that the
average SNR is $\bar{\gamma}=24$ dB, and $\lambda=1550$ nm, the OPs that are
required for acceptable communication quality ($\gamma_{th}=7$ dB) can be
obtained as given in Case 1, Table III.

TABLE III: Outage probabilities for ground level deployment and high ground windy weather deployment scenarios Scenario | $\zeta$ | $N$ | $L_{W}$ | Outage probability
---|---|---|---|---
Case 1 | $40^{\circ}$ | $0.5$ | $3.128\times 10^{-4}$ | $0.8166$
$30^{\circ}$ | $0.5$ | $3.128\times 10^{-4}$ | $0.1619$
$15^{\circ}$ | $0.5$ | $3.128\times 10^{-4}$ | $5.583\times 10^{-4}$
$0^{\circ}$ | $0.5$ | $3.128\times 10^{-4}$ | $2.165\times 10^{-5}$
Case 2 | $40^{\circ}$ | $0.5$ | $3.128\times 10^{-4}$ | 0.0986
$30^{\circ}$ | $0.5$ | $3.128\times 10^{-4}$ | $2.314\times 10^{-4}$
$15^{\circ}$ | $0.5$ | $3.128\times 10^{-4}$ | $2.27\times 10^{-9}$
$0^{\circ}$ | $0.5$ | $3.128\times 10^{-4}$ | $1.28\times 10^{-11}$

####  III-A2 Case 2 - High ground windy weather deployment scenario

In the second setup, all GSs are deployed at $h_{0}=1000$ m above from the
ground level, and $h_{E}=1200$ m above from the mean sea level, where the wind
is blowing at a speed of $11.176$ m/s. Assuming that the average SNR is
$\bar{\gamma}=24$ dB and $\lambda=1550$ nm, the OPs that are required for
acceptable communication quality can be summarized in Case 2, Table 3. We
believe that these deployment scenarios can be used to provide practical
insights about the system architecture, so that a system designer can get a
quick idea about the overall system performance.

Table III shows the OP results for thin cirrus cloud formations, when the the
satellite is on a $500$ km circular orbit. We can see from the table that,
even though the wind speed increases in the high ground deployment, the
overall OP performance of the proposed system enhances as both atmospheric
turbulence and atmospheric attenuation effects reduces due to lower
propagation distance. Furthermore, increasing the zenith angle shows that the
GSs are affected by the attenuation and atmospheric turbulence at higher
levels. Thereby, keeping the zenith angle small can boost the overall
performance and can enhance the overall diversity as can be observed from Fig.
2. Furthermore, when $\mathcal{Z}=600$ and $1000$ constellation of satellites
are used, diversity order can be further improved as can be observed from Fig.
3.

Figure 2: Diversity order versus number of ground stations for the ground
level deployment scenario. Figure 3: Diversity order versus number of ground
stations for a constellation of satellites in the ground level deployment
scenario.

###  III-B Ergodic Capacity Analysis

####  III-B1 First Approximate Bound on the Ergodic Capacity

Ergodic capacity, expressed in bits/channel in use, can be defined with the
aid of (21a) as

$\displaystyle{\mathcal{C}}_{\text{erg}}^{\text{B}_{1}}={\log_{2}(e)}\max_{1\leq
j\leq\mathcal{K}}\Bigg{\\{}\int_{0}^{\infty}\frac{1}{1+\gamma}\bar{F}_{\gamma_{j}}(\gamma)d\gamma\Bigg{\\}},$
(29)

where $\bar{F}_{\gamma_{j}}(\gamma)=1-F_{\gamma_{j}}(\gamma)$ is the
complementary CDF of $\gamma_{j}$, which can be given as

$\displaystyle\bar{F}_{\gamma_{j}}(\gamma)=\sum_{\rho=1}^{\infty}\binom{\alpha_{j}}{\rho}(-1)^{\rho+1}\exp\Bigg{[}-\rho\Bigg{(}\frac{\gamma}{\Big{(}\eta_{j}I_{j}^{(a)}\Big{)}^{2}\bar{\gamma}_{j}^{(t)}}\Bigg{)}^{\frac{\beta_{j}}{2}}\Bigg{]},$
(30)

then, by invoking $\bar{F}_{\gamma_{j}}(\gamma)$ into (29),
${\mathcal{C}}_{\text{erg}}$ can be expressed as

$\displaystyle{\mathcal{C}}_{\text{erg}}^{\text{B}_{1}}$
$\displaystyle={\log_{2}(e)}\max_{1\leq
j\leq\mathcal{K}}\Bigg{\\{}\sum_{\rho=1}^{\infty}\binom{\alpha_{j}}{\rho}(-1)^{\rho+1}\int_{0}^{\infty}\frac{1}{1+\gamma}$
$\displaystyle\times\exp\Bigg{[}-\rho\Bigg{(}\frac{\gamma}{\Big{(}\eta_{j}I_{j}^{(a)}\Big{)}^{2}\bar{\gamma}_{j}^{(t)}}\Bigg{)}^{\frac{\beta_{j}}{2}}\Bigg{]}d\gamma\Bigg{\\}}.$
(31)

To find the closed-form solution of the above integral, we first use the
identity of $\frac{1}{1+\gamma}=G_{1,1}^{1,1}\bigg{[}\gamma\hskip
2.84544pt\bigg{|}\begin{matrix}0\\\ 0\end{matrix}\bigg{]}$, and we can express
${\mathcal{C}}_{\text{erg}}$ as

$\displaystyle{\mathcal{C}}_{\text{erg}}^{\text{B}_{1}}$
$\displaystyle={\log_{2}(e)}\max_{1\leq
j\leq\mathcal{K}}\Bigg{\\{}\sum_{\rho=1}^{\infty}\binom{\alpha_{j}}{\rho}(-1)^{\rho+1}\int_{0}^{\infty}G_{1,1}^{1,1}\bigg{[}\gamma\hskip
2.84544pt\bigg{|}\begin{matrix}0\\\ 0\end{matrix}\bigg{]}$
$\displaystyle\times\exp\Bigg{[}-\rho\bigg{(}\frac{\gamma}{\Omega_{j}}\bigg{)}^{\frac{\beta_{j}}{2}}\Bigg{]}d\gamma\Bigg{\\}},$
(32)

where $\Omega_{j}=\Big{(}I_{j}^{(a)}\eta_{j}\Big{)}^{2}\bar{\gamma}_{j}^{(t)}$
and $G_{c,d}^{a,b}\bigg{[}\cdot\hskip 2.84544pt\bigg{|}\begin{matrix}\cdot\\\
\cdot\end{matrix}\bigg{]}$ denotes the Meijer-G function [35, eqn.
07.34.02.0001.01]. By changing variables in the above integration as
$\chi=\rho\bigg{(}\frac{\gamma}{\Omega_{j}}\bigg{)}^{\frac{\beta_{j}}{2}}$ and
by using the identity of $\exp(-x)=G_{0,1}^{1,0}\bigg{[}x\hskip
2.84544pt\bigg{|}\begin{matrix}-\\\ 0\end{matrix}\bigg{]}$, the ergodic
capacity can be written as

$\displaystyle{\mathcal{C}}_{\text{erg}}^{\text{B}_{1}}$
$\displaystyle={\log_{2}(e)}\max_{1\leq
j\leq\mathcal{K}}\Bigg{\\{}\sum_{\rho=1}^{\infty}\binom{\alpha_{j}}{\rho}(-1)^{\rho+1}\rho^{\frac{\beta_{j}}{2}}\Bigg{(}\frac{2\Omega_{j}}{\beta_{j}}\Bigg{)}$
$\displaystyle\times\int_{0}^{\infty}\chi^{\frac{2}{\beta_{j}}-1}G_{1,1}^{1,1}\bigg{[}\Omega_{j}\Bigg{(}\frac{\chi}{\rho}\Bigg{)}^{\frac{2}{\beta_{j}}}\hskip
2.84544pt\bigg{|}\begin{matrix}0\\\
0\end{matrix}\bigg{]}G_{0,1}^{1,0}\bigg{[}\chi\hskip
2.84544pt\bigg{|}\begin{matrix}-\\\ 0\end{matrix}\bigg{]}d\chi\Bigg{\\}},$
(33)

and with the aid of [35, eqn. 07.34.21.0012.01], the closed form solution of
the above expression can be obtained as

$\displaystyle{\mathcal{C}}_{\text{erg}}^{\text{B}_{1}}$
$\displaystyle={\log_{2}(e)}\max_{1\leq
j\leq\mathcal{K}}\Bigg{\\{}\sum_{\rho=1}^{\infty}\binom{\alpha_{j}}{\rho}(-1)^{\rho+1}\rho^{\frac{\beta_{j}}{2}}\Bigg{(}\frac{2\Omega_{j}}{\beta_{j}}\Bigg{)}$
$\displaystyle\times
H_{2,1}^{1,2}\bigg{[}\Omega_{j}\rho^{\frac{\beta_{j}}{2}}\hskip
2.84544pt\bigg{|}\begin{matrix}(0,1),(1-\frac{2}{\beta_{j}},\frac{2}{\beta_{j}})\\\
(0,1)\end{matrix}\bigg{]}\Bigg{\\}}.$ (34)

where $H_{p,q}^{m,n}\bigg{[}\cdot\hskip
2.84544pt\bigg{|}\begin{matrix}\cdot,\cdot\\\ \cdot,\cdot\end{matrix}\bigg{]}$
denotes the Fox H-function [36].

####  III-B2 Second Approximate Bound on the Ergodic Capacity

Another tight bound on the ergodic capacity can be obtained by using (21b).
First, recall that

$\displaystyle{\mathcal{C}}_{\text{erg}}^{\text{B}_{2}}\approx\log_{2}\Big{(}1+{\mathbb{E}}\Bigg{[}\underbrace{\max_{1\leq
j\leq\mathcal{K}}({\gamma_{j}})}_{\gamma}\Bigg{]}\Big{)}.$ (35)

Thereafter, $\mathbb{E}\left[\gamma\right]$ can be expressed as

$\displaystyle\mathbb{E}\left[\gamma\right]=\int_{0}^{\infty}(1-F_{\gamma}(\gamma))d\gamma,$
(36)

where $F_{\gamma}(\gamma)$ can be obtained very similarly to (23) after
changing $\gamma_{th}$ with $\gamma$, and it can be approximately expressed as

$\displaystyle
F_{\gamma}(\gamma)\approx\sum_{\rho=0}^{\infty}\binom{\mathcal{K}\alpha}{\rho}(-1)^{\rho}\exp\Bigg{[}-\rho\Bigg{(}\frac{\gamma_{th}}{\Big{(}I^{(a)}\eta\Big{)}^{2}\bar{\gamma}^{(t)}}\Bigg{)}^{\frac{\beta}{2}}\Bigg{]}.$
(37)

By substituting (37) into (36), and after few manipulations,
$\mathbb{E}\left[\gamma\right]$ can be obtained as

$\displaystyle\mathbb{E}\left[\gamma\right]=\sum_{\rho=1}^{\infty}\binom{\mathcal{K}\alpha}{\rho}(-1)^{\rho+1}\rho^{-2/\beta}\bar{\gamma}\Big{(}I^{(a)}\eta\Big{)}^{2}\Gamma(1+2/\beta).$
(38)

Finally, substituting (38) into (35),
${\mathcal{C}}_{\text{erg}}^{\text{B}_{2}}$ can be easily obtained.

###  III-C Aperture Averaging

In downlink optical SatCom, the impact of scintillation can be great enough to
limit the performance of GS receivers. To compensate for the impact of
scintillation, an important enabler is aperture averaging. In downlink
communication, aperture averaging takes place especially when the atmospheric
correlation width $\rho_{c_{j}}$, which describes the effective diameter of
the $j$-th aperture, is lower than the aperture diameter of the $j$-th GS
($D_{G_{j}}$), i.e, $\rho_{c_{j}}<D_{G_{j}}$. In this case, scintillation is
spatially averaged over the aperture to reduce the adverse effects of
scintillation. For the proposed setup, $\rho_{c_{j}}$ can be calculated as
[31, Sect. (12)]

$\displaystyle\rho_{c_{j}}\approx\sqrt{\frac{45\times
10^{3}\sec(\zeta_{j})}{k}},\hskip 5.69046pt\sigma_{R_{j}}^{2}<<1,\hskip
5.69046pt0\leq\zeta_{j}<50,$ (39)

where $\rho_{c_{j}}$ represents the diameter of the point-like aperture for
the GS. For example, when $\zeta_{j}=40^{\circ}$, $\rho_{c_{j}}\approx 1.204$
cm shows the point-like aperture size for the optical downlink SatCom
operating at $\lambda=1.55$ nm wavelength. In aperture averaging, the aperture
diameter dependent scintillation index can be expressed as [31, Sect. (12)]

$\displaystyle\sigma_{I_{j}}^{2}=8.7k^{7/6}(H-h_{0})^{5/6}\sec^{11/6}(\zeta_{j})$
$\displaystyle\times$
$\displaystyle\Re\Bigg{\\{}\int_{h_{0}}^{H}C_{n}^{2}(h)\Bigg{[}\Bigg{(}\frac{kD_{G}^{2}}{16L}+i\frac{h-h_{0}}{H-h_{0}}\Bigg{)}^{5/6}-\Bigg{(}\frac{kD_{G_{j}}^{2}}{16L_{j}}\Bigg{)}^{5/6}\Bigg{]}dh\Bigg{\\}},$
(40)

where $\Re\\{\cdot\\}$ represents the real-valued terms. For example,
considering $\zeta_{j}=40^{\circ}$, $H=5\times 10^{5}$ m, and thin cirrus
cloud formations for the ground level deployment scenario, $D_{g}=20$ cm
aperture size can boost the overall outage performance from $6.884\times
10^{-7}$ to $3.441\times 10^{-9}$ at $30$ dB average SNR.

TABLE IV: List of Parameters and Values Parameters | Values
---|---
$\mathcal{K}$ | 10, 20 [t]
$\zeta$ |  $15^{\circ}$, $40^{\circ}$, $50^{\circ}$ [t]
$\lambda$ |  $1.55$ nm [t]
$H$ |  $5\times 10^{5}$ m [t]
$h_{0}$ |  $0$ m, $1000$ m [t]
$h_{E}$ |  $0$ km, $1.2$ km [t]
$v_{g}$ |  $2.8$ m/s, $11.176$ m/s [t]
$\gamma_{th}$ |  $7$ dB [t]
$C_{0}$ |  $1.7\times 10^{-14}$ [t]
Figure 4: Outage probability performance of the proposed scheme for the ground
level deployment scenario. Figure 5: Outage probability performance of the
proposed scheme for the high ground windy weather deployment scenario.

##  IV Numerical Results

In this section, the theoretical results are first verified by a set of
simulations. Then, the ground level deployment scenario is compared with the
high ground windy weather deployment scenario in terms of outage probability
and ergodic capacity. Finally, aperture averaging is illustrated in terms of
outage probability, and important design guidelines are outlined for practical
downlink laser SatComs.

In all figures, the parameters are set as $\alpha_{j}=\alpha$,
$\beta_{j}=\beta$, $\eta_{j}=\eta$, $I_{j}^{(a)}=I^{(a)}$,
$I_{j}^{(t)}=I^{(t)}$ $\zeta_{j}=\zeta$, $\psi_{j}=\psi$,
$C_{n_{j}}^{2}(h)=C_{n}^{2}(h)$, $\sigma_{I_{j}}^{2}=\sigma_{I}^{2}$,
$\sigma_{R_{j}}^{2}=\sigma_{R}^{2}$ for notational brevity, without losing
generality. Furthermore, in the ground level deployment scenario, the
parameters are set as $h_{0}=h_{E}=0$ m, $v_{g}=2.8$ m/s, whereas to
demonstrate the high ground windy weather scenario, the parameters are set as
$h_{0}=1000$ m, $h_{E}=1.2$ km, and $v_{g}=11.176$ m/s. Finally, the SNR
threshold is set to $\gamma_{th}=7$ dB, and two different cloud forms, cirrus
and thin cirrus, are used in the simulations together with three different
$\zeta$ angles: $\zeta=15^{\circ},40^{\circ}$ and $50^{\circ}$. All these
parameters and values are illustrated in Table IV.

###  IV-A Verifications of the Theoretical Expressions

Here we verify the theoretical results with the simulations. As we can see in
Fig. 4 and 5, the theoretical outage probability results, which are shown with
dashed lines, are in good agreement with the marker symbols, which are
generated by the simulations. Furthermore, both figures show that the overall
outage performance of the proposed scheme can be enhanced remarkably by
decreasing the $\zeta$ angle as expected. Finally, asymptotic curves, which
are depicted with solid lines, show the overall diversity order, which is
close to $\mathcal{G}_{d}=10\alpha\beta\approx 70$ at both figures.

Figure 6: Ergodic capacity performance of the proposed scheme for the high
ground windy weather deployment scenario. Figure 7: Comparison of ground level
deployment with high ground windy weather deployment in terms of outage
probability. Figure 8: Comparison of ground level deployment with high ground
windy weather deployment in terms of ergodic capacity. Figure 9: The impact of
aperture averaging on the ground level deployment scenario in terms of outage
probability.

Fig. 6, on the other hand, illustrates the ergodic capacity performance of the
proposed scheme for the high ground windy weather deployment
scenario555Theoretical results that are evaluated by using bound 1 (shown with
dashed lines) can be easily obtained by using Fox-H function in well-known
software programs like MATHEMATICA or MATLAB. However, in a few cases, the
Fox-H function may not work properly due to fractional fading severity values.
For this reason, numerical integrations are used to verify the correctness of
the first bound in certain cases.. As we can see in the figure, theoretical
bound 2, which is shown with solid lines, provides a tight upper bound,
whereas theoretical bound 1, which are indicated with the dashed lines behaves
like a tight lower bound with the marker symbols, which are the results of the
simulations. Furthermore, thin cloud formations can bring up to $11$ bit/s/Hz
capacity gain at $50$ dB average SNR as the atmospheric attenuation increases
due to cirrus cloud formations, and results in losses on the ergodic capacity
performance of the proposed scheme.

###  IV-B Comparison of High Ground Windy Weather Deployment and Ground Level
Deployment Scenarios

Next, the ground level deployment scenario is compared with the high ground
windy weather deployment scenario in terms of outage probability and ergodic
capacity. Fig. 7 compares both schemes in terms of outage probability. As we
can see in the figure, the high ground windy weather scenario outperforms the
ground level counterpart both at thin cirrus and cirrus cloud formations, when
$\zeta=15^{\circ}$. Similarly, Fig. 8 shows the ergodic capacity performance
of high ground windy weather deployment and ground level deployment scenarios
in terms of zenith angle. As we can see from the figure, high ground windy
weather deployment scenario outperforms its counterpart almost for all zenith
angles. However, when $\zeta\rightarrow 60^{\circ}$, the gap between two
deployment scenarios are closing and the ergodic capacity performance of the
proposed scheme goes to $0$. Thereby, in downlink optical SatCom, almost no
capacity gain can be obtained when $\zeta>60^{\circ}$ as the optical signal is
affected more from dense clouds, wind, and atmospheric turbulence induced
fading.

###  IV-C Impact of Aperture Averaging

This section shows the impact of aperture averaging for the ground level
deployment scenario in terms of outage probability. Fig. 9 shows that
increasing the receiver aperture diameter up to $20$ cm increases the
diversity order and enhances the overall performance depending on the zenith
angle and weather conditions. Also, interestingly, choosing the optimum zenith
angle can yield a higher performance gain than aperture averaging.

###  IV-D Constellation of Satellites

This section illustrates the outage probability by considering a constellation
of satellites that are communicating with the best GS by using opportunistic
scheduling. In Fig. 10, outage probability is depicted for
$\mathcal{\zeta}=1000$ constellation of satellites by considering ground level
deployment scenario. As can be observed from the figure, a constellation of
satellites can bring up to $7$ dB average SNR gain at about $10^{-6}$ outage
probability. However, the potential traffic congestion on the inter-satellite
communication and the data delivery latency should be taken into
consideration.

Figure 10: Outage probability performance of the proposed scheme for a
$\mathcal{Z}=1000$ constellation of satellites in the ground level deployment
scenario.

###  IV-E Design Guidelines

Finally, we provide some important design guidelines that can be helpful in
the design of downlink laser SatCom.

  * •

The simulations have shown that placing the GS to higher ground reduces the
adverse atmospheric conditions and enhances the overall performance. So, the
altitude of the GSs is of utmost importance in the design of downlink laser
SatComs.

  * •

The zenith angle has a direct impact on the design of downlink laser SatComs
as the effect of atmospheric turbulence is much lower when the zenith angle is
a small value. For this reason, the zenith angle should be taken into
consideration in the design of downlink laser SatCom.

  * •

Aperture averaging can be an important enabler for enhancing the overall
performance of the downlink laser SatComs, especially in the presence of
adverse weather conditions, and when the zenith angle is higher than
$\zeta=30^{\circ}$.

  * •

The simulations have shown that less than $10^{-10}$ outage probability can be
achieved at $30$ dB when a set of $20$ GSs is used. This shows that creating
site diversity through GS selection can be of utmost importance in downlink
laser SatComs to enhance the overall performance.

  * •

The simulations have shown that almost no capacity gain can be obtained when
the zenith angle is greater than $60^{\circ}$ as the system is exposed to
higher atmospheric turbulence induced fading and atmospheric attenuation due
to increased propagation distance.

  * •

Forming a constellation of satellites can further enhance the performance gain
up to $7$ dB at $10^{-5}$ outage probability and diversity order up to a
thousand times when the constellations consists of $1000$ satellites. However,
the potential traffic congestion on the inter-satellite communication and the
data delivery latency should be taken into consideration in the design of
satellite networks.

##  V Conclusion

This paper has focused on downlink laser SatComs, where the best GS is
selected from among a set of $\mathcal{K}$ candidates to provide fully
reliable connectivity and maximum site diversity. For the proposed structure,
outage probability and ergodic capacity expressions were derived and
asymptotic outage probability was conducted to provide the overall diversity
order. Furthermore, we considered two different ground station deployment
scenarios and investigated the impact of aperture averaging in terms of outage
probability. Finally, important system design guidelines were provided to help
in the design of downlink laser SatComs. The results have shown that the
altitude of the GS and the zenith angle are of utmost importance in downlink
laser SatCom.

## Appendix A: SITE DIVERSITY FOR CONSTELLATION OF SATELLITES

Let us consider that a set of $\mathcal{Z}$ satellites form a constellation
and communicating with the best GS which provides the highest signal-to-noise
(SNR) with the aid of opportunistic scheduling. In that case, the best GS can
be selected as

$\displaystyle\\{z^{*},j^{*}\\}=\arg{\max_{\begin{subarray}{c}{1\leqslant
j\leqslant\mathcal{K},}\\\ {1\leqslant
z\leqslant\mathcal{Z}}\end{subarray}}}\Big{[}\gamma_{j,z}\Big{]},$ (41)

where $j^{*}$ is the selected GS index and $z^{*}$ is the index of the
selected satellite, and the outage probability can be minimized as

$\displaystyle
P_{\text{out}}=\Pr[\gamma\leq\gamma_{th}]=\Pr\bigg{[}{\max_{\begin{subarray}{c}{1\leqslant
j\leqslant\mathcal{K},}\\\ {1\leqslant
z\leqslant\mathcal{Z}}\end{subarray}}}(\gamma_{j,z})\leq\gamma_{th}\bigg{]},$
(42)

By substituting (12) into (42) with the aid of (41), with few manipulations,
the OP can be expressed as

$\displaystyle P_{\text{out}}$
$\displaystyle=\prod_{z=1}^{\mathcal{Z}}\prod_{j=1}^{\mathcal{K}}\sum_{\rho=0}^{\infty}\binom{\alpha_{j,z}}{\rho}(-1)^{\rho}$
$\displaystyle\times\exp\Bigg{[}-\rho\Bigg{(}\frac{\gamma_{th}}{\Big{(}I_{j,z}^{(a)}\eta_{j,z}\Big{)}^{2}\bar{\gamma}_{j,z}^{(t)}}\Bigg{)}^{\frac{\beta_{j,z}}{2}}\Bigg{]},$
(43)

where
$\bar{\gamma}_{j,z}^{(t)}=\frac{P_{S}}{N_{0}}\mathbb{E}\Big{[}\Big{(}I_{j,z}^{(t)}\Big{)}^{2}\Big{]}$
is the average SNR. To obtain the asymptotic behavior of the outage
probability, we invoke the high SNR assumption of $\exp(-x/a)\approx 1-x/a$
into (43) as

$\displaystyle
P_{\text{out}}^{\infty}=\prod_{z=1}^{\mathcal{Z}}\prod_{j=1}^{\mathcal{K}}\Bigg{[}\Bigg{(}\frac{\gamma_{th}}{\Big{(}\eta_{j,z}I_{j,z}^{(a)}\Big{)}^{2}\bar{\gamma}_{j,z}^{(t)}}\Bigg{)}^{\alpha_{j,z}\beta_{j,z}/2}\Bigg{]},$
(44)

and after several manipulations, the diversity order of the considered multi-
site system can be obtained as given in (28). Hereto, the proof is completed.

## References

  * [1] D. R. Cheruku, _Satellite Communication_. I.K. International Publishing House Pvt Ltd, Delhi, India, 2010.
  * [2] P. Timothy, W. Bostian, and J. E. Allnutt, _Satellite Communication_. John Wiley & Sons, New York, USA, 2003.
  * [3] S. L. Kota, K. Pahlavan, and P. A. Leppänen, _Broadband Satellite Communications for Internet Access_. Kluwer Academic Publishers, Norwell, USA, 2003.
  * [4] “Starlink,” https://www.starlink.com, Accessed: 2020-08-20.
  * [5] “Oneweb,” http://www.oneweb.world, Accessed: 2020-08-20.
  * [6] H. Kaushal and G. Kaddoum, “Optical communication in space: Challenges and mitigation techniques,” _IEEE Communications Surveys & Tutorials_, vol. 19, no. 1, pp. 57–96, Aug. 2016.
  * [7] ITU-R, “Prediction methods required for the design of Earth-space systems operating between 20 THz and 375 THz,” International Telecommunication Union, Recommendation P.1622.
  * [8] A. N. Ince, _Digital Satellite Communications Systems and Technologies: Military and Civil Applications_. Springer Science & Business Media, Sep. 2012, vol. 186.
  * [9] D. Giggenbach, B. Epple, J. Horwath, and F. Moll, “Optical satellite downlinks to optical ground stations and high-altitude platforms,” in _Advances in Mobile and Wireless Communications_. Springer, Jul. 2008, pp. 331–349.
  * [10] V. V. Mai and H. Kim, “Mitigation of effects of angle-of-arrival fluctuation and pointing error on airborne free-space optical systems,” in _Optical Fiber Communication Conference_ , Mar. 2019, pp. W2A–40.
  * [11] D. Agarwal and A. Bansal, “Unified error performance of a multihop DF-FSO network with aperture averaging,” _IEEE/OSA Journal of Optical Communications and Networking_ , vol. 11, no. 3, pp. 95–106, Mar. 2019.
  * [12] C. Fuchs and F. Moll, “Ground station network optimization for space-to-ground optical communication links,” _IEEE/OSA Journal of Optical Communications and Networking_ , vol. 7, no. 12, pp. 1148–1159, Dec. 2015.
  * [13] M. S. Net, I. Del Portillo, E. Crawley, and B. Cameron, “Approximation methods for estimating the availability of optical ground networks,” _Journal of Optical Communications and Networking_ , vol. 8, no. 10, pp. 800–812, Oct. 2016\.
  * [14] N. K. Lyras, C. N. Efrem, C. I. Kourogiorgas, and A. D. Panagopoulos, “Optimum monthly based selection of ground stations for optical satellite networks,” _IEEE Communications Letters_ , vol. 22, no. 6, pp. 1192–1195, Jun. 2018\.
  * [15] S. Gong, H. Shen, K. Zhao, R. Wang, X. Zhang, T. De Cola, and J. A. Fraier, “Network availability maximization for free-space optical satellite communications,” _IEEE Wireless Communications Letters_ , vol. 9, no. 3, pp. 411–415, Mar. 2019.
  * [16] P. Gopal, V. K. Jain, and S. Kar, “Performance improvement of FSO satellite downlink using aperture averaging and receiver spatial diversity,” _IET Optoelectronics_ , vol. 10, no. 4, pp. 119–127, Jul. 2016.
  * [17] K. Li, J. Ma, A. Belmonte, L. Tan, and S. Yu, “Performance analysis of satellite-to-ground downlink optical communications with spatial diversity over Gamma–Gamma atmospheric turbulence,” _Optical Engineering_ , vol. 54, no. 12, pp. 7575–7585, Aug. 2015.
  * [18] S. Johari and V. Sundharam, “Performance analysis of IM/DD vs. heterodyne detection techniques of an earth-satellite FSO link for next generation wireless communication,” in _Malaysia International Conference on Communications (MICC)_ , Nov. 2017, pp. 191–196.
  * [19] A. Viswanath, V. K. Jain, and S. Kar, “Analysis of earth-to-satellite free-space optical link performance in the presence of turbulence, beam-wander induced pointing error and weather conditions for different intensity modulation schemes,” _IET Communications_ , vol. 9, no. 18, pp. 2253–2258, Dec. 2015.
  * [20] N. Alshaer, T. Ismail, and M. E. Nasr, “Performance evaluation and security analysis of ground-to-satellite FSO system with CV-QKD protocol,” _IET Communications_ , vol. 14, no. 10, pp. 1534–1542, Jun. 2020.
  * [21] N. Alshaer, T. Ismail, H. Seleem, and M. E. Nasr, “Analysis of beam wander and scintillation in ground-to-satellite FSO system with DPSK,” in _Novel Intelligent and Leading Emerging Sciences Conference (NILES)_ , vol. 1, Nov. 2019, pp. 5–8.
  * [22] H. D. Le, V. V. Mai, C. T. Nguyen, and A. T. Pham, “Throughput analysis of incremental redundancy hybrid ARQ for FSO-based satellite systems,” in _IEEE 90th Vehicular Technology Conference (VTC2019-Fall)_ , Nov. 2019, pp. 1–5.
  * [23] E. Illi, F. El Bouanani, F. Ayoub, and M.-S. Alouini, “A PHY layer security analysis of a hybrid high throughput satellite with an optical feeder link,” _IEEE Open Journal of the Communications Society_ , vol. 1, pp. 713 – 731, May 2020.
  * [24] E. Zedini, A. Kammoun, and M.-S. Alouini, “Performance of multibeam very high throughput satellite systems based on FSO feeder links with HPA nonlinearity,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 9, pp. 5908 – 5923, Jun. 2020.
  * [25] E. Erdogan, “Joint user and relay selection for relay-aided RF/FSO systems over exponentiated Weibull fading channels,” _Optics Communications_ , vol. 436, pp. 209–215, Apr. 2019.
  * [26] M. S. Awan, E. Leitgeb, B. Hillbrand, F. Nadeem, M. Khan _et al._ , “Cloud attenuations for free-space optical links,” in _International Workshop on Satellite and Space Communications_ , Sep. 2009, pp. 274–278.
  * [27] D. Dubey, Y. K. Prajapati, and R. Tripathi, “Performance enhancement of hybrid-sim for optical wireless downlink communication with aperture averaging and receiver diversity,” _IET Communications_ , vol. 14, no. 18, Nov. 2020.
  * [28] Z. Ghassemlooy, W. Popoola, and S. Rajbhandari, _Optical Wireless Communications: System and Channel Modelling with Matlab®_. CRC Press, 2019\.
  * [29] R. Barrios and F. Dios, “Exponentiated Weibull distribution family under aperture averaging for Gaussian beam waves,” _Optical Express_ , vol. 20, no. 12, pp. 13 055–13 064, May 2012.
  * [30] R. A. Barrios, “Exponentiated Weibull fading channel model in free-space optical communications under atmospheric turbulence,” Ph.D. dissertation, Universitat Politècnica de Catalunya (UPC), May 2013.
  * [31] L. Andrews and R. Phillips, _Laser Beam Propagation Through Random Media_ , ser. SPIE Press monograph. SPIE Press, 2005\.
  * [32] ITU-R, “Propagation data required for the design of Earth-space systems operating between 20 THz and 375 THz,” International Telecommunication Union, Recommendation P.1621-1.
  * [33] J. Choi and B. L. Evans, “Analysis of ergodic rate for transmit antenna selection in low-resolution ADC systems,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 1, pp. 952–956, Oct. 2018.
  * [34] Z. Wang and G. B. Giannakis, “A simple and general parameterization quantifying performance in fading channels,” _IEEE Transactions on Communications_ , vol. 51, no. 8, pp. 1389–1398, 2003.
  * [35] “From Wolfram Research: The mathematical functions site,” http://functions.wolfram.com, Accessed: 2020-06-18.
  * [36] A. M. Mathai, R. K. Saxena, and H. J. Haubold, _The H-Function: Theory and Applications_. Springer Science & Business Media, 2009.

|  Eylem Erdogan is currently an Associate Professor in the Department of
Electrical and Electronics Engineering, Istanbul Medeniyet University. He was
a Post-Doctoral Fellow in Electrical Engineering department, Lakehead
University, Thunder Bay, ON, Canada from March 2015 to September 2016 and a
visiting professor in Carleton University, Ottawa, Canada during summer 2019.
His research interests are in the broad areas of wireless communications, the
performance analysis of cooperative relaying in cognitive radio networks,
unmanned aerial vehicle communications and networks, free space optical
communications and optical satellite networks.
---|---
|  Ibrahim Altunbas is a Professor in the Electronics and Communication
Engineering Department, Istanbul Technical University. He was a Visiting
Researcher at Texas A&M University (College Station, TX, USA) in 2001. He was
a Postdoctoral Fellow and Visiting Professor at the University of Ottawa (ON,
Canada) in 2002 and 2015 summer, respectively. His current research interests
include spatial modulation, non-orthogonal multiple access, physical layer
security, drone communications, and reconfigurable intelligence surface-based
communication.
---|---
|  Gunes Karabulut Kurt received the B.S. degree with high honors in
electronics and electrical engineering from the Bogazici University, Istanbul,
Turkey, in 2000 and the M.A.Sc. and the Ph.D. degrees in electrical
engineering from the University of Ottawa, ON, Canada, in 2002 and 2006,
respectively. From 2000 to 2005, she was a Research Assistant with the CASP
Group, University of Ottawa. Between 2005 and 2006, she was with TenXc
Wireless, Canada. From 2006 to 2008, she was with Edgewater Computer Systems
Inc., Canada. From 2008 to 2010, she was with Turkcell Research and
Development Applied Research and Technology, Istanbul. Since 2010, she has
been with Istanbul Technical University, where she currently works as a
Professor. She is a Marie Curie Fellow. She is an adjunct research professor
at Carleton University, also currently serving an Associate Technical Editor
(ATE) of the IEEE Communications Magazine.
---|---
|  Michel Bellemare received is B.Eng. degree from Universitéde Sherbrooke,
Sherbrooke, Québec, Canada in Communications Engineering in 1986. He has
gained experience in terrestrial and space wireless communications in various
companies such as Nortel Networks, Ultra Electronics, SR Telecom, etc. He is
now a Space Systems Architect at MDA corporation, Sainte-Anne-de-Bellevue,
Québec, Canada.
---|---
|  Guillaume Lamontagne obtained his B.Eng and M.Eng degrees from École de
Technologie Supérieure (ÉTS), Montréal, Canada in 2007 and 2009 respectively.
His experience in satellite communications started through internships and
research activities at the Canadian Space Agency (CSA) in 2005 and the Centre
national d’études spatiales (Cnes), France, in 2006 and 2008. He joined MDA in
2009 and held various communication systems engineering and management
positions before being appointed Director of Technology, Payloads in 2019.
Through this role, he is leading MDA’s R&D activities for satellite
communications as well as establishing the related long term development
strategy.
---|---
|  Halim Yanikomeroglu (F’17) received the B.Sc. degree in electrical and
electronics engineering from the Middle East Technical University, Ankara,
Turkey, in 1990, and the M.A.Sc. degree in electrical engineering and the
Ph.D. degree in electrical and computer engineering from the University of
Toronto, Canada, in 1992 and 1998, respectively. He is a Professor with the
Department of Systems and Computer Engineering, Carleton University, Ottawa,
Canada. His research interests cover many aspects of 5G/5G+ wireless networks.
His collaborative research with industry has resulted in 37 granted patents.
Dr. Yanikomeroglu is a fellow of the IEEE, Engineering Institute of Canada
(EIC), and the Canadian Academy of Engineering (CAE). He is a Distinguished
Speaker for IEEE Communications Society and IEEE Vehicular Technology Society.
---|---