{'fraction_non_alphanumeric': 0.04539815460959913, 'fraction_numerical': 0.015933938125145382, 'mean_word_length': 4.35054967848994, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Effects of Long-Range Nonlinear Interactions in Double-Well Potentials

We consider the interplay of linear double-well-potential (DWP) structures and nonlinear longrange interactions of different types, motivated by applications to nonlinear optics and matter waves. We find that, while the basic spontaneous-symmetry-breaking (SSB) bifurcation structure in the DWP persists in the presence of the long-range interactions, the critical points at which the SSB emerges are sensitive to the range of the nonlocal interaction. We quantify the dynamics by developing a few-mode approximation corresponding to the DWP structure, and analyze the resulting system of ordinary differential equations and its bifurcations in detail. We compare results of this analysis with those produced by the full partial differential equation, finding good agreement between the two approaches. Effects of the competition between the local self-attraction and nonlocal repulsion on the SSB are studied too. A far more complex bifurcation structure involving the possibility for not only supercritical but also subcritical bifurcations and even bifurcation loops is identified in that case.

by Refs. [28][29][30], as well as the interactions with Gaussian and exponential kernels, motivated by nonlinear-optical models [34]. In particular, we consider the effect of the range of the interaction, with the objective to consider a transition from the contact interactions to progressively longer-range ones. We conclude that the phenomenology of the short-range interactions persists, i.e., the earlier discovered SSB phenomena [12,14,15,35] still arise in the present context. However, the critical point of the SSB transitions features a definite, monotonically increasing, dependence on the interaction range, which is considered in a systematic way. More elaborate scenarios can be detected in the case where in addition to the long-range interactions, there is a competing short-range component. In such a case, we identify not only the earlier symmetry-breaking bifurcations but also reverse, "symmetry-restoring" bifurcations, as well as the potential for symmetry-breaking to arise (for the same parameters) both from the symmetric and from the antisymmetric solution branch. Both of these are phenomena that, to the best of our knowledge, have not been reported previously. Interestingly, the only example where subcritical bifurcations have been previously discussed in the double well setting for GP equations is the very recent one of extremely (and hence somewhat unphysically) high nonlinearity exponents in [36].
The presentation is structured as follows. In section II, we present the model and the quasi-analytical two-mode approximation, which clearly reveals the system's bifurcation properties. In section III, we corroborate these analytical predictions by full numerical results for different types of the long-range kernel. In section IV, we consider a more general case, in which the long-range nonlinear repulsion competes with the local attraction. There, we illustrate how the competition strongly affects the character of the SSB in the DWP setting. Finally, in section V, summarize our findings and present our conclusions.

II. THE LONG-RANGE MODEL AND THE ANALYTICAL APPROACH
In the quasi-1D setting, the normalized mean-field wave function ψ(x, t) obeys the scaled GP equation, where µ is the chemical potential, and is the usual single-particle energy operator, which includes the confining DWP withΩ the normalized harmonic-trap's strength;Ω ≪ 1 in a quasi-1d situation in BECs. The nonlinear term with coefficient s accounts for the long-range interatomic interactions, s = ±1 corresponding to the repulsion and attraction, respectively. Note that the contact interaction is not taken into regard in Eq. (1), as we aim to focus on the effect produced by the long-range nonlinearity (in the gas of 52 Cr atoms, the contact interaction may be readily suppressed by means of the Feshbach resonance [21]). In this work, we consider mainly symmetric spatial kernels in Eq. (1), that are positive definite, with the following three natural forms chosen for detailed analysis: the Gaussian, the exponential, and the cut-off (CO) (alias generalized Lorentzian) kernel, The width of the kernels, σ, determines the degree of the nonlocality. All three kinds of the kernels go over into the δ-function as σ approaches zero, in which case Eq. (1) turns into the usual local NLS/GP equation. All the kernels are normalized as the δ-function, i.e., +∞ −∞ K(x)dx = 1, and the norm of the wave function will be used in the usual form, In what follows below, we adopt typical physically relevant values of the scaled parameters, namely,Ω = 0.1, V 0 = 1 and W = 0.5, in which case the two lowest eigenvalues of linear operator L with potential (3) are numerically found to be ω 0 = 0.1328 and ω 1 = 0.1557.

A. The two-mode approximation
The spectrum of the underlying linear Schrödinger equation (s = 0) consists of the ground state, with wave function u 0 (x), and excited states, u l (x) (l ≥ 1). In the weakly nonlinear regime, the Galerkin-type two-mode approximation is employed to decompose the wave function ψ(x, t) over the minimum basis constituted by the ground state u 0 and the first excited state u 1 , associated to the eigenvalues ω 0 and ω 1 , respectively. For this purpose, it is more convenient to use a transformed orthonormal basis composed by wave functions centered at the left and right wells, viz., Without the loss of generality, ϕ L,R are both chosen to be positive definite. Thus, the two-mode approximation for the wave function is defined as where c L and c R are complex time-dependent amplitudes. Substituting this into Eq. (1), we obtain where the asterisk and overdot stand for the complex conjugate and time derivative, while Ω ≡ (ω 0 + ω 1 )/2 and ω ≡ (ω 1 − ω 0 )/2 are linear combinations of the two lowest eigenvalues. Next, we project Eq. (8) onto the single-well states ϕ L,R , which involves the following overlap integrals: The equations hold if subscripts L and R are swapped, or the variable x and x ′ are interchanged, due to the symmetry of kernel K. In Fig. 1 we show the values of the four integrals, η 0,1,2,3 , as functions of parameter σ, for the three types of kernels defined in Eqs. (4)- (6). We note that η 2 and η 3 remain negligible for any value of σ, and η 0,1 are both positive when the kernel K is positive definite. Naturally, when σ is small, which corresponds to the limit of a nearly local nonlinearity, η 0 is much larger that the other three overlap integrals, due to the weak overlapping of the single-well states ϕ L,R . As σ increases, η 0 decreases while η 1 increases and, finally, they tend to become equal.
Regarding the values of these integrals, the upcoming analysis is conducted in two situations: η 0 much larger than all others, and η 0 being on the same order of magnitude as η 1 . Some value of the σ b is to be fixed as a boundary between the two cases. Since η 0 is always large, the key point is to set up a rule for comparing η 1 with η 2,3 . To this end, we define η rel = η 1 − max(|η 2 |, |η 3 |) and the criterion is stated as follows: if η rel ≥ 0.01, η 1 is taken into regard in the analysis of the two-mode approximation; otherwise, η 1 is insignificant, and only η 0 is kept. For the kernels considered in this context, the criterion yields σ b = 2.96, 1.56 and 1.91 for the Gaussian, exponential and Lorentzian kernels, respectively. Thus, with σ < σ b , all integrals η 1,2,3 are omitted, and only η 0 is retained. Then, the projection of Eq. (8) onto the two-mode set, ϕ L,R , leads to the following ODE system: We then introduce the Madelung representation, c L,R = ρ L,R e iθL,R with real time-dependent ρ L,R and θ L,R , and derive from Eq. (10) a set of equations for ρ L,R and θ L,R : where θ ≡ θ L − θ R is the relative phase between the two modes, and the equations for ρ R and θ R are obtained by interchanging subscripts L and R, and θ with −θ, in Eq. (11) directly. We focus on steady solutions to this system, i.e.ρ L,R =θ L,R = 0. Then, for solutions with nonzero amplitudes, θ may only take values 0 or π, which correspond, respectively, to equal or opposite signs of real stationary solutions for c L and c R . Through a straightforward algebra, three stationary solutions are thus found: the symmetric solution, with θ = 0 and ρ 2 L,R = (µ − ω 0 )/sη 0 , existing when µ > ω 0 (µ < ω 0 ) for s = 1 (s = −1), i.e., for the repulsive (attractive) long-range interactions; the antisymmetric solution, with θ = π and ρ 2 L,R = (µ − ω 1 )/sη 0 , existing when µ > ω 1 (µ < ω 1 ) for s = 1 (s = −1); and an asymmetric solution, with ρ 2 As we assume η 0 > 0, for the repulsive case (s = 1) the asymmetric state exists only if θ = π, i.e., it bifurcates from the antisymmetric solution when µ > Ω + 2ω. On the contrary, in the attractive case, it emerges from the symmetric state (θ = 0) when µ < Ω − 2ω. These conclusions agree with the general principles of the SSB theory, according to which the attractive/repulsive nonlinearity breaks the symmetry/anti-symmetry of solutions with equal numbers of particles in the two wells [12,14,15,35].
The two-mode approximation is a powerful means for identifying different steady states of the underlying problem. The solutions found above are also used as initial conditions in solving the full NLS equation (1), as reported in section III.

B. The bifurcation analysis
The two-mode approximation strongly facilitates the qualitative analysis of the SSB bifurcation, as well as exploring the system's dynamics. To proceed, we define the population imbalance between the two wells, where N L,R = |c L,R | 2 ≡ ρ 2 L,R , hence the total norm is N = N L + N R . Together with the relative phase, θ = θ L − θ R , we eventually derive the following dynamical equations: This form of the equations is relevant for both cases, when only η 0 or both η 0 and η 1 dominate, as discussed before. Note that η stands for η 0 in the former case, and for ∆η in the latter one. Equations (15) take the Hamiltonian form, with Hamiltonian Equations (15) possess the stationary solutions, (z 1 , θ 1 ) and (z 2 , θ 2 ), with z 1 = z 2 = 0, θ 1 = 0, θ 2 = π, which represent the symmetric and the antisymmetric solutions, respectively. Besides that, the asymmetric stationary solutions may exist when N ≥ N cr , with N cr = |2ω/η|, taking the form of for s = 1 (s = −1). The asymmetric solution emerges from the antisymmetric (symmetric) one through a pitchfork bifurcation, in the case of the defocusing (focusing) nonlinearity. Since η = η 0 for σ < σ b and η = η 0 − η 1 for σ > σ b , Fig. 1 suggests that η is, generally, a decreasing function of σ, and consequently N cr is increasing with respect to σ. This way, the bifurcation takes place at larger value of N when the the nonlocality range is wider, which is consistent with the conclusion concerning µ cr . Next, from Eq. (15) we derive the equation of motion, which leads to the system ż = p, When N < N cr , there is a unique stationary solution, p = z = 0, which is a fixed point of the center type. When N > N cr , the origin, (0, 0), becomes a saddle, with another pair of fixed points (centers), p = 0, z = ± 1 − 4ω 2 |η| 2 N 2 , representing the asymmetric solutions. Figure 2 shows the phase space of system (20), along with the linearization near the fixed points for the example of the Gaussian kernel with σ = 5.0 and N = 0.5, in which case N cr = 0.33.

III. NUMERICAL APPROACH
We first consider the repulsive interaction case (s = 1). Branches of stationary solutions to the full partial differential equation (1) are explored for all three kernels defined in Eqs.  where operators L 1 , L 2 are defined as for any function φ. The stationary state is called unstable if there exist any eigenvalues λ with ℜ(λ) = 0, otherwise it is stable (i.e. all corresponding eigenvalues are purely imaginary). We present the results of Gaussian kernel in detail as an example. In each panel (Fig. 3), the solid blue line with highest value of N for any µ among the three branches is the symmetric stationary solution; the continuation of it to the linear limit (N → 0) shows that it starts from µ = ω 0 (the eigenvalue associated to the ground mode of the underlying linear system), and it is stable for any µ. The branch with a part of it denoted by dashed red line is the antisymmetric solution, arising from the linear limit at µ = ω 1 , corresponding to the first excited state. It starts as a stable state from the linear mode (denoted by solid blue line) and, as µ increases, after some critical point µ cr , it is destabilized due to the emergence of the asymmetric branch through a supercritical pitchfork; notice that the asymmetric branch remains stable after it occurs. This bifurcation was predicted by the two-mode analysis in the previous section.
The four panels of Fig. 3 are the bifurcation diagrams of the model with the Gaussian kernel, for σ = 0.1, 1, 5 and 10. The main parts of the panels display the numerically found solutions of the full system, while the small plot in each panel shows the corresponding analytical solutions predicted by the two-mode approximation, as obtained in Section II. The numerical and analytical branches of the solutions demonstrate good agreement in all the four cases. Recall that we have obtained the critical values at which the bifurcation takes place, µ cr (or equivalently N cr ), in an analytical form. As for the Gaussian kernel (with the chosen border at σ b = 2.96), σ = 0.1 and σ = 1 are categorized as belonging to the first case, with solely η 0 taken into account, among all the overlap integrals. The analytically predicted value is µ cr = 0.1671 (since µ cr = Ω + 2ω), while its numerically found counterpart is µ cr = 0.1684 for both σ = 0.1 and 1. On the other hand, σ = 5 and σ = 10 pertain to the the second case, in which both η 0 and η 1 are kept. In this case, the analytical prediction is µ cr = Ω + 2ωη 0 /∆η, which depends on η 0,1 and thus varies for different values of σ. The predicted values of µ cr is 0.1757 and 0.2120, for σ = 5 and 10, respectively, while the respective numerical values are µ cr = 0.1745 and 0.2075. Despite a small discrepancy between the two groups of the values (numerical versus analytical), the analytical results still succeed in predicting the trend of the behavior of µ cr , viz., µ cr increases with the growth of σ, or, in other words, the SSB bifurcation takes place at higher values of N cr , as shown in Fig. 4.
Next we consider cases with the long-range interaction based on the other two kernels, viz., the exponential and Lorentzian ones, given by Eq. (5) and (6), respectively. The corresponding bifurcation diagrams, similar to the case of the Gaussian kernel, are presented in Fig. 5 and Fig. 6, with the same notations as adopted above.
Lastly, we also briefly discuss the self-attractive case, with s = −1 in Eq. (1). Two examples are selected, using the Gaussian kernel with σ = 1 and 5, to help illustrating the similarities and differences with the self-repulsive case. The complete bifurcation diagrams are plotted, for this case, in Fig. 7, with the notation similar to that introduced above.
In this case, the two branches arise from the linear modes at µ = ω 0 and ω 1 , pertaining to the symmetric and antisymmetric stationary solutions, which exist for µ < ω 0 or µ < ω 1 , respectively. This time, the supercritical SSB pitchfork bifurcation occurs on the symmetric branch, leading to the emergence of the asymmetric state. The pitchfork is found at µ cr = 0.1213 and 0.1128 for σ = 1 and 5 respectively, in good agreement with the critical values predicted by the two-mode approximation, which are 0.1212 and 0.1104. In each panel of Fig. 7, the diagram produced by the two-mode approximation, corresponding to its numerical counterpart, is displayed in the top right corner, exhibiting a very good agreement between the two. Here and in the next figure, the value of σ is chosen as σ = 1 (left), σ = 5 (right), and the notation is the same as in Fig. 3.

IV. COMPETING SHORT-AND LONG-RANGE INTERACTIONS
In this section, we consider a model illustrating effects of the competition of the long-range interactions with local ones, based on the following GP equation, where coefficient g accounts for the contact nonlinearity. The interactions are fully attractive or repulsive if s and g are both negative or positive, respectively. However, in such cases the results turn out to be similar to those reported above for the nonlocal nonlinearity. A more interesting case arises for sg < 0, i.e., for the competing interactions. Effects of the competition on 1D solitons were recently studied in detail in Ref. [29]. We use the analysis presented in the previous section as a starting point for our considerations here. Thus, we fix the long-range interaction coefficient s = ±1, and then vary the local interaction coefficient, g from 0 to −1 (g < 0 implies the local attraction), using Gaussian kernel (4) with σ = 5 and s = 1 as a case example. As the attractive local interactions grow stronger, interesting phenomena emerge in the development of the corresponding bifurcation diagram, as is shown in Fig. 8.
When g is close to zero, the three branches of stationary states are similar to those in the case of pure long-range interactions, as shown in the bottom left panel of Fig. 3. An example for g = −0.2 is presented in the top left panel of Fig. 8, where the norm N of each branch is larger at any value of µ, in comparison to the case of g = 0. Also, the bifurcation takes place at a larger critical value, µ cr = 0.1805.
As g grows more negative, we notice that the symmetric and the antisymmetric branches, while monotonically increasing with µ at small values of N , switch to become monotonically decreasing functions of µ at a large value of N . Simultaneously, we observe the bifurcation of two asymmetric states, each from a different branch. While it  −1), as a function of µ. The Gaussian kernel is taken, with σ = 1 (left) and σ = 5 (right). The notation is the same as in Fig. 3. is commonly known that attractive and repulsive interactions favor bifurcations from symmetric and anti-symmetric branches, respectively, here we encounter the first (to our knowledge) example that features bifurcations from both branches. This unusual situation extends up to g = −0.32. The complete bifurcation diagram is displayed in the top right panel of Fig. 8, for the case of g = −0.3. In this case, the symmetric and antisymmetric branches arise, as usual, from their linear limits at µ = ω 0 and µ = ω 1 , respectively. A supercritical pitchfork bifurcation occurring (as before) on the antisymmetric branch gives rise to an asymmetric state at µ = −0.195, destabilizing the antisymmetric one. What is, however, different here is that this asymmetric state "survives" only within a narrow parametric interval before it merges back into the antisymmetric branch through another subcritical ("symmetry-restoring") pitchfork at µ = 0.3263. This bifurcation loop is reminiscent of that reported previously for two-component solitons in local 1D and 2D models with competing self-focusing cubic and self-defocusing quintic nonlinearities [37]. For higher µ, another asymmetric branch arises from the symmetric one through a subcritical pitchfork at µ = 0.3465; it is especially important to highlight that this pitchfork, which destabilizes the symmetric branch, is a subcritical one, hence the asymmetric state is unstable too. It is observed that, after its emergence, the norm of the asymmetric branch decreases as µ decreases (dN/dµ > 0), before starting to rise at the fold point at µ = 0.2768, featuring dN/dµ < 0 after that. The stability of the state changes at the turning point, which is in agreement with the well-known Vakhitov-Kolokolov (VK) criterion [38], according to which the slope of the branch determines its stability. In Fig. 9, a blowup provides a clearer view of these bifurcations. The symmetric and the antisymmetric branches continue to increase their norms with the growth of µ, each of them featuring its own turning point, at µ = 0.391 and µ = 0.3853 respectively, and change their monotonicity thereafter. However, this change on the antisymmetric branch is not accompanied by a stability change. To demonstrate that this is not an exception to the VK criterion, we define the linearization operators, using notation L 1,2 defined in Eq. (22). It is known from Ref. [39] that the instability arises, when the slope condition of the VK criterion is violated, if |n(L + ) − n(L − )| = 1, n standing for the number of negative eigenvalues of each operator. The state under consideration is generically unstable if |n(L + ) − n(L − )| > 1. In our case, the antisymmetric branch has n(L − ) ≡ 1 (due to its single-zero-crossing configuration profile, which is a zero mode of L − ). Meanwhile, n(L + ) = 1 is true before the branch turns to the left, and consequently |n(L + ) − n(L − )| = 0. Hence, even though the slope condition is violated (dN/dµ > 0), the relevant theorem suggests that the VK criterion does not actually apply here (and no stability change occurs). After the turning point, the stability is expected since n(L + ) = 2, and the slope condition is satisfied, which is in accordance with our observations. Next, as the local interactions continue to grow stronger, another regime is explored, as shown in the bottom left panel of Fig. 8, with g = −0.35 as an example. In this case, the phenomenology of the symmetric branch is similar to that in the previous situation; however, the two sequential bifurcations taking place on the antisymmetric branch no longer arise. Hence, the entire antisymmetric branch remains dynamically stable. The asymmetric branch still emerges from the symmetric one through a subcritical pitchfork, and later becomes stable past the fold point due to the change of slope, in accordance with the VK criterion.
The above regime persists until g = −0.38. After that, the pitchfork bifurcation occurring on the symmetric branch switches from subcritical to supercritical, leading to the emergence of a stable asymmetric branch, as shown in the last panel of Fig. 8. This shape of the bifurcation diagram, consisting of the three branches, persists as far as g decreases to −1. Note that during this process, the turning points of the symmetric and antisymmetric branches keep getting lower (with respect to N ), finally leading to a situation where the turning points do not exist, and the two branches immediately go left starting from their linear limits, i.e., the symmetric and the antisymmetric solutions only exist for µ < ω 0 , ω 1 when g < −0.69, −0.6, respectively.

V. CONCLUSIONS AND FUTURE CHALLENGES
In this work, we have presented a systematic study of the interplay between the long-range nonlinear interactions (of either sign, repulsive or attractive) and linear double-well potential (DWP) in the 1D setting. The two-mode approximation has been developed, that accounts for the nonlocality but retains the structure similar to that established before in the case of the contact interaction. This conclusion demonstrates that the fundamental phenomenology of the spontaneous symmetry breaking (SSB) in the DWP persists in the presence of the longer-range interactions, although the critical properties themselves (e.g., the critical values of the chemical potential and norm), at which the SSB bifurcation occurs, giving rise to the asymmetric steady states, are sensitive to the precise interaction range. In particular, our analysis has revealed a monotonous increase of the critical values as a function of the interaction range, in the case of the self-repulsion. A considerably more elaborate phenomenology was revealed by the analysis in the context of competing long-and short-range interactions. In that case, a delicate interplay between the strengths of the nonlocal repulsion and local attraction gives rise, in addition to predominantly attractive and predominantly repulsive regimes, to mixed ones with complex bifurcation phenomena. On the one hand, symmetry-breaking effects were shown to arise from each of the relevant solution branches; in some cases, they are accompanied by reverse bifurcations, to form closed loops. On the other hand, some of the previously studied supercritical bifurcations, such as the one occurring on the symmetric branch, could become subcritical, being subsequently coupled to additional fold bifurcations. All of these effects are absent, to our knowledge, in models with a single cubic nonlinear term, being consequences of the competition.
These results may be a motivation for studies in a number of future directions. In particular, it would be relevant to consider the generalizations of the DWP setting to higher dimensions, such as, the four-well configuration, which was recently demonstrated to yield a much richer phenomenology in the case of the contact interactions [40]. In the same connection, it is relevant to note that the few-mode approach, used in this work for the consideration of the existence and stability of the symmetric, antisymmetric, and asymmetric states, as functions of the interaction strength, may be applied to other structures, such as multi-dimensional bright [25], dark [41] and vortex [31] solitons, especially in cases where such structures are stabilized by nonlocal nonlinearities. A natural objective of such an analysis may be to unravel effects of the interaction range on structural stability of nonlinear waveforms. Such studies are presently in progress and will be reported elsewhere.