Open Access
22 September 2018 Laser space debris cleaning: elimination of detrimental self-focusing effects
Alexander M. Rubenchik, Irina A. Vaseva, Mikhail P. Fedoruk, Sergei K. Turitsyn
Author Affiliations +
A ground-based laser system for space debris cleaning requires pulse power well above the critical power for self-focusing in the atmosphere. Self-focusing results in beam quality degradation and is detrimental for the system operation. We demonstrate that, for the relevant laser parameters, when the thickness of the atmosphere is much less than the focusing length (that is, of the orbit scale), the beam transit through the atmosphere produces the phase distortion only. The model thus developed is in very good agreement with numerical modeling. This implies that, by using phase mask or adaptive optics, it may be possible to eliminate almost completely the impact of self-focusing effects in the atmosphere on the laser beam propagation.



The proliferation of satellites in Earth orbit, increasing in both number and economic value, makes the problem of collision with orbital debris very pressing. A promising solution for this problem is debris removal with the help of a ground-based pulsed laser. In this approach, laser pulses ablate debris material, changing the debris velocity and moving the debris to a lower orbit where natural burn-up takes place.1,2 One problem with energy delivery to the orbit is self-focusing in the atmosphere. The nonlinear effects typically are not strong enough to destroy the laser beam in the atmosphere, but the acquired phase distortion is sufficient to filament the beam during propagation to debris (see Fig. 1), thus degrading the efficiency of the debris cleaning.3

Fig. 1

The normalized intensity distribution I(r,z)/I(0,0) {with [I(r,z)=|A(r,z)|2]}. (a) P/Pcr=50, (b) P/Pcr=1000, (c) P/Pcr=2500, and (d) P/Pcr=5000. These conditions correspond to a focusing distance L=1000  km for linear propagation.


We start the analysis following4 with the basic requirements for the laser pulse on the target. Then, we discuss beam propagation and focusing, completing the requirements for the laser pulse. We will see that in typical situations the pulse power must be hundreds of times larger than the critical power for self-focusing. Then we consider the thin window (TW) model5,6 for a description of nonlinear propagation through the atmosphere. Within this model, the propagation through the atmosphere produces only phase distortions. Our numerical modeling6 demonstrates that the TW model describes surprisingly good beam propagation through the atmosphere even in the case of beam filamentation. The applicability of the TW model yields simple expressions for the average beam description on the target: the Strehl ratio and the beam quality M2.

If propagation through the atmosphere produces only phase distortion, we suggest that preimposing the phase correction on the beam with the help of the phase mask will compensate for the detrimental effect of the atmosphere. Our modeling demonstrates that a complete compensation of the detrimental nonlinear effects is possible.

The phase correction to compensate nonlinear effects depends on the density profile in the atmosphere, which is not well defined, the laser intensity, which can fluctuate, etc. As a result, it is important to present the detailed analysis of the compensation sensitivity. During the long, free propagation to the orbit even small phase disturbances can greatly modify the beam. The study of this sensitivity is the main subject of the paper.

The big problem with sensitivity studies is to find a good figure of merit (FOM). We selected Strehl ratio and the beam quality as FOM. These are the parameters characterized laser intensity and spot size on the target. Also, in Ref. 6, we derived the simple relation that calculated these quantities in terms of laser field existing in the atmosphere, solving the laser propagation problem to the orbit. The FOM approach provides the quantative measure of the sensitivity to phase self-modulation, change of atmospheric density profile, laser altitude, etc.

We conclude that the suggested phase mask method is very robust and that its implementation is an attractive way to increase the laser system performance.


Interaction of Laser Radiation with Debris

High intensity pulsed laser radiation incident on debris vaporizes surface material, creating recoil momentum that changes the debris velocity. It is clear that an optimal laser intensity exists for any specified pulse duration. At low intensity, surface temperature and evaporation are low and recoil momentum is small. At high intensity, a large fraction of laser energy is used to create plasma, which contributes little to momentum change of the debris. A crucial parameter for pulsed laser debris removal is the coupling coefficient Cm, the ratio of momentum imparted to the target, to the incident laser energy, Cm=ΔP/E. A review of data of Cm dependence on intensity for different materials is presented in Refs. 1 and 7. The experimental data from different groups demonstrate that for broad ranges of wavelength, pulse duration, and pulse energy, the coupling coefficient maximum is reached at the intensity

where τ(ns) is the pulse duration in nanoseconds. This numerical coefficient is valid for Al alloys but does not change much for other materials and wavelengths. As a function of laser intensity, Cm is peaked not far from the vaporization threshold where plasma starts to be generated and absorptivity increases rapidly. This explains the weak sensitivity to target material. Near the maximum, Cm is not sensitive to the intensity. Typical values of Cm are 1 to 10 dyne/W.7 The latest experimental results and data discussion can be found in Ref. 8.

The fluence corresponding to optimal coupling is given as

Eq. (1)

F=2.5  J/cm2τ(ns).

Let us derive the requirement for laser pulse energy that corresponds to delivering the optimal fluence to debris targets. The energy E delivered by the laser to the target is required to be E=Fπr2, where r is the radius of the beam in the focal plane. For beam radius that accounts for beam quality and beam diffraction,

Eq. (2)

where M2 is a factor describing the beam quality in comparison with an ideal Gaussian beam, λ is the laser wavelength, L is the path length from the beam director to the target, and D is the diameter of the beam director.

The required laser pulse energy E for delivering the pulse fluence for optimal coupling is found by combining Eqs. (1) and (2) giving4


We now consider a specific example, in which λ=1μm, L=1000  km, D=3  m, and M2=2, a value of M2 that can be achieved for high-energy lasers by using spatial filters and adaptive-optic systems. For this case, r34  cm, the required pulse energy is E=11τ(ns)kJ, and pulse power P=1.1×1013/τ(ns). For the optimal point of view from laser technology,4 pulse duration 3  ns the required power P=6  TW. This is much higher than the critical power for self-focusing in atmosphere, Pcr=4  GW, which corresponds to P/Pcr=1500. The required power P scales as PM4L2/D2 and the ratio of P/Pcr for practical values of parameters is above 1000.


Laser Beam Propagation Modeling

Specific modeling is needed to describe the long distance propagation between the atmosphere and the orbit, which necessitates high accuracy of nonlinear effects modeling. We will use a code applied and validated for modeling of this problem.3

The propagation of the laser beam is described by the nonlinear Schrodinger equation,2,7 i.e.,

where Δ is the two-dimensional Laplacian operator. The analysis in Ref. 4 demonstrates that the optimal pulse length, based on physics and engineering considerations, is on the order of a few nanoseconds. For this order of pulse length, temporal dispersion can be neglected in the main order of the considered effects. The inhomogeneity of the density must be taken into account in the nonlinearity only.2,9 Here, we consider a laser beam propagating vertically (relative to the ground). This is not very different from the optimal angle for the interaction with debris, which is 30  deg from the vertical.4 The assumption of perpendicular propagation is not critically important, but it simplifies the presentation.

The resulting normalized equation in dimensionless variables has the form6

Eq. (3)

where h=Z0/LD. For R0=1  m and the parameters given above, we have LD=11,855  km, P0=0.339  GW, and Pcr=4πP0=4.258  GW for a Gaussian input beam.

It is possible to show that the fastest growing perturbations resulting in filamentation are axisymmetric,10 and that filamentation, at least initially, breaks the beam into ring-like structures. The formation and breakup of the ring structure are the well-documented pathway to Gaussian beam filamentation,11 and beam propagation can be described in the main approximation within the axisymmetric version of Eq. (3).

Let us consider the propagation of an initially Gaussian laser beam. On the surface at z=0

Eq. (4)


Here, Pin is the normalized input power of the laser beam, and the dimensionless parameter C=k0·R02/F=LD/2F is the initial beam prefocusing parameter. I0=I(0,0)=Pin[W]/(πP0[W])=4Pin/Pcr. F has the meaning of a focal distance that in this case is the debris height L. Therefore, dimensional initial prefocusing is given by C=LD/2  L.

We would like to stress that the problem under consideration, though similar in terms of the basic equations to numerous self-focusing studies (see, e.g., Refs. 3, 12 and references therein), is rather different in terms of underlying physics. We consider here light propagation over a finite distance (the thickness of the nonlinear layer), with the focusing point located beyond this region, where the propagation is linear. In this case, the self-focusing effect compresses the beam, but without the catastrophic collapse of all the energy into a small volume. Yet even the small phase disturbances after the atmospheric pass during the long propagation to the target can greatly modify the beam. The TW model developed in Ref. 6 is a powerful tool to evaluate the beam structure.


Thin Window Model

The nonlinear effects are important only in a thin layer of the atmosphere (in comparison with the distance to the target); see Fig. 2.6 After passage through the atmosphere, we have free propagation wherein the effects acquired due to nonlinear interaction beam can modify the beam. The TW makes evaluation of this process quite simple.

Fig. 2

Illustration of the TW model.


For a Gaussian initial beam, the laser field after propagation through the atmosphere has the explicit description6

Eq. (5)

where ϕ=b×exp(r2), b=I0h[1exp(z1/h)], and I0=I(0,0)=Pin[W]/(πP0[W])=4Pin/Pcr. The parameter b has the meaning of the nonlinear phase shift scale and is the analog of the B integral widely used by laser designers for a laser systems with uniform optics.13

The value of z1 is approximately a few times h; the exact value must be evaluated after comparison with numerical modeling results. It is clear that there must exist an optimal value, since for small z1 we cut out the part of the atmospheric propagation, and if z1 is too large free propagation will modify the solution (the window is no longer thin). Due to the exponential dependence in Eq. (5), the optimal value is about a few atmospheric thicknesses. The detailed studies presented in Ref. 6 give the optimal value of z1 as 5h. Following results will be presented for this z1 value.

After the beam exits the atmosphere, we have free linear propagation (see Fig. 2). From Eq. (5), we see that the phase is not quadratic, and the ensuing beam propagation is not described by the simple formulae available for the focused Gaussian input. We see that the curvature of the phase corresponds to additional focusing, and the atmosphere serves as a focusing astigmatic lens. As a result, the maximal field intensity is reached before the linear focal plane. By varying the prefocusing parameter, we can partially compensate for propagation through the atmosphere.3

The free propagation with complex enough phase distribution is not easy to evaluate. But for average beam characteristics in the focal plane (the peak intensity and the average beam spot size), it is possible to get compact expressions using the fact that the field in the focal plane is proportional to the Fourier transform of the field exiting the atmosphere.6

It is convenient to use the Strehl ratio, the ratio of the intensity in the center after propagation of the field, to that produced by the linear propagation of the Gaussian beam [Eq. (5)]. The electric field in the center of the focal spot is proportional to U(r)rdr giving for S6


Recalling that U(r)=I0exp(r22+iϕ), UG(r)=I0exp(r22), we can calculate the Strehl ratio explicitly getting the analytical expression for S6

Eq. (6)


The integral estimate of the square of the spot size on the target can be calculated in a similar way, using the Parseval formula6 as follows:

Eq. (7)

πr2=r2|A(r,L)|2rdr|A(r,L)|2rdr=4  L2|Ur|2rdr|U|2rdr.

Similarly, it is convenient to calculate the beam quality parameter M2, the ratio spot size squared [Eq. (7)] to the value calculated for the Gaussian beam. For the exit from the atmosphere field distribution [Eq. (5)], M2 can be evaluated analytically, as follows:6

Eq. (8)


The above-presented expressions for S and M2 give a quantitative description of the beam degradation due to the self-phase modulations. We got this advance in beam evaluation because, different from laser system studies, we are interested in the intensity distribution in the focal plane only.

The derived equations can have multiple applications. One can estimate quantitatively the detrimental effect of self-focusing in a case without compensation, or evaluate the effect of the finite aperture.

The fact that all effects of nonlinear interaction with atmosphere produce only phase distortion suggests that if we impose the initial phase equal ϕ [Eq. (5)] with the opposite sign, the propagation through the atmosphere will compensate the nonlinear phase distortion. During the long propagation to the orbit, even a small deviation from the TW model can produce a large deviation. We proceed to evaluate in detail the applicability of the model.


Results and Discussion

Figures 3Fig. 4Fig. 5Fig. 67 show that the results obtained by using the simple TW model are in good agreement with the direct nonlinear Schrödinger equation (NLSE) simulations. Here h=0.000506, C=k0R02/L=5.93, L=1000  km, R0=1  m, λ0=1.06  μm, and Pin/Pcr varies from 100 to 7500, which corresponds to b=0.2 and b=15, respectively. Let us note that the nonuniformity is important: even for Pin/Pcr few thousands, the b values are about 10 and we do not expect the massive beam filamentation.

Fig. 3

(a, b) Comparison of Strehl ratio S and beam quality parameter M2 computed by the NLSE [Eq. (3) (black line)] and TW model [Eqs. (6), (8) (red line)]. Green points correspond to the values computed by NLSE [Eq. (3)] with laser elevation at the height z0=3  km. (c, d) Errors calculations, with error(S)=|SNLSESTW| and error(M2)=|MNLSE2MTW2|/MTW2.


Fig. 4

(a, b) Pin/Pcr=2500—conditions correspond to Fig. 1(c); (c, d) Pin/Pcr=5000—conditions correspond to Fig. 1(d). (a, c) Intensity distribution along z. (b, d) Radial distribution at the initial focal point L=1000  km, C=5.93. Black line: solution of NLSE [Eq. (3)] with initial condition [Eq. (4)]. Green line: solution of TW model [Eq. (5)]. Red line: NLSE with preimposed phase [Eq. (9)]. Blue line: propagation of Gaussian beam in linear case.


Fig. 5

(a, c) Intensity distribution along z. (b, d) Radial distribution at the initial focal point L=1000  km, δ=Δ/b. (a, b) b=2.01   (Pin/Pcr=1000)—corresponds to Fig. 1(b). (c, d) b=10.05 (Pin/Pcr=5000)—corresponds to Fig. 1(d).


Fig. 6

The normalized intensity distribution I(r,z)/I(0,0) {with [I(r,z)=|A(r,z)|2]}. (a) Propagation of Gaussian beam in linear case, (b) δ=0, (c) δ=0.1, and (d) δ=0.1. δ=Δ/b. C=5.93 and b=10.02—the conditions correspond to Fig. 1(d).


Fig. 7

(a) Strehl ratio as a function b and δ=Δ/b.Strehl ratio varies from 0.4 (blue) to 1 (red). (b) Beam quality M2 as a function b and δ=Δ/b. M2 varies from 0.997 (red) to 6.5 (blue). Red color means that the Strehl ratio or M2 is close to 1, i.e., the beam is close to Gaussian.


Let us discuss the calculations of the beam average parameters in the focal plane, i.e., Strehl ratio S and beam quality M2. For the TW model both quantities can be calculated analytically, as in Eqs. (6) and (8). The comparison of the Strehl ratio computed with the NLSE solution and the TW model is shown in Fig. 3(a). We see that the TW model is very close to the NLSE solution, reproducing even nonmonotonic S behavior. Let us stress again that the TW model is accurate for the calculation of intensity in spot center even at large b when the beam is far from Gaussian. The calculations of beam quality M2 are presented in Fig. 3(b), together with analytical expression in Eq. (8). It is important to note that this result is also valid for beams quite different from Gaussian. We see that the TW model provides an excellent description at modest b, and slightly overestimates beam quality for large b, when the beam is already completely destroyed.

Let us mention that in the TW model the atmosphere characterized only by one parameter b. Within the NLSE, we can start with different values of intensities and atmospheric heights but the same values of b. The coincidence of these runs additionally validates the TW model. In Fig. 3, we present the comparison of the propagation starting from the sea level and from 3-km elevation. It is a practically interesting situation. To minimize the effect of atmospheric turbulence, it is attractive to place the laser as high as possible.1 One can see the very good TW applicability. With increase of laser elevation, b value decreases, decreasing the self-phase modulation effects.

In Fig. 4, we present the radial intensity distribution and the peak intensity distribution along z for the TW model of Eq. (5) and the corresponding NLSE solution. Here Pin/Pcr=1500 and 5000, with the corresponding factor nonlinear phase shift b=3.02 or b=10.05, respectively.

Figure 4 shows an excellent agreement between the TW model and the full numerical simulations based on the NLSE, both for evolution of the peak intensity with distance and for the radial beam intensity distribution in the focal plane. One can see that the TW model approximates well the exact solution of the NLSE even in the situation with well-developed filamentation [see Figs. 1(c), 1(d), and 4] with the field distribution being very far from the Gaussian beam. Note that in the TW model, the solution depends on the dimensionless parameter b=I0h[1exp(z1/h)], only, which simplifies the system optimization.

The excellent results of the TW model give hope that the initial phase predistortion can be used to compensate the nonlinear phase changes. As a result, one can have an almost perfect Gaussian beam at the atmospheric exit and the detrimental effects of self-focusing can be eliminated to a great extent. The new initial condition with the corrected phase is

Eq. (9)

δ=0 means complete compensation of the nonlinear phase with the help of phase mask, δ different from zero is the single parameter describing the incomplete compensation due to the variation of intensity, laser altitude, or effective atmospheric density.

We compare the solution of the NLSE with the initial condition [Eq. (4)] and chirp C=5.93, which corresponds to linear focusing at L=1000  km; the solution of the NLSE [Eq. (3)] with a preimposed phase [Eq. (9)]; and the solution of the linear problem with the initial condition [Eq. (4)]. The result is presented in Fig. 4. We see that the initial phase modification compensates nonlinear effects and the solution of the NLSE is very close to the linear one. Note, that as expected, the solution of the NLSE with preimposed phase correction [Eq. (9)] preserves the Gaussian shape in the transversal direction.

The parameter b in the phase is not well defined; we used the simple exponential model of the atmosphere, while the real one is more complex and depends on many parameters. The phase mask is design for specific intensity which can vary from shot to shot. This uncertainty we model with the help of parameter δ.

Hence, it is very important to understand how much our results change, if the bb+δb parameter in the preimposed phase will be different from the b parameter in the acquired phase [Eq. (9)]. We did an extensive investigation of the sensitivity of the results to the value of δb. The results of the modeling are presented in Figs. 5 and 6. One can see that, even for high intensity pulse P/Pcr=5000, a variation of b of about 10% (both positive and negative) does not affect the smooth radial distribution in the focal plane and slightly reduces the peak intensity in the focal point.

The result of the incomplete compensation is a function of both δb and b. In Fig. 7, we present the Strehl ratio and beam quality M2 as a function of both parameters. We see that the suggested scheme is robust, and nonsensitive to δb, especially for smaller beam powers. The reason is that we do not need the complete phase compensation; it is sufficient to have the efficient B integral b small enough. We see for an example that even for b=10 variation of b by 10% changes Strehl number by 5% only. We see also that detrimental effects of atmosphere reduces with laser elevation (smaller b).

The required phase pattern can be preimposed with the help of a phase mask. The spatial phase distribution is independent of the intensity, which simplifies the mask production. The amplitude of the phase variation (mask thickness) is proportional to the intensity. The above studies of δb indicate that the results are not very sensitive to the phase amplitude. In practice, a few masks will be sufficient to cover the range of intensities, making the proposed method more practical. The above results can be useful for the material processing applications.

Femtosecond material processing has become more common due to the development of new, compact lasers, high precision processing, and small collateral damage. In many cases, the processing takes place in vacuum in order to eliminate laser breakdown in air. The laser light leaving the entrance window can be distorted due to self-focusing and can attain peak intensity away from the initial focal point. This focal point displacement can be important for three-dimensional memory writing or volume Bragg gratings. The phase distortion in the window distorts the beam, increasing the spot size—an effect detrimental for the precise machining.

Let us demonstrate that this situation is similar to that considered above. For fused silica, the critical power is about 1.5 MW. For typical fs processing (pulse energy 1 mJ and pulse duration 1 ps), the pulse power is about 1000Pcr and the situation is similar to that discussed above. Thus we can directly use the results from the previous sections, for example to estimate the spot size and the Strehl ratio. The use of a phase mask to compensate the detrimental windows effects also looks promising.

We have demonstrated that the nonlinear effect of self-focusing in the atmosphere for space debris cleaning can be described with good accuracy within the TW model. Within this model, the nonlinearity produces phase front distortion, serving as a high aberration focusing lens. Optical phase distortion results in displacement of the focusing point and beam filamentation, degrading the system performance.

We demonstrated that the TW model describes beam focusing with strong aberrations, with field distributions in the focal plane far from Gaussian. The pattern of laser field is determined by a single dimensionless parameter b similar to the B integral used in laser design to control the self-phase modulation. The dependence on only a single parameter greatly simplifies the optimization of the beam prefocusing arrangements.

The description of linear propagation after exit from the atmosphere can be simplified using the fact that the field in the focal plane is proportional to the Fourier transform of the field exiting atmosphere. As a result, we obtained simple expressions for the peak intensity (Strehl ratio) and beam quality M2 which can be calculated in terms of the exiting field. The use of these quantities as FOM provides a convenient quantative description of the beam structure on the target.

Because of the high accuracy of the TW model, one can compensate for self-focusing in the atmosphere by a preimposed phase distribution, which will cancel the nonlinear phase acquired during propagation in the atmosphere. Our modeling demonstrated that the detrimental effects of self-focusing even for P/Pcr5000 can be almost completely eliminated by a preimposed phase calculated within the TW model.


This work was supported by the Russian Science Foundation (Grant No. 17-72-30006). Part of the work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344



C. R. Phipps et al., “Removing orbital debris with lasers,” Adv. Space Res., 49 (9), 1283 –1300 (2012). ASRSDW 0273-1177 Google Scholar


C. R. Phipps et al., “Orion: clearing near-earth space debris using a 20-kw, 530-nm, earth-based, repetitively pulsed laser,” Laser Part. Beams, 14 (1), 1 –44 (1996). LPBEDA 0263-0346 Google Scholar


A. M. Rubenchik, M. P. Fedoruk and S. K. Turitsyn, “The effect of self-focusing on laser space-debris cleaning,” Light: Sci. Appl., 3 e159 (2014). Google Scholar


A. M. Rubenchik, A. C. Erlandson and D. Liedahl, “Laser system for space debris cleaning,” in AIP Conf. Proc., 448 –453 (2012). Google Scholar


J. H. Marburger, “Self-focusing: theory,” Prog. Quantum Electron., 4 (1), 35 –110 (1975). PQUEAH 0079-6727 Google Scholar


I. A. Vaseva et al., “Light self-focusing in the atmosphere: thin window model,” Sci. Rep., 6 30697 (2016). SRCEC3 2045-2322 Google Scholar


C. Phipps et al., “Impulse coupling to targets in vacuum by KrF, HF, and CO2 single-pulse lasers,” J. Appl. Phys., 64 1083 –1096 (1988). JAPIAU 0021-8979 Google Scholar


R. A. Lorbeer et al., “Experimental verification of high energy laser-generated impulse for remote laser control of space debris,” Sci. Rep., 8 8453 (2018). SRCEC3 2045-2322 Google Scholar


A. M. Rubenchik, M. P. Fedoruk and S. K. Turitsyn, “Laser beam self-focusing in the atmosphere,” Phys. Rev. Lett., 102 (23), 233902 –233904 (2009). PRLTAO 0031-9007 Google Scholar


V. E. Zakharov and A. M. Rubenchik, “Instability of waveguides and solitons in nonlinear media,” Sov. Phys. JETP, 38 (3), 494 –500 (1974). SPHJAR 0038-5646 Google Scholar


G. Fibich et al., “Control of multiple filamentation in air,” Opt. Lett., 29 (15), 1772 –1774 (2004). OPLEDP 0146-9592 Google Scholar


Self-focusing: Past and Present. Fundamentals and Prospects, Springer, New York (2009). Google Scholar


A. E. Seigman, Lasers, University Science Books, Mill Valley, California (1986). Google Scholar

Biographies for the authors are not available.

© 2018 Society of Photo-Optical Instrumentation Engineers (SPIE) 0091-3286/2018/$25.00 © 2018 SPIE
Alexander M. Rubenchik, Irina A. Vaseva, Mikhail P. Fedoruk, and Sergei K. Turitsyn "Laser space debris cleaning: elimination of detrimental self-focusing effects," Optical Engineering 58(1), 011003 (22 September 2018).
Received: 20 June 2018; Accepted: 9 August 2018; Published: 22 September 2018 Logo
Cited by 2 scholarly publications.
Get copyright permission  Get copyright permission on Copyright Marketplace
Atmospheric propagation

Atmospheric modeling


Laser beam propagation

Pulsed laser operation


Beam propagation method

CHORUS Article. This article was made freely available starting 22 September 2019

Back to Top