The interaction of ultraintense short laser beams with plasmas has attracted a great deal of attention for fundamental research and technological applications, such as particle acceleration^{[1, 2]} and terahertz (THz) radiation.^[3] Especially, there have been a lot of activities in the field of pair (electron– positron) plasma, ^[4] which usually behaves as a fully ionized gas consisting of electrons and positrons. Electron– positron plasmas have received significant attention due to its importance not only in plasma astrophysics (e.g., in the solar atmosphere, ^[5] polar regions of neutron stars, ^[6] in the active galactic nuclei, ^[7] in magnetospheres of pulsars, ^{[8, 9]} in the early universe), but also in intense laser fields, ^[10] and laboratory experiments.^{[11, 12]} The pair plasmas produced by the interaction between intense lasers with solid targets have been archived in laboratory.^{[13– 15]} Many of the astrophysical plasmas contain ions, in addition to electrons and positrons. In fact, the propagation of intense short laser beams in plasmas can also lead to pair production, resulting in a three-component electron– positron– ion (EPI) plasma.^{[10, 16– 18]} The annihilation, which takes place in the interaction of matter (electrons) and anti-matter (positrons), usually occurs at much longer characteristic time scales compared with the time in which the collective interaction between the charged particles takes place.^[19] Thus, the study of the dynamics of the nonlinear wave motions in an EPI plasma is important to understand the behavior of both astrophysical and laboratory plasmas.^{[10, 20– 22]}

There are some investigations which have been devoted to the effect of the plasma inhomogeneity on the interaction between laser pulse and plasmas.^{[16, 23– 27]} Plasma inhomogeneity can lead to pancake-like pulse acceleration in a nonuniform cold-electron– positron plasma.^[23] Periodic solutions have been found in inhomogeneous EPI plasma when the positron density is above a critical value, which strongly depends on both the inhomogeneity nature and the phase velocity of the wave.^[24] It is found that the presence of positrons and inhomogeneity results in strong modulational and filamentational instabilities, which induce strong nonlinear interactions between the laser beam and the inhomogeneous EPI plasma.^[16]

Wave self-modulation and soliton formation are two of the most remarkable features of the overall plasma dynamics. It is found that if the strength of the electromagnetic field amplitude is below the wave breaking limit in pure pair plasmas, the beam can enter the self-trapped regime, resulting in the formation of stable, self-guided 2D solitonic structures.^[28] The beam trapping owes its origin to the thermal pressure.^[28] The dynamics of the electromagnetic pulses is described by a nonlinear Schrö dinger equation (NLSE).^[29] To describe analytically the dynamics of the localized solutions of the NLSE, three methods developed are appropriate: the moment theory, ^[30] the paraxial approach, ^{[31, 32]} and the variational approach.^{[33, 34]}

In this paper, taking into account the effects of the background density inhomogeneity and ion concentration, we study self-trapping of laser beams in inhomogeneous EPI plasmas. Three kinds of inhomogeneities along the axial direction are discussed. The complex dynamics of a beam governed by the NLSE is analyzed by the variational approach. We conclude that the beam can be trapped in the potential well when its parameters are suitable. The detailed dynamics of the arbitrary field distribution are studied by simulations of the NLSE. The plasma inhomogeneity has great influence on the dynamics of laser beam propagating in it.

We investigate the propagation of a circularly polarized electromagnetic with a frequency ω ₀ and wave number k(z) along the axial direction in a smooth inhomogeneous EPI plasma. The vector potential can be written as

The system’ s evolution in terms of the vector (A) and the scalar (φ ) potentials is governed by E = − c^{− 1}∂ A/∂ t – ∇ φ , B = ∇ × A. We start with Maxwell equations complemented by the electron and positron fluid equations^{[20, 35]}

Here n⁻ (n⁺ ) is the electron (positron) number density, υ ⁻ (υ ⁺ ) is the electron (positron) velocity in the electromagnetic fields, p⁻ (p⁺ ) is the electron (positron) momentum, and n_i0(z) is the unperturbed number density of the ions. We assume that the ion density is homogenous in the radial r direction, but inhomogeneous in the axial z direction. The local equilibrium state of the three-component system is characterized by the dimensionless local charge neutrality

The coefficient denotes the ratio of the ion density to the electron density at z = 0. The ions do not respond to the dynamics under consideration and provide a neutralizing background due to their relatively large inertia. Equations  (1)– (6) are dimensionless with the following scalings: , , υ ^± ∼ υ ^± /c, p^± ∼ p^± /m₀c, n_i0(z) ∼ n_i0(z)/n_i0(0), ,  ∼  e/m₀c², where  ≡ [A; φ ], is the electron plasma frequency, e is the magnitude of the electron charge, and m₀ is the rest mass of electrons.

Since the dimensionless quiver velocity is given by υ ^± = p^± /γ ^± , where γ ^± = (1+ (p^± )²)^1/2, the transverse momentum (Eq.  (4)) is satisfied by . Hence, the electron and positron velocities are , and . The longitudinal motion of plasma is determined by

Using the assumption (the plasma density is homogenous in radical direction) and to eliminate , the continuity is rewritten as

We introduce new variables, τ ̃ = t – z/v_g(z), r̃ = r, and z̃ = z, where v_g(z) is the group velocity of the electromagnetic wave packet. The wave number of the electromagnetic wave beam satisfies the plasma dispersion relation, , where ω _p is the positron plasma frequency.^{[21, 36]} The non-dimensional dispersion relation implies v_g(z) ≤ 1 for a plasma, for which ω ₀ ≫ 1.^{[20, 36]} Therein, we assume the unperturbed electrons and positrons obey the same distribution in the axial direction, viz., . In terms of A and φ , equations  (7) and (8) are integrated to obtain , and n^± = 0.5 + 0.5(1+ | A| ²)/(1 ∓ φ )². Using Poisson equation and the equation for conservation of charge, then taking algebraic manipulations, wave equation  (2) becomes^[37]

Inserting Eq.  (1) and expressions for and n^± into Eq.  (9) and using the variables mentioned above, we obtain

where terms in ∂ ²/(∂ z̃ ∂ t̃ ) and ∂ ²/(∂ z̃ ²) are neglected, since the amplitude A in terms of z̃ is assumed to be slower than the transverse one (r̃ ).^[25] The term in , since v_g ≤ 1. By substituting the expression for n^± into Eq.  (3) and assuming that the (normalized) spatial extension of the pulse is large, i.e., L ≫ (1 + | A| ²)^1/2, ^[21] equation  (3) implies

Substituting Eq.  (11) into Eq.  (10) and neglecting terms of φ ³ and higher orders, we obtain

We linearize Eq.  (12) and obtain a modified dispersion relation

Since we consider a weak inhomogeneous EPI plasma, k(z) can be considered as a constant. Using Eq.  (13) and making the self-evident renormalization of variables ξ = z̃ (α ′ /2k), , (α ′ = α ²/(4(2 − α ))) equation  (12) can be rewritten as

where Ψ = [1 − (1 + | A| ²)^{− 2}][2n_i0(ξ ) − 1]². The function Ψ (| A| ², ξ ) represents the generalized nonlinearity.

The Lagrangian density corresponding to Eq.  (14) is

where . An appropriate variation of Lagrangian δ L/δ A* = 0 yields Eq.  (14) as the Euler– Lagrange equation. We also find that the Hamiltonian conserves in Eq.  (14), and the power . In the optimization procedure, the first variation of the variational function must vanish on a suitably chosen trial function. We will use the Gaussian-shaped beam as a trial function,

Such a choice greatly simplifies computation. The amplitude (Λ ), the beam radius (a), the wave front curvature (b), and the phase (ϕ ), as the unknown functions of the propagation coordinate ξ , will be further used in order to make the variational functional an extremum. In order to make the time-dependant problem tractable, an averaging over radial coordinate for the Lagrangian is carried out in terms of these parameters of the trial function (16)

and , where p = ρ /a. The expression for the function K(Λ ², ξ ) is

The plasma inhomogeneity is expressed by the term n_i0(ξ ). When n_i0(ξ ) = 1, the system is reduced to the homogeneous case. Demanding that the variation of the spatially averaged Lagrangian with respect to each of these parameters is zero, we obtain the corresponding set of Euler– Lagrange equations,

Equation  (19) implies that during the electromagnetic beam evolution, its power is conversed, , where Λ ₀ and a₀ are, respectively, the initial amplitude and the initial “ radius” of the EM beam at ξ = 0. Using Eq.  (19), integrating Eq.  (20) leads to

with

Equation  (23) plays the role of an effective potential for the evolution of the radius a. It has been assumed that the initial beam has a plane front (or zero curvature) (da/dξ | _{ξ = 0} = 0 = b(0)). Choosing the initial beam radius a₀ to be equal to the equilibrium radius a_e, a stationary solution of Eq.  (20) is obtained if ∂ V/∂ a| _{a= ae} = 0 and the equilibrium radius of the beam is

The expression can help us to obtain equilibrium parameters: equilibrium radius a_e and the corresponding initial power.

We conclude that the “ effective particle” (i.e., the beam) can be trapped in a potential well depending on Eqs.  (18) and (23). If the initial radius is appropriate, the beam will be trapped in self-guiding regime of propagation. We can find V(a, ξ )| _{a→ 0} → ∞ in Eq.  (23), and conclude that the beam radius is bounded from bottom in terms of the saturating nonlinearly, so the so-called collapse does not occur.

We discuss the effective potential in different inhomogeneous cases. In order to investigate the effects of ion density distribution, three types of inhomogeneous cases are discussed: a decreasing linear axial density distribution n_i0(ξ ) = 1 + bξ , b < 0, where b is a characteristic inhomogeneity parameter (Fig.  1); an increasing linear axial density distribution n_i0(ξ ) = 1 + bξ (b > 0) (Fig.  2); and a Gaussian density distribution (Fig.  3). The potential of V(a, ξ ) versus a and ξ in inhomogeneous case with b = − 0.001 is shown in Fig.  1(a). The depth of the potential well decreases with the axial ξ , due to the density variation in the axial direction. If ξ is large enough, the potential well cannot exist. The potential versus beam radius a for b = 0, − 0.001, and − 0.005 at fixed position ξ = 10 are illustrated in Fig.  1(b). We can find that if the density decreases quickly (b = − 0.005), the potential well becomes very shallow at the position.

Fig.  1.

	Figure Option ViewDownloadNew Window
	Fig.  1. (a) Potential profile of V(a, ξ ) versus a and ξ in inhomogeneous case for ni0(ξ ) = 1 + bξ , b = − 0.001, and N = 11. (b) Potential profile of V(a, ξ ) versus a for different characteristic inhomogeneity parameters b = 0, − 0.001, − 0.005 at fixed position ξ = 10.

Fig.  2.

	Figure Option ViewDownloadNew Window
	Fig.  2. (a) Potential profile of V(a, ξ ) versus a and ξ in inhomogeneous case for ni0(ξ ) = 1 + bξ , b = 0.001, and N = 11. (b) Potential profile of V(a, ξ ) versus a for different characteristic inhomogeneity parameters b = 0, 0.001, 0.005 at fixed position ξ = 10.

Fig.  3.

	Figure Option ViewDownloadNew Window
	Fig.  3. (a) Potential profile of V(a, ξ ) versus a and ξ in inhomogeneous case for ni0(ξ ) = exp[− (ξ − d)2/40000, d = 100, and N = 11. (b) Potential profile of V(a, ξ ) versus a for different characteristic inhomogeneity parameters d = 0, 50, 100 at fixed position ξ = 10.

When b > 0, the plasma density increases along with the axial direction. We plot the potential profile versus radius a and position ξ for b = 0.001 and N = 11 in this case in Fig.  2(a). We find that the depth of the potential well increases with increasing ξ , since the density increases with ξ . In Fig.  2(b), when b = 0.005 the depth of the potential well is almost twice the depth of that with b = 0.001. At the high density area, the depth of potential increases.

The potential well with a Gaussian density distribution n_i0(ξ ) = exp[− (ξ − d)²/40000] is shown in Fig.  3, where d is another different characteristic inhomogeneity parameter, and it means that the density maximum of a Gaussian density distribution is located at the position ξ = d. When d = 100, we can find the deepest potential well at the position ξ = 100 (see Fig.  3(a)). It can be predicted that the strongest self-trapping of beam can take place at the positron ξ = 100. Therefore, the potential profile has different features at fixed position ξ = 10 due to different values of d. Figures  1– 3 indicate that the inhomogeneous character has a strong influence on the potential depth and width, and will influence the beam behavior in the plasma.

In the above section, we applied a variational approach involving a Gaussian trial function. The main limitation of this approach, like the paraxial momentum, is that it is valid only in the aberration-less approximation. However, it can provide guideline for simulation. In order to account for the structural changes in the beam shape, we solve Eq.  (14) numerically for the beam dynamics of arbitrary field distribution both in homogeneous case and inhomogeneous cases. Corresponding laser parameters are:^[26] FWHM duration is 320  fs, λ = 0.526  μ m, ω ₀ = 3.6 × 10¹⁵  rad· s^{− 1}, and I = 7.9 × 10¹⁷  W· cm^{− 2}. The corresponding critical density . Since we make a theoretical analysis in this work, the plasma density can be changed as long as the electron density at the initial position satisfies n_e < n_c. Generally, the positron density is smaller than the electron density.

We take the initial profile of the beam to be Gaussian . In homogeneous case when the initial profile of the beam (solid lines in Figs.  4– 6) is close to the stable equilibrium one, the beam quickly attains the profile of the ground-state soliton and propagates for a long distance without distortion of its shape.^[38]

Fig.  4.

	Figure Option ViewDownloadNew Window
	Fig.  4. The ξ propagation of the beam fields | A(ρ = 0, ξ )| with the same initial power N = 11 and the same amplitude Λ 0 = 0.40, but with different characteristic inhomogeneity parameters (b = 0, − 0.001, − 0.005) in inhomogeneous case ni0(ξ ) = 1 + bξ .

Fig.  5.

	Figure Option ViewDownloadNew Window
	Fig.  5. The ξ propagation of the beam fields | A(ρ = 0, ξ )| with the same initial power N = 11 and the same amplitude Λ 0 = 0.40, but with different characteristic inhomogeneity parameters b = 0, 0.001, 0.005 in inhomogeneous case ni0(ξ ) = 1 + bξ .

Fig.  6.

	Figure Option ViewDownloadNew Window
	Fig.  6. The ξ propagation of the beam fields | A(ρ = 0, ξ )| with the same initial power N = 11 and the same amplitude Λ 0 = 0.40, in homogeneous case ni0(ξ ) = 1 and inhomogeneous case ni0(ξ ) = exp[− (ξ − d)2/40000] for d = 0, 100.

Figures  4– 6 show the numerical simulations in inhomogeneous cases. In all these cases the parameters of the beam (i.e., Λ ₀ and a₀) are chosen to ensure the same energy and the same initial amplitude Λ ₀. For a decreasing linear axial density distribution n_i0(ξ ) = 1 + bξ (b < 0), the background density decreases with the lager axial ξ . Compared with the oscillation of the field amplitude in homogeneous case (b = 0, see solid line in Fig.  4), the beam amplitude decreases and the beam radius increases with the axial ξ for b = − 0.001 (see dashed line in Fig.  4). When density decreases quickly (b = − 0.005), the potential well almost does not exist and the trapped light finally defocuses. Figure  4 shows that the beam defocusing can be seen in low-density plasma, in agreement with the prediction of variational approach. If b > 0 (see Fig.  5), the plasma density increases with the axial ξ . The red dashed line in Fig.  5 shows that the beam amplitude increases and radius decreases with the axial ξ . For b = 0.005, we obtain a larger beam amplitude increment than those with b = 0.001. It is found that the trapped light focuses in the high-density area. We can conclude that higher plasma density could strength the self-trapping processes. For the Gaussian types of density distribution n_i0(ξ ) = exp[− (ξ − d)²/40000], the beam field submits nonlinearly oscillating regime in Fig.  6. The beam amplitude reaches a maximum value at the highest density position ξ = d. We find that the behavior of beam can be controlled in an inhomogeneous plasma.

To summarize, the self-trapping of laser beams in inhomogeneous EPI plasmas is investigated in the paper. Starting from the Maxwell equations and electron– positron fluid equations, we derive two coupled equations describing the dynamics of nonlinearly interacting arbitrary large amplitude electromagnetic pules and plasma slow motions. Then we show that the system of the governing equations can be reduced to an NLSE with density inhomogeneity effect. The complex dynamics of a beam governed by the NLSE is analyzed by variational approach. We conclude that the beam can be trapped in the potential well when its parameters are suitable. In the inhomogeneous plasma, the depth of potential well trapping the laser beam increases with the plasma density. The detailed dynamics of arbitrary field distribution are studied by simulations of the NLSE. The beam around the equilibrium state propagates for a long distance without distortion of its shape in homogeneous case. However, the beam defocusing and self-focusing in inhomogeneous plasma can take place. The plasma inhomogeneity has great influence on the laser beam dynamics. The present results should be useful in understanding the dynamics of intense laser pulses in plasma.