Correlations of photon trajectories in the problem of light scintillations
aa r X i v : . [ qu a n t - ph ] S e p Correlations of photon tra jectories in theproblem of light scintillations
R. A. Baskov and O. O. Chumak ∗ Institute of Physics of the National Academy of Sciencespr. Nauki 46, Kyiv-28, MSP 03028 UkraineSeptember 24, 2018
Abstract
A distribution function approach is applied to describe the dynam-ics of the laser beam in the Earth atmosphere. Using a formal solutionof the kinetic equation for the distribution function, we have devel-oped an iterative scheme for calculation of the scintillation index ( σ ).The problem reduces to obtaining the photon trajectories and theircorrelations. Bringing together theoretical calculations and many-foldcomputer integrations, the value of σ is obtained. It is shown thata considerable growth of σ in the range of a moderate turbulence isdue to the correlations of different trajectories. The criteria of appli-cability of our approach for both the coherent and partially coherentlight are derived. Basic principles of the radiation transfer theory were formulated in the sem-inal paper of Schuster [1] as early as the beginning of the 20th century. The ∗ Corresponding author: [email protected] C n is that it cannot be reliably computed from first principles. The refractive-index fluctuations arise from temperature inhomogeneities of the air. Theinhomogeneities cause turbulent eddies which give rise to a random distribu-tion of the air density([3]-[5]). This results in random spatial variations ofthe refractive index.The turbulent eddies are described by a wide range of characteristiclengths of inhomogeneities. These lengths cover the interval from few mil-limeters (the inner radius, l ) to hundred meters (the outer radius, L ).Therefore, various types of beam scattering are observed. The scatteringby large-size eddies results in random redirections of the beam as a whole.This process is known in the literature as a ”wandering” or ”dancing” of thebeam [6],[7]. On the other hand, the scattering by small-size eddies causesspreading of the beam. For a long-distance propagation or a strong turbu-lence, the beam radius becomes greater than the characteristic sizes of theinhomogeneities. In this case the probability of the beam to be redirectedbecomes small and the relative value of the wandering radius decreases [8].The beam wandering and broadening can be considered as the specificmanifestations of a more general phenomenon, namely the intensity fluctua-tions (i.e. scintillations) caused by the atmosphere turbulence. The scintil-lations have a tendency of saturating for a long-distance propagation [9],[10](the regime of a strong turbulence). This is because in the course of propaga-tion the radiation acquires the properties of the Gaussian statistics when thesignal-to-noise ratio (SNR) tends to unity. The asymptotic behavior of thescintillation index, σ →
1, was explained in Refs. [11]-[13]. Moreover, it wasshown quite generally that this property stays unchanged for any refractiveindex distribution, provided the response time of the recording instrument isshort compared with the source coherence time. This result was confirmedanalytically in [14].At the same time, calculations, performed by different methods in [15]2nd [16], show a possibility of significant suppression of the scintillations.To this end partially coherent laser beams with the coherence time shorterthan the detector integration time (a slow detector) can be used. The caseof a partial coherence was also studied in [17],[18]. Recent theoretical andexperimental developments on propagation of partially coherent beams in aturbulent atmosphere were discussed in [19].There are several analytical approaches explaining behavior of the scin-tillation index in the case of strong turbulence [16, 20, 21]. Their analysis isbased on the physical picture where four waves, forming the second momentof the intensity, conserve only pair correlations in course of long-distancepropagation. Two different pairs of the photon trajectories contribute intothe square of the photon density at the detector. Dashen used the Feynmanpath integrals to prove that in a convincing manner [20].The recent interest to beam propagation was awakened by the develop-ment of quantum communication in the free atmosphere [22], [23]. Detailedstudies of the effect of the turbulence-induced losses on the quantum state ofthe light in the course of satellite-mediated communication and for realiza-tion of the entanglement transfer in the atmosphere were reported in Refs.[24] and [25].The formalism of the photon distribution function (the photon density inthe coordinate-momentum space [26]), is also applicable to the problem ofscintillations [8, 16, 27, 28]. The mentioned papers are based on a physicalpicture, which is similar to the described above. The method of photon distri-bution function is used for description of both the classical and the quantumlight including propagation of single-photon pulses (see, for example, Refs.[27, 29, 30]). Solution of the kinetic equation for the operator of photondensity is based on the method of characteristics. The assumption of weakdisturbances of photon momenta by the atmosphere (the paraxial approxi-mation) reduces the problems of scintillations to the problem of obtainingphoton trajectories and their correlations. A slowly varying fluctuating force,deflecting photon trajectories from straight lines, describes the effect of theatmospheric eddies.In this work we study the scintillation index for moderate and strongturbulences, when correlation of trajectories of only two photons is required.Accuracy of the calculations depends on the accuracy of obtaining the tra-jectories. Using high-order iterations and bringing together analytical andnumerical procedures, we calculate the scintillation index. Our main interestis to analyze the range of moderate turbulence strengths where previous the-3ries do not ensure a reliable description. Comparison of the obtained resultswith those represented in [16] helps indicate the range of turbulence where asimplified approach should be corrected by high-order iterations. Also, ourstudies describe more realistically the effect of partial coherence.
The photon distribution function is defined by analogy with distributionfunctions in solid state physics. In particular, it is similar to the phonondistribution function. Both of them are defined as [26],[31] f ( r , q , t ) = 1 V X k e − i kr b † q + k / b q − k / , (1)where b † q and b q are the bosonic creation and annihilation operators of pho-tons or phonons with the momenta q , and V ≡ L x L y L z is the normalizingvolume. Polarization of the corresponding modes is not specified in (1).In the paraxial approximation, assumed here, the initial polarization of thebeam remains almost unchanged even for a long-distance propagation (see,for example, Ref. [32]).The operator f ( r , q , t )) describes the photon (phonon) density in thephase ( r,q ) space. Usually, the characteristic sizes of spatial inhomogeneitiesof the radiation field are much greater than the wave-length. In this casethe sum in Eq. (2) can be restricted by small k . Here and in what followswe consider that k < k ≪ q , where q is the wave vector corresponding tothe central frequency of the radiation, ω = cq . At the same time k shouldbe taken sufficiently large to provide a required accuracy of the beam profiledescription.The evolution of the Heisenberg operator f ( r , q , t ) is determined by thecommutator ∂ t f ( r , q , t ) = 1 i ~ [ f ( r , q , t ) , H ] , (2)where H = X q ~ ω q b † q b q − X q , k ~ ω q n k b † q b q + k , (3)is the Hamiltonian of photons in a medium with a fluctuating refractive index n ( r ) ( n k is its Fourie transform), ~ ω q = ~ cq and c q = ∂ω∂ q are the vacuumvalues of the photon energy and velocity, respectively.4ssuming the characteristic values of the photon momentum to be muchgreater than the wave vectors of turbulence, the kinetic equation for thephoton distribution function can be written as { ∂ t + c q ∂ r + F ( r ) ∂ q } f ( r , q , t ) = 0 , (4)where F ( r ) = ω ∂ r n ( r ) is the random force originating from the atmosphericturbulence. The general solution of Eq. (4) is given by f ( r , q , t ) = φ ( r − Z t dt ′ ∂ r ( t ′ ) ∂t ′ ; q − Z t dt ′ ∂ q ( t ′ ) ∂t ′ ) , (5)where the function φ ( r , q ) is the ”initial” value of f ( r , q , t ), i.e. φ ( r , q ) = 1 V X k e − i kr ( b + q + k / b q − k / ) | t =0 ≡ X k e − i kr φ ( k , q ) . (6)The derivatives ∂ r ( t ′ ) ∂t ′ and ∂ p ( t ′ ) ∂t ′ should satisfy the equations ∂ r ( t ′ ) ∂t ′ = c [ q ( t ′ )] ∂ q ( t ′ ) ∂t ′ = F [ r ( t ′ )] , (7)completed with the boundary conditions r ( t ′ ) = r and q ( t ′ ) = q for t = t ′ .As we see, Eqs. (7) coincide with the classical (the Newton) equations ofmotion of a point particle moving with the velocity c q and affected by anexternal force F ( r ). Formal solutions of Eqs. (7) can be written as q ( t ′ ) = q + Z t ′ t dt ′′ F [ r ( q , t ′′ )] (8)and r ( t ′ ) = r − c q ( t − t ′ ) − cq Z t ′ t dt ′′ ( t ′′ − t ′ ) F [ r ( t ′′ )] , (9)Eqs. (8) and (9) allow us to rewrite the expression (5) as f ( r , q , t ) = φ ( r − c q t + cq Z t dt ′ t ′ F [ r ( q , t ′ )]; q − Z t dt ′ F [ r ( q , t ′ )] ) . (10)5f F ( r ) is a known function, an approximate value for f ( r , q , t ) can be ob-tained by inserting the term r ( q , t ′ ) ≈ r − c q ( t − t ′ ) into Eq. (10). In thiscase the argument of the fluctuating force F [ r ( q , t ′ )] is replaced by a straightline, that is correct only in the absence of the turbulence. Improvement ofthe theory can be achieved if the argument of F accounts for the turbulence.It follows from Eq. (10) that statistical properties of the radiation dependnot only on the turbulence but also on the initial distribution function φ ( r , q ).This function is determined by the source field. Its explicit form is determinedin the course of ”sewing” of the near-aperture and the atmospheric fields [16]given by the amplitudes b q ( b † q ). We consider the light propagation in the z -direction. The source field is assumed to be described by the Gaussianfunction, Φ( r ) = (2 /π ) / r e − r ⊥ /r . Then the propagating amplitudes aregiven by b q ⊥ ,q ( t = 0) = b (2 π/S ) / r e − q ⊥ r / , (11)where b is the near-aperture amplitude of the laser field, index ( ⊥ ) means theperpendicular to the z -axis components, and S = L x L y .We will take into account the effect of the phase diffuser by multiplyingthe distribution Φ( r ) by the phase factor e − i ar ⊥ where the quantity a is arandom variable. In this case Eq. (11) should be modified by substitutingin its right-hand side q ⊥ + a ≡ q a for q ⊥ . Such a simple modeling of thephase diffuser is justified if (i) the detection time is much longer than thecharacteristic time of the variation of a (slow detector) and (ii) there is alarge root mean square of the phase fluctuations. (More detailed analysis ispresented in [28].) This case corresponds to the Gaussian distribution of a : P ( a x,y ) = λ π / e − a x,y λ / , (12)with a covariance h a x,y i = λ − and the transverse correlation function of theoutgoing field (at t = 0) is given by h E ( r ⊥ ) E ( r ⊥ + ∆ ) i a = E e − [ r ⊥ +( r ⊥ + ∆ ) ] r − e − ∆ λ − . (13)Here E = E ( r ⊥ = 0 , ∆ = 0 , t = 0) and the notation h ... i a means averag-ing over distribution P ( a x,y ). The radiation, whose correlation propertiesare described by function (13), is referred to as the Gaussian Shell-modelfield. The parameter λ in the exponential factor describes the decrease ofthe transverse correlation length. It can also be said that this parametergenerates a new characteristic length, 1 /r , in the momentum distribution6i.e., in the q -domain). This is seen from the explicit term for φ ( k , q ) whichafter averaging over the fluctuations of a reduces to h φ ( k , q ) i a = 2 π b † bV S r e − q ⊥ r − k ⊥ r , (14)where r = r (cid:18) r λ − (cid:19) − , and variables q z and k z are omitted.It is seen from Eq. (14) that q ⊥ is distributed in the range of the orderof √ /r that is greater than the one for coherent beam. In contrast, thecharacteristic value of ˜ k depends only on the initial size of the beam ( ˜ k ∼√ /r ).In the course of light propagation, the diffraction phenomena and scatter-ing by atmospheric inhomogeneities broaden the beam resulting in decreaseof ˜ k . At the same time, the value of ˜ q increases with the distance. This isbecause of the Brownian-like motion of photons in the q ⊥ -domain (see Ref.[16]). Such a simple physical picture, elucidating evolution of the beam ge-ometry, is, however, not applicable to the description of scintillations. Thephenomenon of scintillations is more complecated and can be described interms of spatio-temporal correlations of four waves. The photon distribution function is used here to obtain the scintillation index σ . The definition of σ is given by σ = h I ( r ) i − h I ( r ) i h I ( r ) i . (15)The photon density I ( r , t ) is expressed in terms of the distribution functionas I ( r , t ) = X q f ( r , q , t ) = 2 π b † br SV X q , k e − i k [ r − c ( q ) t + cq R t dt ′ t ′ F ( r ( q ,t ′ ))] − Q a r − k r , (16)where Q a ≡ Q + a = q + a − R t dt ′ F [ r ( q , t ′ ). The summation is taking over q ⊥ and k ⊥ components, while q z and k z are considered to be fixed: q z = q and k z = 0. The exponential term originates from the solution (10) of thekinetic equation (4). 7o obtain h I ( r , t ) i , three independent averagings are required. One ofthem concerns the source variables. In the case of a coherent state of thesource, | β i , we have h b † b i = | β | . The second averaging over a random phaseof the diffuser should be carried out as explained by Eq. (14). The thirdaveraging deals with the fluctuating force F . These three actions can beperformed independently that facilitates the analysis. Also, the calculationsare simplified if we use the identity e − Q r / ≡ Z d p πr e i pQ − p / r . (17)Because of Eq. (17), the term in the exponent of Eq. (16) reduces to the lin-ear in F form. Then, considering F as a random Gaussian variable, the valueof h I ( r , t ) i can be easily obtained in a manner, explained in Ref. [16]. Tocalculate h I ( r , t ) i , an explicit form of the refractive-index correlation func-tion, h n ( r ) n ( r ′ ) i , is required. In a statistically homogeneous atmosphere itcan be written as h n ( r ) n ( r ′ ) i = Z d g e − i g ( r − r ′ ) ψ ( g ) . (18)A widely used the von Karman approximation for the spectrum, ψ ( g ), isgiven by ψ ( g ) = 0 . C n exp[ − ( gl / π ) ][ g + L − ] / , | g | ≡ g, (19)where the vector g is defined in the three dimensional domain.The ”source” part of h I ( r ) i , given by h b † bb † b i , is approximately equal to h b † b † bb i = | β | , when the condition | β | >> | β | is satisfied. This inequalityimplies that the initial laser radiation is in a multiphoton coherent state. Theaveraging over independent random quantities a and a ′ can be used insteadof the time averaging of the diffuser state. Then we have h I ( r , t ) i = (cid:12)(cid:12)(cid:12)(cid:12) πβ r V S (cid:12)(cid:12)(cid:12)(cid:12) X q , k , q ′ , k ′ h e − i k [ r − c q t + cq R t dt ′ t ′ F ( r q ( t ′ ))] − i k ′ [ r − c q ′ t + cq R t dt ′ t ′ F ( r q ′ ( t ′ ))] × (cid:8) e − ( Q a + Q a ′ + k k ′ ) r + e − [( Q a + k ) +( Q ′ a − k ′ ) +( Q a ′ − k ) +( Q ′ a ′ + k ′ ) ] r (cid:9) i . (20)There are two terms in the braces of Eq. (20). They appear only if theinitial four-wave correlation reduces to the pair correlation [16]. Such a mod-ification of the statistical properties of the radiation occurs when the waves8ropagate for a long time which is sufficient for randomization of the trans-verse photon momentum. A more general case, which includes the regime offast detection, was analyzed in Ref. [28].The averaging of Eq. (20) over a and a ′ results in h I ( r , t ) i = (cid:12)(cid:12)(cid:12)(cid:12) πβ r V S (cid:12)(cid:12)(cid:12)(cid:12) X q , k , q ′ , k ′ h e − i { k [ r − c ( q ) t ]+ k ′ [ r − c ( q ′ ) t ]+ cq R t dt ′ t ′ [ kF ( r ( q ,t ′ ))+ k ′ F ( r ( q ′ ,t ′ ))] } × (cid:8) e − ( Q + Q ′ ) r / − ( k + k ′ ) r / + e − [( Q − Q ′ ) +( k + k ′ ) / r / − [( Q + Q ′ ) +( k − k ′ ) / r / (cid:9) i . (21)In the absence of a phase diffuser, r = r , the summands in the last bracescontribute equally into (21).Similarly to Eq. (17), the factor e − ( Q + Q ′ ) r in (21) can be expressed inthe integral form as e − ( Q + Q ′ ) r = Z d p d p ′ (2 πr ) e i pQ + i p ′ Q ′ − ( p + p ′ ) / r . (22)As we see, the exponent in the left-hand side is represented as a linear formof the force F . A similar transform is applicable to the second term in thelast braces of (21). As a result, the fluctuating force enters the right-handside of (21) only via the common multiplier, M , given by M = e − i R t dt ′ { ( p + k t ′ c/q ) F [ r ( q ,t ′ )]+( p ′ + k ′ t ′ c/q ) F [ r ( q ′ ,t ′ )] } . (23)Obtaining of the average value of I reduces to averaging of M with many-foldintegration. Assuming the exponent in (23) as a Gaussian random variable,we can write h M i = e − (cid:10)(cid:0) R t dt ′ { ( p + k t ′ c/q ) F [ r ( q ,t ′ )]+( p ′ + k ′ t ′ c/q ) F [ r ( q ′ ,t ′ )] } (cid:1) (cid:11) ≡ e − ( φ PP +2 φ PP ′ + φ P ′ P ′ ) . (24)Two types of the correlation functions determine h M i : φ P P ′ = Z t Z t dt ′ dt ′′ ( p + k t ′ c/q ) · h F [ r ( q , t ′ )] F [ r ( q ′ , t ′′ )] i · ( p ′ + k ′ t ′′ c/q )] , (25) φ P P = Z t Z t dt ′ dt ′′ ( p + k t ′ c/q ) · h F [ r ( q , t ′ )] F [ r ( q , t ′′ )] i · ( p + k t ′′ c/q )] , (26)9here symbols P and P ′ denote sets of three vector variables P = { q , p , k } and P ′ = { q ′ , p ′ , k ′ } . The correlation functions of the forces along different( q = q ′ ) and coinsiding ( q = q ′ ) trajectories enter Eqs. (25) and (26),respectively. The former can be rewritten as h F α [ r ( q , t ′ )] F β [ r ( q ′ , t ′′ )] i = h F α [ r ( q , t ′ ) − r ( q ′ , t ′′ )] F β [0] i , (27)where the notations α and β stand for the x and y - components. Theexpression for (26) follows from Eq. (27) by setting q = q ′ .The right-hand side of Eq. (27) is assumed to be a function of the coor-dinate difference, r ( q , t ′ ) − r ( q ′ , t ′′ ). It is so if the atmosphere is statisticallyhomogeneous. In the course of averaging, dependence of the coordinate dif-ference on the fluctuating force should be also taken into account. Thisdependence is given by the relation r ( q , t ′ ) − r ( q ′ , t ′′ ) = ( e z c + c q ′ )( t ′ − t ′′ ) − c q − q ′ ( t − t ′ )+ cq Z t ′′ t ′ dt ( t ′ − t ) F [ r ( q ′ , t )]+ cq Z t ′ t dt ( t ′ − t ) { F [ r ( q , t )] − F [ r ( q ′ , t )] } , (28)which follows from Eq. (9). The distance | r ( q , t ′ ) − r ( q ′ , t ′′ ) | should be ofthe order or less than the outer radius, L , of the turbulence. Taking intoaccount that c >> | c q − q ′ | , | c q ′ | , we infer that | t ′ − t ′′ | ≤ L /c . This meansthat in the right-hand side of Eq. (28) c q ′ in the first term and the thirdterm, which is proportional to ( t ′ − t ′′ ) , can be omitted. Then Eq. (28)reduces to r ( q , t ′ ) − r ( q ′ , t ′′ ) = e z c ( t ′ − t ′′ ) − c q − q ′ ( t − t ′ )+ cq Z t ′ t dt ( t ′ − t ) { F [ r ( q , t )] − F [ r ( q ′ , t )] } . (29)The last two terms in Eq. (29) describe the displacement of two photonsfrom each other because of the difference of their initial velocities. The term − c q − q ′ ( t − t ′ ) describes the divergence of two straight-line trajectories. Thelast term accounts for the different actions of the atmosphere on the particlesmoving in different spatial regions.Obtaining of the average values in Eq. (27), which depend on the wave-vectors r ( q , t ′ ) and r ( q ′ , t ′′ ), seems to be challenging because of the presenceof the fluctuating force in r ( q , t ′ ) and r ( q ′ , t ′′ ). Nevertheless the analysis10implifies if we neglect the correlations between the forces F α or F β and theforces entering r ( q , t ′ ) or r ( q ′ , t ′′ ). This simplification can be justified by thefollowing reasonings. The explicit value of the α -force is given by F α [ r ( q , t ′ )] = F α (cid:2) r − c q ( t − t ′ ) − cq o Z t ′ t dt ( t − t ′ ) F [ r ( q , t )] (cid:3) = F α (cid:2) r ⊥ − c q ⊥ ( t − t ′ ) + c e z t ′ − cq o Z t ′ t dt ( t − t ′ ) F [ r ( q , t )] (cid:3) , (30)where the relation z = ct is used.If the correlation exists, the distance | r ( q , t ) − r ( q , t ′ ) | can be estimatedby the value c ( t − t ′ ) ≤ L . In this case, the integral in Eq. (30) is propor-tional to ( L /c ) . Hence, the correlation between F α [ r ( q , t ′ )] and F [ r ( q , t )]can be neglected. This approximation implies the physical picture where thevariation of the photon momentum on the correlation length, L , is muchsmaller than q . Therefore, the averaging h F α F β i can be performed in twosteps. Firstly, we obtain h F α F β i considering the arguments of F α and F β to be fixed. After that, the averaging of the forces, entering the arguments,should be performed. For example, the term (25) is expressed as φ P P ′ = ω Z t Z t dt ′ dt ′′ Z d g ψ ( g ) g · (cid:0) p + k t ′ cq (cid:1) g · (cid:0) p ′ + k ′ t ′′ cq (cid:1) h e − i g [ r ( q ,t ′ ) − r ( q ′ ,t ′′ )] i , (31)where the first-step averaging results in appearance of the spectral density ψ ( g ). The second-step averaging is shown in (31) by the angle brackets. Tosimplify the derivation of σ , the authors of [16] represented the average ofthe exponential function in Eq. (31) as a product, h e − i g [ r ( q ,t ′ ) − r ( q ′ ,t ′′ )] i ≈ h e − i gr ( q ,t ′ ) ih e i gr ( q ′ ,t ′′ ) i , (32)neglecting the correlation of the photon displacements r ( q , t ′ ) and r ( q ′ , t ′′ ).Further analysis explains how this correlation can be accounted for.First of all, it should be noted that we can integrate Eq. (31) over t ′ − t ′′ because of the presence of the term e z c ( t ′ − t ′′ ) in r ( q , t ′ ) − r ( q ′ , t ′′ ) [see Eq.(29)]. The corresponding fast oscillating function, e i e z g c ( t ′ − t ′′ ) , appears in thelast factor of Eq. (31). Integration of this factor results in Z ∞−∞ d ( t ′ − t ′′ ) e i e z g c ( t ′ − t ′′ ) = 2 πc δ ( g z ) . (33)11he lower and the upper limits of the integration over t ′ − t ′′ are replaced by ∓∞ . This can be approved when the propagation time, t , is much greaterthan L /c . In other factors in Eq. (31), the substitution t ′′ = t ′ is used.The relation (33) means, that only the g x,y - components enter Eq. (31). Inparticular, the Fourier-transform ψ ( g ) should be considered as a function ofthe two-dimensional vector g ⊥ : ψ = ψ (cid:18)p g x + g y (cid:19) . This observation corre-sponds to the known Markov approximation [3] where it is assumed that theindex-of-refraction fluctuations are delta-function correlated in the directionof propagation. In fact, our derivation, based on the paraxial approximation,supports the validity of the Markov approach which at first sight seems tobe doubtful.Using Eqs. (29) and (33), the expression (31) is simplified to φ P P ′ = 2 πω c Z t dt ′ Z d g ψ ( g ) g · (cid:0) p + k t ′ cq (cid:1) g · (cid:0) p ′ + k ′ t ′ cq (cid:1) e i gc q − q ′ ( t − t ′ ) × h e − i g cq R tt ′ dt ( t − t ′ ) { F [ r ( q ,t )] − F [ r ( q ′ ,t )] } i , (34)where all the vectors have only the x − and y -components, and c q − q ′ = c ( q − q ′ ) /q .As we see from Eq. (34), to obtain φ P P ′ one needs to calculate theaverage value of the exponential function which is similar to the function in(23). Following the previous procedure, this average can be rewritten as (cid:28) exp (cid:26) − i g cq Z tt ′ dt ( t − t ′ ) { F [ r ( q , t )] − F [ r ( q ′ , t )] } (cid:27)(cid:29) = exp (cid:26) − πc Z tt ′ dt ( t − t ′ ) Z d g ′ ψ ( g ′ )( g · g ′ ) (cid:2) − h e − i g ′ · [ r ( q ,t ) − r ( q ′ ,t )] i (cid:3)(cid:27) . (35)Again, the same function appears in the exponent of the right-hand side ofEq. (35) after using the trajectories (28). Similar steps can be undertakenmany times. In this way, the time hierarchy, 0 < t ′ ≤ t ... ≤ t i ≤ t , isgenerated. If the photon-turbulence interaction time, t − t i , is short, thedisturbance of the trajectory is small and vanishes when t i → t . In this caseboth values, r ( q , t i ) and r ( q ′ , t i ), approach the value of r irrespective of theinitial momenta q and q ′ . Therefore we substitute the quantity12 (cid:10)(cid:0) g ′ · [ r ( q , t ) − r ( q ′ , t )] (cid:1) (cid:11) (36)12nstead of 1 − h e − i g ′ · [ r ( q ,t ) − r ( q ′ ,t )] i (37)assuming the exponent in Eq. (37) to be small. The linear in g ′ term inthe expansion of the exponential factor is ignored because of its zero-valuecontribution into the integral over g ′ in Eq. (35). Then the term (37) reducesto 12 (cid:10)(cid:0) g ′ · [ r ( q , t ) − r ( q ′ , t )] (cid:1) (cid:11) ≈ ( t − t ) (cid:0) c q − q ′ · g ′ (cid:1) + πc
30 ( t − t ) Z d g ′′ ψ ( g ′′ ) (cid:0) c q − q ′ · g ′′ (cid:1) ( g ′ · g ′′ ) . (38)To obtain Eq. (38), the approximate relation, F α [ r ( q , t ))] − F α [ r ( q ′ , t ))] ≈ c q − q ′ ( t − t ) ∂ r F α [ r + c q ( t − t )] , (39)where t ≤ t ≤ t , was used. This approximation is in the spirit of the previ-ous step, where the turbulence effect was assumed as a small perturbation.Substitution of Eq. (38) into the right-hand side of Eq. (35) and integra-tion over variables g ′ , g ′′ and t result inexp ( − . · − C n l ′ − / c c q − q ′ ( t − t ′ ) g " C n l ′ − / c ( t − t ′ )
560 ++ cos 2 θ C n l ′ − / c ( t − t ′ ) · ! , (40)where l ′ = l / π , and θ is the angle between the two-dimensional vectors g and q − q ′ .After substitution of (40) into (35), (35) into (34) and (34) into (25), wecalculate h I ( r , t ) i . Many-fold integrations over the variables q , q ′ , p , p ′ , k , k ′ , θ, and t ′ are performed mainly numerically with employing a computer cluster.In the course of integration, we have used the Tatarskii modification of therefractive index spectrum which is derived from the von Karman form (19)by setting L − = 0. The results for σ are shown in Figs. 1-3. Figs. 1-3 can be used to illustrate the importance of the correlations ofdifferent trajectories. To simplify our argumentations, we consider a coherent13aser beams, i.e., the case r = r . Two terms in the last braces of Eq.(21) contribute equally into h I ( r , t ) i . Moreover, if one sets φ P P ′ = 0 in Eq.(24), thus ignoring the correlations of photons with different initial momenta,we obtain h I ( r , t ) i = 2 h I ( r , t ) i . The scintillation index, σ , is equal tounity here. This physical picture is realized for a long-distance propagation (r /r ) =1z, kmr =1cm C =10 -13 m -2/3 Figure 1: Scintillation index of a coherent and partially coherent beams inthe atmosphere versus propagation distance z . Dashed curves correspondto the multiplicative approximation (32) for the photon correlations; solidcurves are obtained within the present paper’s approach [see Eqs. (35 -40)]. C n = 10 − m − / , r = 0 . m , l π = 10 − m , and q = 10 m − . The uppertwo curves correspond to the coherent beam.( t → ∞ ) when the oscillating factor e i gc q − q ′ ( t − t ′ ) confines the effective volumeof the integration over g and q − q ′ to zero [see Eq. (34)]. For finite valuesof t , the contribution of φ P P ′ becomes quite sizeable that is seen in Figs. 1-3where the values of σ are greater than unity. There is a positive contributionof φ P P ′ term into the last exponent in Eq. (24) when the vectors p , k and p ′ , k ′ have opposite signs and the difference | q − q ′ | is not too large. Themost favorable conditions are realized when p = − p ′ , k = − k ′ , q = q ′ . (41)14 z, kmr =1cm C =2.5x10 -14 m -2/3 (r /r ) =10.50.20.05 Figure 2: The same as in Fig. 1 but for a weaker turbulence strength: C n = 2 . × − m − / . 15n this case the sum φ P P + 2 φ P P ′ + φ P ′ P ′ is equal to zero . Eqs. (41) canbe interpreted as the ”super-correlation” conditions under which the valueof M is equal to unity and does not depend on the turbulence.The dependence of σ on the initial radius r can be explained as follows.The characteristic values of the initial momentum, ˜ q ∼ √ /r , is greater forsmall r . Hence the volume of integration over q − q ′ is also greater. At thesame time the corresponding increase of φ P P ′ occurs only for short distances, z , where time intervals t are sufficiently small and the oscillating factor in Eq.(34) is close to unity. Therefore, when r decreases, there is an increase of σ accompanied with the displacement of the region with enhanced fluctuationstowards small z . This is clearly seen in Fig. 3.In a similar way we can explain a considerable difference of σ found forthe plane-wave and spherical-wave models of radiation in Ref. [34] (Figs. 1and 2 there). It follows from the above reasonings that this effect arises dueto very different initial q -volumes in the two models.Also, the calculations of σ in the Ref. [15] should be mentioned where asimplified model of the turbulence was used (see Fig. 1 there). The resultsof Ref. [15] well correlate with ours.Comparing the results of the present paper and those, based on the ap-proximation of uncorrelated trajectories (32) (respectively, solid and dashedlines in Figs. 1 and 2), we see a more pronounced growth of σ at a moder-ate turbulence in the former case. Figures 1-3 illustrate that this holds truefor the distances of 1 − σ → r /r when z → ∞ .The phase diffuser with a short characteristic time (a high-frequency dif-fuser) does not change qualitatively the physical picture described above. Atthe same time, both approaches reveal an ability of the diffuser to suppressscintillations which is favorable for communication performances.The effect of the phase diffuser is explained as follows. The initial phaserelief, introduced by the diffuser, varies in time. The photon trajectoriesdepend on the initial state of the radiation and varies synchronously withthe diffuser state. A “slow” detector integrates the contribution of thesephotons. Although the atmosphere stays almost frozen during the integrationtime, the diffuser provides a better averaging of the propagating radiationover the refractive-index relief. Therefore, the fluctuations of the detectedsignal decrease. 16 =1 cm (r /r ) =0.2 z, kmC =10 -13 m -2/3 (r /r ) =1r =1 cm3 cm Figure 3: Scintillation index versus propagation distance z for different initialradii of the beam: r = 0 . m , 0 . m , 0 . m . The rest of the parametersare the same as in Fig 1.This is not a unique way to suppress fluctuations. For example, theauthors of Ref. [35] proposed to use asymmetric optical vortices. The rangeof a weak and moderate turbulence was studied. Numerical simulations ofthe beam propagation showed promising results. It should be emphasizedthat in this case the experimental setup does not require a high-frequencyphase diffuser. Our analysis is based on Eq. (21) obtained within the concept of photontrajectories. To consider photons as particles, whose density in the ( r , q )17omain is defined by the distribution function f ( r , q , t ), the uncertainty ofthe momentum, q , should be small. The value of the uncertainty can beestimated from the definition of the distribution function (1) as ˜ k/
2. Itfollows from Eq. (14) that close to the source and in the absence of thediffuser the ratio ˜ q ˜ k/ ∼ h q i / h k / i / = (2 /r ) / (2 /r ) / = 1. Hence in the vicinity of thesource, our calculations of σ are not applicable if the light is in a coherentstate.The situation changes drastically for a remote detector. With increase ofthe propagation path, z, the value of ˜ q increases. The corresponding gain ofthe photon momentum, ∆ q , is generated by a random force, F . Hence theaverage value, h ∆ q i , is equal to zero while the nonzero mean-square value isgiven by [16] h ∆ q i = 0 . π Γ(1 / q l ′ − / C n z. (42)In contrast to ˜ q , the value of ˜ k decreases because of the broadening of thebeam. The mean-square of the beam radius is given by [6, 16] R = r (cid:2) z q r r + 8 z Tr (cid:3) , (43)where T = 0 . l − / C n . When the last term in square brackets dominates,the ratio ˜ q / (˜ k/ can be estimated as h ∆ q i R ≈ · q l − / C n z , (44)where h ∆ q i is assumed to be of the order of ˜ q thus ignoring the square ofthe initial momentum 2 /r .Substituting z = 10 m , q = 10 m − , and l = 2 π · − m into Eq. (44),we obtain ˜ q ˜ k/ ∼
21 that provides adequacy of our approach for the wholerange of z variations shown in Figs. 1-3. This range concerns not onlycoherent, but also partially coherent beams. For partially coherent beams,the minimum z can be even smaller than for coherent beams. This is becauseof an additional diffuser-caused growth of ˜ q , which is estimated by the value∆˜ q diffuser ∼ /r . The paper continues the studies presented in Refs. [8, 16]. Using the ap-proach of the distribution function, the problem of obtaining of σ reduces to18alculation of the correlations between different photon trajectories. Assum-ing the outer radius of turbulent eddies much smaller than the propagationdistance, the iterative procedure for calculations of these correlations is de-veloped. The modified approach makes it possible to extend applicability ofthe theory to a wider range of the propagation distances. This range includesa strong turbulence as well as a considerable part of a moderate turbulencewhere the scintillation index tends to reach its maximum value. The cri-terium, derived in Sec. 5, imposes the restriction on our theory from the sideof short distances (weak turbulences). The authors thank V. Bondarenko, V. Gorshkov and A. Semenov for usefuldiscussions and comments.
References [1] A.P. Schuster, Astrophys. J., (1905).[2] P.W. Milonni, J.H. Carter, J.C. Peterson, and R.J. Hughes, J. Opt. B:Quantum Semiclass. Opt. , S742 (2004).[3] V.I. Tatarskii, The effect of the Turbulent Atmosphere on Wave Prop-agation. Springfield, VA: National Technical Information Service, U.S.Department of Commerce, (1971).[4] L.C. Andrews and R.L. Phillips, Laser Beam Propagation Through Ran-dom Media. Bellingham, WA: SPIE Press (1998).[5] L.C. Andrews, R.L. Phillips, and C.Y. Hopen, Laser Beam Scintillationwith Applications. Bellingham, WA: SPIE Press (2001).[6] R.L. Fante, Proc. IEEE, (1975).[7] X. Liu, F. Wang,and Y. Cai, Opt. Lett., , 3336, (2014).[8] G.P. Berman, A.A. Chumak, and V.N. Gorshkov, Phys. Rev.E ,056606 (2007). 199] M.E. Gracheva and A.S. Gurvich, Izv. Vysch. Uchebn. Zaved. Radiofiz., , 717 (1965); M.E. Gracheva, A.S. Gurvich, and M.A. Kallistratova,Radiophys. Quantum Electron, , 40 (1970).[10] R.L. Fante, Proc. IEEE, , 1424 (1980).[11] S. Wang, M. Plonus, and C. Ouyang, Appl. Optics, , 1133 (1979).[12] R.L. Fante, IEEE Tranc. Antennas Propagat., AP-25 , 266 (1977).[13] M. Lee, J. Holmes, and J. Kerr, J. Opt. Soc. Amer., , 1279 (1977).[14] V.A. Banakh and V.M. Buldakov, Opt. Spectrosk., , 707 (1983).[15] V.A. Banakh, V.M. Buldakov, and V.L. Mironov, Opt. Spectrosk., ,1054 (1983).[16] G.P. Berman and A.A. Chumak, Phys. Rev. A, , 013805 (2006).[17] O. Korotkova, L.C. Andrews, and R.L Phillips, Proceedings of SPIE, , 98 (2002).[18] O. Korotkova, L.C. Andrews, and R.L Phillips, Opt.Eng., , 330(2004).[19] F. Wang, X. Liu, and Y. Cai, Prog. in Electromagn. Res., , 123(2015).[20] R. Dashen, J. Math. Phys., , 894 (1979).[21] I.G. Yakushkin, Radiophys. Quantum Electron., , 270 (1976).[22] A. Fedrizzi, R. Ursin, T. Herbst, M. Nespoli1, R. Prevedel, Th. Scheidl,F. Tiefenbacher., Th. Jennewein, and A. Zeilinger, Nature Phys. , 389(2009).[23] R. Hughes, J. Nordholt, D Derkacs and Ch. Peterson, New J. Phys. ,43.1 (2002).[24] J. Bourgoin, E. Meyer-Scott, B. L. Higgins, B. Helou, C. Erven, H.Hubel, B. Kumar, D. Hudson,I. D’Souza, R. Girard, R. Laflamme, andT. Jennewein, New J.Phys, , 023006 (2013).2025] A. Semenov and W. Vogel, Phys Rev. A , 023835 (2010); D. Vasylyev,A. Semenov, and W. Vogel, Phys Rev. Lett., , 220501 (2012); V.Usenko et al., New J. Phys, , 093048 (2012).[26] A.I. Rarenko, A.A. Tarasenko, and A.A. Chumak, Ukr. J. Phys., ,1577 (1992); O. Chumak and N. Sushkova, Ukr. J. Phys., , 30 (2012).[27] G.P. Berman and A.A. Chumak, Proc. of SPIE, (2007).[28] G.P. Berman and A.A. Chumak, Phys. Rev. A, , 063848 (2009).[29] O.O. Chumak and E.V. Stolyarov, Phys. Rev. A, , 013855 (2013).[30] O.O. Chumak and E.V. Stolyarov, Phys. Rev. A, , 063832 (2014).[31] A. A. Tarasenko and A. A. Chumak, JETP , 625 (1977)[32] J. Strohbehn and S. Clifford, IEEE Trans. Antennas Propag., AP-15 ,416 (1967).[33] Yu.A. Kravtsov, Rep. Prog. Phys.,bf 55, 39(1992).[34] L.C. Andrews and R.L. Phillips, SPIE, , 90 (1999).[35] G.P. Berman, V.N. Gorshkov, and S.V. Torous, J. Phys. B: At. Mol.Opt. Phys.44