Gross-Pitaevskii-Poisson model for an ultracold plasma: Density waves and solitons

We introduce one-and two-dimensional (1D and 2D) models of a degenerate bosonic gas composed of ions carrying positive and negative charges (cations and anions), under the condition of the electroneutrality. The system may exist in the mean-ﬁeld condensate state, enabling the competition of the Coulomb coupling, contact repulsion, and kinetic energy of the particles, provided that their effective mass is reduced by means of a lattice potential. The respective model combines the Gross-Pitaevskii (GP) equations for the two-component wave function of the cations and anions, coupled to the Poisson equation for the electrostatic potential mediating the Coulomb interaction. In addition to its direct introduction, the contact interaction in the GP system can be derived, in the Thomas-Fermi approximation, from a system of three GP equations, which includes the wave function of heavy neutral (buffer) atoms. In the system with fully repulsive contact interactions, we construct stable spatially periodic patterns (density waves, akin to ionic crystals). The transition to the density wave is identiﬁed by analysis of the modulational instability of a uniformly mixed neutral state. The density-wave pattern, which represents the system’s ground state (GS), is accurately predicted by a variational approximation. In the 2D case, a stable pattern is produced too, with a quasi-1D shape. The 1D system with contact self-attraction in each component produces bright solitons of three types: neutral ones, with fully mixed components, dipoles,


I. INTRODUCTION AND THE MODEL
For Bose-Einstein condensates (BECs), created in ultracold atomic gases [1], an exceptionally accurate dynamical model is provided by the Gross-Pitaevskii (GP) equation [2].The GP equation was derived for BEC with contact interparticle interactions, and extended for the gas of dipolar atoms, with long-range interactions between them.Various aspects of theoretical and experimental studies of dipolar condensates are summarized in Refs.[3][4][5][6].Further, the analysis was developed, chiefly in the context of astrophysical and cosmological models, for BEC bound by gravity forces [7][8][9][10][11][12][13]. It was also demonstrated that similar long-range attractive interactions can be artificially induced in atomic condensates by means of specially designed laser illumination [14,15].A natural model for the BEC with gravitational or Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
pseudogravitational interactions is a system of the GP equation for the mean-field atomic wave function and Poisson equation for the gravitational potential.A model of BEC in the gas of particles carrying electric dipole moments is also based on the GP-Poisson (GPP) system, making use of two different mean-field approximations: one for the particles' wave function, and another for the interaction of the particle's dipole moment with the electrostatic potential, which is created, as per the Poisson equation, by the distribution of the polarization density in the gas [16].
The current work with ultracold plasmas [17][18][19][20][21][22] suggests a possibility to consider them in the state of quantum degeneracy [23][24][25][26][27][28].Although presently available experimental techniques make it possible to cool ions down to the lowest level of ∼10 μK [17,[29][30][31][32][33], which remains ca.two orders of magnitude above the temperature of the BEC transition, the ongoing progress in the experimental studies makes it plausible that a plasma in the BEC state will be eventually created.
The theoretical analysis of quantum plasmas was developed for electron-ion mixtures [23][24][25].They are modeled by the GPP system, including the nonlinear Schrödinger (NLS) equation for the wave function of electrons, , coupled to the electrostatic potential.In turn, the potential is created by the respective density, | | 2 , pursuant to the Poisson equation.The NLS equation for includes a self-repulsion term ∼| | 4/D , where D is the spatial dimension; thus it is similar to the effective quintic nonlinearity of the Tonks-Girardeau gas in the 1D setting [34], usual cubic nonlinearity of BEC [2] in 2D, and the nonlinear term known in the density-functional model of Fermi gases in 3D [35,36], which may be applied to experimentally available settings [37].
A less common type of plasmas represents a mixture of positive and negative ions (cations and anions), containing a negligibly small density of free electrons.Such plasmas were created with different [38,39] and equal [40] cation and anion masses.Theoretically, an ultracold plasma composed of cations and anions with equal masses was addressed in Ref. [41], by means of the molecular-dynamics methods, i.e., numerically solving coupled equations of motion for all particles in the plasma.In particular, known results [42] suggest that a candidate for achieving this purpose may be a mixture of deuterium anions and cations.
Assuming that an electroneutral cation-anion plasma may be, eventually, cooled down into the BEC state, we here introduce a model for it, which includes a two-component GP equation for mean-field wave functions ψ ± of the positive and negative bosonic ions, and the Poisson equation for the electrostatic potential, φ, induced by the distribution of the local charge density.The objective is to produce stable GPP states in the form of spatially periodic density waves and bright solitons.Chiefly, the results are reported for the effectively 1D setting, and some findings for density waves are obtained for the 2D system too.
Before introducing the model written in a scaled form, it is relevant to estimate possible physical parameters of the setting under the consideration.The state of degenerate plasmas is controlled by the Coulomb-coupling strength C, viz., ratio of the energy of the electrostatic interaction between charge carriers to their kinetic energy [43] (the latter energy term is restricted by the condition of the degeneracy of the gas).For a high-density plasma (up to 10 18 cm −3 ) of light cations and anions, an estimate yields C ∼ 10 4 (very strong Coulomb interaction is actually a reason impeding the creation of quantum plasmas in the experiment).At C 10 3 , the strongly coupled plasma is expected to build states such as Wigner crystals [44], different from those predicted by the GP theory, as the gradient terms becomes negligible in it.The plasma may be kept in the regime of moderate Coulomb coupling if the effective ion mass, m * , is made essentially smaller than its bare value, by means of a lattice (spatially periodic) potential.Lattice potentials were employed for such purposes in many theoretical and experimental studies, making use of the fact that the effective mass may be strongly reduced close to edges of the band gap in the respective spectrum [45,46].In particular, it was demonstrated theoretically [47] and experimentally [48] that one may readily reduce m * by up to two orders of magnitude, against the bare value.
As concerns the comparison of the Coulomb and kinetic energies, it is relevant to stress that the patterns reported below are actually produced by the competition of the electrostatic interaction not with the quantum dispersion (which corresponds to the kinetic energy in the GP equation), but, chiefly, with contact interactions between the ions (solutions for the density waves and solitons do not exist in the absence of the contact interactions).For the same typical values of physical parameters which produce the above-mentioned estimate, C ∼ 10 4 , the ratio of the contact-interaction and kinetic energies is estimated as being ∼10 2 -10 3 , if the bare atomic mass is used.Thus the introduction of the reduced effective mass may make magnitudes of the latter energy terms comparable, which helps to explain the possibility of the creation of the spatial patterns considered below.
According to what is said above, the GPP system is written, in the scaled form, as where m * ± and q ± are scaled effective masses and absolute values of charges of the cations and anions, g ++ , g −− , g +− are coefficients of the contact nonlinearity, and U ± (r) are trapping potentials, if any (if the lattice potential is used, as outlined above, to renormalize the effective masses, it is not explicitly included, as its effect is represented by m * ± ).Equations ( 6)-( 8) conserve two numbers of ions, where dr = dx, dx dy, or dx dy dz, in the 1D, 2D, and 3D cases, respectively, and the electroneutrality condition implies Below, we focus on the consideration of the symmetric system in free space (U ± = 0), with equal effective masses which are fixed, by rescaling the coordinates, to be m * + = m * − ≡ 1; further, the symmetry implies g ++ = g −− ≡ g, g +− ≡ G, and q + = q − ≡ q.Then, q + = q − implies N + = N − ≡ N, as per Eq. ( 5).Finally, making use of the remaining scale invariance of the GPP system, we set g = ±1 for the repulsive and attractive signs of the contact self-interaction, respectively, unless g = 0 .Thus the symmetric version of the GPP system of Eqs.(1)-(3) takes the form of In fact, the same system can be derived in the framework of a more general setting, which includes, in addition to the cations and anions, a buffer component of heavy neutral bosonic atoms with wave function 0 (the neutral atoms may be also used for cooling the ions by means of the sympathetic method [31,33]).Thus, in the limit case when the contact interactions are much weaker than the electrostatic coupling mediated by potential φ (which is a natural situation, as mentioned above), two GP equations ( 6) and ( 7) are replaced by three, the extra one being the equation for 0 : where M is the relative mass of the buffer atoms, while f and F are coefficients accounting for, respectively, the buffer-ion and self-buffer interactions.For heavy atoms with large M, the kinetic-energy term in Eq. ( 11) may be neglected, which is tantamount to the Thomas-Fermi approximation [2], applied to this equation.Then, looking for a solution as 0 = exp (−iμ 0 t )ψ 0 (x, y, z; t ) and ψ ± ≡ exp [−i( f /F )μ 0 t] ψ± , where μ 0 is the buffer's chemical potential and ψ 0 is a slowly varying function of t in comparison with exp(−iμ 0 t ), one can use Eq. ( 11) to eliminate the buffer's density, The substitution of this approximation in Eqs. ( 9) and ( 10) leads back to Eqs. ( 6) and ( 7) for wave functions ψ± , with effective coefficients Spatial patterns are expected to exist if the Coulomb coupling between cations and anions is balanced by the contact repulsion between them; therefore, we adopt G > 0 in Eqs. ( 6) and (7) [in particular, this implies the choice of F < 0, i.e., self-attraction of the buffer atoms, if the nonlinearity coefficients are produced by Eq. ( 13)].
It may be interesting to consider the system which includes, on a par with the cations and anions, polar molecules formed as their bound states, which should be represented by a separate wave function.The respective model should be based on a system of three GP equations, including various "reactions," such as merger of colliding ions into the molecule [however, the merger may be suppressed by the strong contact repulsion with G > 0 in Eqs. ( 6) and ( 7)], and breakup of the molecule due to collisions.In the present paper, we do not aim to address such a system.
Note also that Eqs. ( 6) and ( 7) with g = 0 [or with g given by Eq. ( 13), if the ion-ion interaction is mediated by the buffer atoms] may model a fermionic plasma [49], in which the Pauli principle forbids direct self-interaction; cf.Ref. [50].Furthermore, the derivation of the NLS-Poisson system for the quantum plasma dominated by the electron component [23][24][25] suggests that, in the case of g = +1 and imaginary G = −i , the 2D version of Eqs. ( 6)-( 8) may serve as a model of the degenerate electron-positron plasma, with coefficient > 0 representing annihilation losses, although the consideration of this possibility is beyond the scope of the present work.
The rest of the paper is organized as follows.Spatially periodic density-wave solutions, in both 1D and 2D forms, are considered below in Sec.II.That section includes an exact analytical investigation of the modulational instability (MI) of the uniformly mixed (locally neutral) state, and an analytical approximation for the density-wave patterns.Solutions for 1D solitons, which may exist in the form of neutral localized states, with fully mixed cations and anions, or dipole and quadrupole ones, are addressed in Sec.III.In particular, the transition from neutral solitons to dipoles is predicted analytically.The chart of the system's ground states (GSs), represented by the neutral, dipole, or quadrupole solitons, is produced in a numerical form.Collisions between moving dipole solitons are considered too, by means of direct simulations.The paper is concluded by Sec.IV.

A. Stationary equations
Stationary 1D solutions of Eqs. ( 6)-( 8) with chemical potentials μ ± of the two components are looked for as with real wave functions u ± (x) obeying the stationary version of the GPP system in free space (U ± = 0): Note that Poisson equation ( 8) can be solved by means of its Green's function, which yields (cf.Ref. [51], where the Green's function was used to solve the Poisson equation for microwave field coupling two different states of neutral atoms, represented by two components of a spinor wave function).Accordingly, Eqs. ( 15) and ( 16) may be replaced by a system of integro-differential equations: Strictly speaking, the Poisson equation should be taken in the 3D form, even if the atomic components are confined to the effectively 1D setting, as usual, with the help of a tight confining potential applied in the perpendicular plane [52,53].However, the logarithmic form of the fundamental solution of the 2D Poisson equation in the orthogonal plane implies that the transverse component of the electric field will give rise to a relatively weak force acting on the ions, and their transverse motion will be suppressed by the confining potential.

B. Modulational instability (MI) of the uniformly mixed neutral state
First, the GPP system of Eqs. ( 15)-( 17) gives rise to obvious uniformly mixed (locally neutral) states, with These states may be subject to MI against perturbations demixing the two components.The well-known condition for the instability of the mixed state against the phase separation in the absence of the Coulomb interaction is G > g [54].In the present case, one may expect that the instability threshold is shifted to larger values of G, as the Coulomb attraction between the components tends to enhance the trend to their mixing.
To investigate the MI, we follow the usual approach, substituting, in the 1D version of GP equations ( 6)-( 8), the amplitude-phase form of the wave functions, [55].Then, equations are linearized for modulational perturbations, δu ± = u ± − u 0 , χ ± , and φ.Solutions for eigenmodes of the small perturbations are looked for as where (δu (0) ± , χ (0) ± , φ (0) ) and p are amplitudes and an arbitrary wave number of the perturbation, while γ is the respective MI gain.
The substitution of perturbations (22) in the linearized GPP system leads to the dispersion equation, relating γ and p, which is written in the form of the corresponding determinant in the Appendix.Finally, a straightforward calculation yields four branches of γ (p): In this section, we consider the case of repulsive selfinteractions, with g > 0.Then, the branch of the dispersion relation given by Eq. ( 23) is purely imaginary and does not lead to MI.On the other hand, branch (24) produces real values of γ , which represent MI, when the density of the uniform neutral state exceeds a critical value (i.e., the nonlinearity is strong enough in comparison with the Coulomb attraction): At u 2 0 = (u 2 0 ) cr , the MI emerges at wave numbers with The largest value of the squared MI gain, is attained at p 2 max = 2(G − g)u 2 0 , which is tantamount to the value given by Eq. ( 26) at the MI-onset point, u 2 0 = (u 2 0 ) cr .Note that p max does not depend on the charge, q, being the  6)-( 8) with G = 2 and g = q = 1, in the domain of size L = 10 with periodic boundary conditions, starting from input (28).(b) The evolution of same as in the case of the MI against demixing of the two components in the absence of the Coulomb interactions [54].On the other hand, the critical value of the squared wave number, given by Eq. ( 26), does not depend on coefficients G and g of the contact interactions.
The predicted stability of the uniform neutral state against modulational perturbations at u 2 0 (u 2 0 ) cr is readily confirmed by direct simulations of Eqs. ( 6)-( 8) in 1D.As an example, Fig. 1(a) displays the evolution initiated by input with parameters G = 2, g = 1, q = 1, in the domain of size L = 10.In this case, u 2 0 = 10 is definitely smaller than critical value (25), (u 2 0 ) cr = 8π .For this reason, Fig. 1(a) demonstrates stable propagation of perturbations, which may be considered as ion-acoustic waves in the cation-anion plasma; cf.Ref. [56]. Figure 1(b) shows the evolution of |ψ + | at x = L/2.The numerically found temporal period of the wave is T ≈ 0.28, which is very close to that predicted by Eq. ( 24) at the same values of the parameters, 2π/|γ | = 0.283.

C. Density waves
Stable spatially periodic solutions were obtained as solutions of Eqs. ( 6)-( 8), produced by means of the imaginarytime integration method [57], in the 1D domain of size L with periodic boundary conditions.Figure 2(a) displays a typical example, with L = 10 and the density-wave's period l = L/12.
The pattern displayed in Fig. 2(a) demonstrates phase separation (demixing) of the two components, in the case when their contact repulsion, accounted for by G > 0 in Eqs.(15) and (16), is stronger than the Coulomb attraction between the cations and anions.Naturally, for fixed values of the system's parameters, viz., g, G, q in Eqs. ( 6)- (8), and fixed size L, the destabilization of the uniform state takes place, with the increase of the effective nonlinearity strength, when the density of each component exceeds a certain critical value: N ± /L ≡ N/L > n cr ; see Eq. ( 38) below.Numerical findings for the transition from the uniform neutral state to the density wave are summarized in Fig. 3(a), which displays  15), (16), and (18) (the density wave) with components u + (x) and u − (x) (solid green and dashed blue lines, respectively) in the domain of size L = 10 with periodic boundary conditions.The chemical potential and norms of the numerical solution are μ ± = 94.85 and N ± = 325.(b) The fit of the numerical solution for u + (x) (the continuous line) to the analytical ansatz, u 0 + A cos(2π x/l + α) (the dashed line), with u 0 = 5.49 and A = 2.17 determined numerically as the spatial average value of u + (x), and the square root of the average of [u + (x) − u 0 ] 2 .The best-fit value of the phase shift is α = 0.68.The system's parameters are g = 1, G = 2, and q = 1.
the modulation depth of the arising pattern, defined as the half-difference between local maxima and minima of φ + (x), [it is the same for φ − (x)] as a function of N ± ≡ N, with the same system's parameters as in Fig. 2. The striped pattern may be considered as an analog of an ionic crystal [58] (unlike the above-mentioned Wigner crystal, which is a pattern emerging in strongly coupled degenerate plasmas, including the quasi-1D setting [43,44]).The critical point, N = N cr ≡ n cr L, at which the periodic pattern appears, and the shape of the pattern can be predicted in an analytical form.To this end, we note that the energy of  29), vs the norm of each component, N ± = N, for the same system's parameter as in Fig. 2, i.e., g = 1, G = 2, q = 1, and L = 10.The dashed line shows A(N ) as predicted by analytical approximation (35) with l = L/12.(b) The chain of rhombuses and dashed line depict, severally, the numerically found chemical potential of the same family of the density waves, μ(N ), and its analytically predicted counterparts, given by Eq. (37).In addition, the chain of crosses shows μ as predicted, in a simple form, by Eqs.(36) and ( 32), (35).
the GPP system, represented by Eqs. ( 6)- (8), is the last term being the electrostatic energy.Then, the demixed pattern with spatial period l may be approximated by a simple variational ansatz, where u 0 is the mean value, and modulation depth A is the same as defined by Eq. ( 29), the respective density of particles being Figure 2(b) displays the fit of a typical numerically obtained component u + (x) to ansatz u 0 + A cos(2π x/l + α).Further, the substitution of ansatz (31) in Eq. ( 17) and consideration of the balance condition for the fundamental harmonic, cos (2π x/l ), yields the corresponding approximation for the electrostatic potential: Next, substituting components ( 31) and ( 33) of the ansatz in expression (30), and using Eq. ( 32) to eliminate u 2 0 in favor of the density of particles, we find the corresponding energy density, In the framework of the variational method [59], the modulation depth is determined by minimization of the energy density with respect to A, i.e., dE/d (A 2 ) = 0, which yields The chemical potential of the density wave can be predicted by taking Eqs. ( 15) and ( 16) at points where cos (2π x/l ) = 0; hence Eqs. ( 31) and (33) give u + = u − = u 0 and d 2 u ± /dx 2 = φ = 0: Figure 3(b) shows the numerically found chemical potential μ (rhombuses) with μ predicted by Eq. (36) (crosses) and, additionally, the chemical potential calculated, as usual, as the derivative of the energy density with respect to the total density (2n): obtained from Eq. ( 34) and depicted by the dashed line.The result given by Eq. ( 35) makes sense at A 2 0, i.e., at At the critical point, n = n cr , the MI is driven by perturbations with spatial period In fact, Eqs. ( 38) and (39) yield exact results, which are identical, respectively, to Eqs. ( 25) and ( 26), with l cr ≡ 2π/p cr .This is explained by the fact that the ansatz based on Eqs. ( 31) and (33) with infinitesimal A reproduces the exact MI eigenmode which leads to Eqs. ( 23)- (27).In particular, at q = 1 and L = 10, Eq. ( 38) yields N cr ≡ Ln cr ≈ 251.3, which precisely agrees with the numerically found critical value of N in Fig. 3(a).For the same case, q = 1, Eq. ( 39) yields l cr = √ π/2 ≈ L/11.284 for L = 10 [the calculation of the critical value (39) is not subject to the condition that ratio L/l cr must be an integer].Picking up a close integer, 12, i.e., l = L/12, the dashed line in Fig. 3(a) shows A, as predicted by Eq. ( 35), vs N for l = L/12 and fixed q = 1, G = 2, and g = 1.The figure demonstrates very accurate agreement of the analytically predicted A(N ) dependence with the numerically found counterpart.
It is relevant to compare energies E of different stationary states sharing the same value of N. The calculation of E as per Eq.(30) demonstrates that, at N > N cr , the energy is always lower for the density-wave pattern than for the uniform (neutral) state existing at the same value of N. Therefore, it is plausible that the spatially periodic pattern realizes the system's GS.
In addition to the above family of spatially periodic solutions with two free parameters, l and n, found in the numerical and approximate analytical forms, it is possible to find a family of exact analytical solutions, in the form given by Eq. ( 14), with The norms of solution (40) are Amplitude U 0 is a single free parameter of the family, while the spatial period takes the single value, l 0 .On the contrary to the above numerical solutions, the exact ones (40) are unstable, as illustrated by simulations of the evolution displayed in Fig. 4 for the exact solution with amplitude U 0 = 2 √ 5.The simulations were run in the domain of size L = 8l 0 , to make l 0 compatible with the periodic boundary conditions.The instability, which leads to establishment of an apparently turbulent state, is explained by the fact that the energy of exact solution (40) is larger than that of the uniform neutral state with the same norm: i.e., the exact solution represents an excited state of the system.Note that, for these values of the parameters, norm ( 41) is N = 80 √ 2π ≈ 200.53, being smaller than the respective critical value, N cr = 8π q 2 L ≈ 503.99, as given by Eqs. ( 25) and (38); hence the uniform neutral state with the same N is indeed stable, representing the system's GS.
Numerical solution of the 2D version of Eqs. ( 6)-( 8) readily produces stable spatially periodic patterns in the form of quasi-1D stripes, see an example in Fig. 5 6)-( 8) in a square of size 5 × 5 with periodic boundary conditions; see the definition of black stripes given by Eq. ( 43).The parameters are G = 2 and g = q = 1.
N + = N − = 650; hence the corresponding average densities are n = N ± /25 = 26.Note that the above MI analysis pertains equally well to all spatial dimensions.For the present values of the parameters (G = 2, g = q = 1), the critical density, given by Eqs. ( 25) and (38), is n cr = 8π ≈ 25.133.The existence of the stable striped pattern is natural, as n = 26 exceeds n cr , hence the uniform neutral solution should be replaced by the striped one, as the respective GS.

A. Dipole modes: Analytical and numerical results
The underlying system of Eqs. ( 6)-( 8) may give rise to bright solitons in the case of the attractive sign of the intracomponent contact nonlinearity, i.e., g < 0. It is expected that the Coulomb attraction will tend to keep the two components together, in competition with the contact repulsion between them, accounted for by G > 0. Thus creation of a dipoleshaped soliton is expected.
Figure 6 shows a set of stable solitons, obtained by means of the imaginary-time-integration method in the domain of size L = 10 with zero-flux (Neumann's) boundary conditions and fixed norms, N ± = 10.It is seen that L = 10 is sufficient to produce completely localized states.For fixed parameters g = −1, q = 0.3 and varying G, the two components remain fully overlapped at G 0.05, splitting at G > 0.05.In the state with coinciding components, an obvious exact soliton solution with the center set at x = L/2 (see Fig. 6) is provided that G < |g|.Here, x ≡ x − L/2, and A is an arbitrary amplitude of the soliton.
To approximate solutions with small separation 2ξ between the components, we adopt a variational ansatz FIG. 6. Stationary profiles of u + (x) and u − (x) (solid and dashed lines, respectively) of stable dipole-shaped solitons, produced as solutions of Eqs. ( 6)-( 8) in the domain of size L = 10 with Neumann's boundary conditions.The profiles are shown at a set of decreasing values of the contact-repulsion strength, G = 0.15, 0.125, 0.1, 0.075, 0.05, and 0.025.All solitons have norms N ± = 10.Other parameters in Eqs. ( 6)-( 8) are g = −1 and q = 0.3.
suggested by solution (49): with norms The dipole moment of the weakly split state represented by ansatz (50 (52) The respective Poisson equation ( 8) takes the form of Straightforward integration of Eq. ( 53) yields The Coulomb energy of the soliton with the slightly split components can be readily evaluated, using approximation (54): where the expansion is built for small ξ 2 .FIG. 7. Numerically found values (rhombuses) of the halfdistance between centers of the |ψ + | 2 and |ψ − | 2 components in the soliton.The dashed line shows the analytical prediction for the same separation, produced by Eq. ( 58), with A 2 expressed in terms of N ± by means of Eq. ( 51).Parameters are g = −1 and q = 0.3.The norms are N ± = 10 and the size of the domain of the numerical solution is L = 10.
Further, the energy of the contact interaction of the two components in ansatz ( 50) is found as Then, the total energy of the interaction of the two weakly separated components is Thus the value of separation ξ is predicted as one for which the total interaction energy attains a minimum, at A 2 > 5π q 2 /G and ξ 2 = 0 at A 2 > 5π q 2 /G.The dependence of ξ on G, as produced by Eq. ( 58), along with its numerically found counterpart, is shown in Fig. 7.The weak splitting of the components, with small ξ , is well predicted by the analytical approximation, while the simple approximation used above becomes irrelevant at larger values of the separation.In particular, for given G < |g|, the splitting takes place provided that where Eq. ( 51) is used to eliminate A 2 in favor of the physically relevant norm parameter.At q 2 = q 2 cr , Eq. ( 58) gives ξ = 0. Equation ( 59) may be easily inverted, to produce a critical value of G at which the splitting sets in, for given q.

B. Quadrupole solitons
The system under the consideration gives rise not only to dipole solitons, but also to quadrupoles, in which one component is located at the center, while the other one splits in two side lobes.The respective quadrupole moment is cf. Eq. ( 52).The simplest ansatz which may approximate quadrupole solitons is with small ε > 0. The norm of ansatz (61) keeps value (51) at order ε.The quadrupole moment produced by the substitution of ansatz (61) in Eq. ( 60) is QM = An example of a numerically obtained stable quadrupole soliton is shown in Fig. 8(a) for G = 0.9, g = −1, and q = 0.3.Note that ansatz (61) gives rise to a minimum of the ψ − component at x = 0 and a pair of adjacent maxima at ε > 1/2 (the quadrupole's shape is somewhat similar to that of dark solitons predicted for a quantum electron plasma in Ref. [25]).With these parameters, quadrupoles are stable in the interval of values of their norms 4.4 N ± 5.6.To identify the system's GS, Fig. 8(b) shows the energy difference between the dipole and quadrupole solitons and the neutral one, given by Eq. ( 49).The solitons with equal norms, N ± = 5.25, are compared here.In Fig. 8(a), the energy differences are shown as functions of the strength of the intercomponent repulsion, G.It is seen that neutral, dipole, and quadrupole solitons realize the GS (energy minimum) at G < 0.4, 0.4 < G < 0.7, and G > 0.7, respectively.
The results are summarized in the chart displayed in Fig. 8(c), which identifies the system's GS as neutral (N), dipole (D), or quadrupole (Q) solitons, in the parameter plane of N ± and G, with fixed g = −1 and q = 0.3.In the figure, the neutral solitons lose their stability and give rise, by splitting, to the quadrupole and dipole solitons at boundaries between the N and Q or D areas, respectively (at the boundaries, the solitons of different types have equal energies).It is natural that the increase of the repulsion strength, G, drives splitting of neutral solitons into quadrupole and dipole ones in Fig. 8(c).It is natural too that further increase of G leads to additional fragmentation, replacing the dipole solitons by quadrupole ones as the GS.
No region of coexistence of stable solitons of different types could be found.In particular, in the part of area D above the short-dashed line in Fig. 8(c), quadrupole solitons exist but are unstable, spontaneously transforming into dipole counterparts; see an illustration in Fig. 9(a).In the Q area, dipole  6)-( 8) in the domain of size L = 20, with Neumann's boundary conditions.The parameters are g = −1, G = 0.9, and q = 0.3.The quadrupole moment of the soliton [see Eq. ( 60)] is QM = −13.54.(b) The energy difference, E − E 0 , between quadrupole and dipole solitons (solid and dashed lines, respectively) and the neutral one (49), for fixed norms, N ± = 5.25, vs G, keeping g = −1 and q = 0.3 fixed too.The GS (ground state) is identified as one with the lowest energy.(c) The chart which identifies GS as different soliton species (Q, D, and N: quadrupole, dipole, and neutral ones, respectively) in the plane of N ≡ N ± and G, for fixed g = −1 and q = 0.3.In region D above the short-dashed line in Fig. 8(c), unstable quadrupole solitons also exist, in addition to the stable dipoles.solitons exist too, but they are unstable against spontaneous transformation into quadrupole ones; see Fig. 9(b).

C. Collisions between traveling dipole solitons
The Galilean invariance of the underlying system of Eqs. ( 6)-( 8) makes it possible to set stable solitons in motion with velocity c, by applying a kick to them, ψ ± → ψ ± e icx .This way, collisions between solitons moving with velocities ±c may be simulated.We have considered two types of collisions of dipole solitons, with identical or opposite dipole moments, i.e., (DM, DM) or (DM, −DM).Obviously, the dipole-dipole interaction force is attractive in the former case, and repulsive in the latter one.Figure 10(a G = 0.9, L = 20, and N ± = 10.The collision is elastic, with the solitons readily passing through each other.On the other hand, Fig. 10(b) demonstrates that, for the same parameters, the solitons with opposite dipole moments undergo multiple collisions (ca.four in this picture), bouncing back and colliding again, until merger of the dipole-antidipole pair into a stable quadrupole soliton.It features a central ψ − component sandwiched between two ψ + side lobes; cf.Fig. 8(a).

IV. CONCLUSION
The objective of this work is to introduce 1D and 2D models of the "ultracold plasma," composed of two atomic species, which carry positive and negative charges (cations and anions), interacting via the electrostatic field governed by the Poisson equation.An estimate of parameters for a relatively dense degenerate plasma composed of light cations  6)-( 8) are G = 0.9, g = −1, and q = 0.3.The dipole-dipole collision is elastic, while the dipole-antidipole interaction leads to recurring collisions and, eventually, merger of the dipole-antidipole pair into a stable quadrupole soliton.and anions demonstrates that, in the "bare form," it would be too strongly coupled by the Coulomb interaction; however, reduction of the effective mass by means of a lattice potential may bring the system into the mean-field state described by the GP equations for wave functions of the ionic species, coupled to the Poisson equation for the electrostatic field.The GP equations include, along with the Coulomb terms, contact interactions, namely, the repulsion between the species and self-repulsion or attraction in each one.The contact interactions may be included directly, or induced by a buffer species of heavy neutral atoms (which may also be used for sympathetic cooling of the ions).The estimate demonstrates that, in the above-mentioned regime with the reduced effective mass, the contact interactions provide for relevant competition with the Coulomb forces and kinetic energy of the particles.For the system with full contact repulsion, the main objective is to predict stable spatially periodic density-wave patterns, which are akin to ionic crystals in solid-state physics.The transition from the uniform neutral state to the density wave is exactly identified by means of the analytical consideration of the MI (modulational instability) of the uniform state.The emerging 1D pattern is accurately predicted by means of the variational approximation, and found in the numerical form.It represents the system's GS above the onset of the MI.In the 2D setting, a stable density wave is found with the quasi-1D shape.
The 1D system with contact self-attraction in each ionic species gives rise to bright solitons.If the contact repulsion between the species is strong enough, neutral solitons split into dipole or quadrupole states, which represent the GS.The transition from the neutral solitons to dipole ones is accurately predicted by an analytical approximation.Different types of the solitons may coexist as stationary states, but only one of them is stable.Collisions between moving dipole solitons were simulated too.The result is that the dipole-dipole collision is elastic, while the dipole-antidipole pair features multiple collisions, eventually merging into a stable quadrupole soliton.
This work can be developed in other directions.In particular, as mentioned above, it may be relevant to add a wave function representing polar molecules built as cationanion bound states.The respective GP equations should then include "reactions" between colliding ions and molecules.A challenging extension is to address the 3D version of the present system, and to study the 2D case more systematically.In particular, it may be relevant to construct domain-wall patterns between two sets of stable stripes, such as those displayed in Fig. 5, but with different orientations [60,61].

FIG. 2 .
FIG. 2. (a) Stable spatially periodic numerical solution of Eqs.(15),(16), and (18) (the density wave) with components u + (x) and u − (x) (solid green and dashed blue lines, respectively) in the domain of size L = 10 with periodic boundary conditions.The chemical potential and norms of the numerical solution are μ ± = 94.85 and N ± = 325.(b) The fit of the numerical solution for u + (x) (the continuous line) to the analytical ansatz, u 0 + A cos(2π x/l + α) (the dashed line), with u 0 = 5.49 and A = 2.17 determined numerically as the spatial average value of u + (x), and the square root of the average of [u + (x) − u 0 ] 2 .The best-fit value of the phase shift is α = 0.68.The system's parameters are g = 1, G = 2, and q = 1.

FIG. 3 .
FIG. 3. (a)Modulation depth of the density wave, defined as per Eq.(29), vs the norm of each component, N ± = N, for the same system's parameter as in Fig.2, i.e., g = 1, G = 2, q = 1, and L = 10.The dashed line shows A(N ) as predicted by analytical approximation(35) with l = L/12.(b) The chain of rhombuses and dashed line depict, severally, the numerically found chemical potential of the same family of the density waves, μ(N ), and its analytically predicted counterparts, given by Eq.(37).In addition, the chain of crosses shows μ as predicted, in a simple form, by Eqs.(36) and (32),(35).
FIG. 5. Stable 2D striped pattern with norms N + = N − = 650, produced by the imaginary-time solution of Eqs.(6)-(8) in a square of size 5 × 5 with periodic boundary conditions; see the definition of black stripes given by Eq.(43).The parameters are G = 2 and g = q = 1.

FIG. 8 .
FIG. 8. (a) Stationary profiles of |ψ + | and |ψ − | − (solid green and dashed blue lines, respectively) for a stable quadruple soliton with norms N ± = 5, obtained by means of the imaginary-time method as a solution of Eqs.(6)-(8) in the domain of size L = 20, with Neumann's boundary conditions.The parameters are g = −1, G = 0.9, and q = 0.3.The quadrupole moment of the soliton [see Eq. (60)] is QM = −13.54.(b) The energy difference, E − E 0 , between quadrupole and dipole solitons (solid and dashed lines, respectively) and the neutral one(49), for fixed norms, N ± = 5.25, vs G, keeping g = −1 and q = 0.3 fixed too.The GS (ground state) is identified as one with the lowest energy.(c) The chart which identifies GS as different soliton species (Q, D, and N: quadrupole, dipole, and neutral ones, respectively) in the plane of N ≡ N ± and G, for fixed g = −1 and q = 0.3.In region D above the short-dashed line in Fig.8(c), unstable quadrupole solitons also exist, in addition to the stable dipoles.
FIG. 9. (a) Spontaneous transformation of an unstable quadrupole soliton with norm N ± = 5.25 into a stable dipole one, at g = −1, q = 0.3, and G = 0.65.(b) The transformation of an unstable dipole into a stable quadrupole at the same values of q, g, and N ± , but with larger G = 0.75.In terms of the chart in (c), parameters corresponding to cases (a) and (b) belong to areas D (above the internal short-dashed boundary) and Q (slightly above its bottom boundary), respectively.The evolution is displayed by means of profiles of |ψ ± (x, t )|.

FIG. 10 .
FIG. 10.Panels (a) and (b) display, severally, typical examples of dipole-dipole and dipole-antidipole soliton-soliton collisions, with velocities c = ±1.In both cases, norms of the solitons are N ± = 10 and the size of the domain is L = 20.The parameters in Eqs.(6)-(8) are G = 0.9, g = −1, and q = 0.3.The dipole-dipole collision is elastic, while the dipole-antidipole interaction leads to recurring collisions and, eventually, merger of the dipole-antidipole pair into a stable quadrupole soliton.