{'fraction_non_alphanumeric': 0.03729815019975157, 'fraction_numerical': 0.01356296370900057, 'mean_word_length': 4.5162962962962965, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{\phi}^4 Solitary Waves in a Parabolic Potential: Existence, Stability, and Collisional Dynamics

We explore a {\phi}^4 model with an added external parabolic potential term. This term dramatically alters the spectral properties of the system. We identify single and multiple kink solutions and examine their stability features; importantly, all of the stationary structures turn out to be unstable. We complement these with a dynamical study of the evolution of a single kink in the trap, as well as of the scattering of kink and anti-kink solutions of the model. We see that some of the key characteristics of kink-antikink collisions, such as the critical velocity and the multi-bounce windows, are sensitively dependent on the trap strength parameter, as well as the initial displacement of the kink and antikink.


Introduction
Models of the sine-Gordon (sG) [4] and more generally of the nonlinear Klein-Gordon type, such as φ 4 [5] have been studied intensely via a combination of analytical and numerical techniques for well over three decades now. Part of the reason for this widespread appeal concerns the diverse set of physical applications for which such models have been argued to be of relevance. These start from the simplest (coupled torsion pendula, motion of dislocations etc. in the case of the discrete sG [6]) and extend far and wide. Case examples include, but are not limited to structural phase transitions [7,8], domain walls arising in cosmological models [9,10], simple polymeric models, as well as models of uniaxial ferroelectrics; see, e.g., Ch. 9 in [11].
However, there is another source of the fascination, especially from an applied mathematics perspective, in non-integrable models such as φ 4 . This has to do with their remarkable properties related to the collisional dynamics of kinks and antikinks. There, resonant phenomena between translational, internal and extended modes have been identified and and a phenomenological "resonance energy exchange mechanism" [2] has been used to characterize features such as the "multi-bounce windows", the fractal emergence of such windows, and the universal nature of these phenomena not only in φ 4 , but also in other models such as the double, and the parametric sG equations [2,12,10]. While recent years have seen both a more rigorous mathematical analysis [13] and an experimental realization of a system reminiscent of the effective kink-antikink interactions [14], there have also been further studies causing a renewed interest in the phenomenological approach and towards posing new questions. For instance, particular variations of the φ 6 model have been shown numerically to possess multi-bounce windows without possessing internal modes [15]. It has also been found that by parametrically deforming the φ 4 model, one can introduce more internal modes and suppress the two-bounce windows [16]. Lastly, the recent work of [17,18] has performed a careful bookkeeping of the collective coordinate approaches intended to provide an effective description of the collision events. In so doing, a number of nontrivial problems regarding both computing with such reduced models and the conclusions drawn from them have arisen. Thus, despite its time-honored tradition, the study of kink-antikink collisions in non-integrable models remains a surprisingly active field of investigation to this day [19,20,21,22,23].
In the present work, we examine a different variant of the φ 4 problem, considering the scenario where an external potential has been imposed on the field theory. In particular, motivated by the extensive work on atomic Bose-Einstein condensates (BEC), we consider a parabolic trap, which in the latter context emulates the role of a generic magnetic confinement [24,25,27]. This possibility turns out to alter dramatically the properties of the kink waveforms. In particular, even the single kink turns out to be dynamically unstable in the emerging energy landscape. We examine the case of the single kink and its existence, stability and dynamic properties, as well as that of the kink-antikink state. In the latter, we examine the kink-antikink collisions and attempt to capture them by means of low dimensional collective coordinate models. We highlight both the favorable, as well as the unfavorable traits of these models with respect to their ability to capture the full PDE dynamics and raise some relevant possibilities towards future studies.

The Model, Calculational Approach and Ground State
Our starting point for the considerations that follow is the classical φ 4 equation with an added x-dependent trapping term: where the potential function is given by Ω is the trap strength constant and typically (in settings related to our motivating example of atomic BECs [24,25,27]) assumes values in (0, 1). However, we mostly restrict our attention to the dynamics for Ω ≪ 1, using Ω = 0.15 for most of the data runs 4 . The case Ω = 0 has been well-studied in the past and in its analysis, the perturbation theory and collective coordinate methods have been used [10,5,2,3,17]. We extend these methods to study our φ 4 variant bearing a trapping in different parameter regimes.
To explore the properties of Eq. (1) we employ the following numerical scheme. First, we use a fourth-order central difference in space (see e.g. [10,28] for this and other similar difference schemes). The resulting system of ordinary differential equations is then evolved in time using fourth-order Runge-Kutta methods, ensuring the accuracy of the conservation of the energy. The average value of the energy is of O(1), while the deviations from the mean are no more than O(10 −4 ). We apply free boundary conditions, although as we shall see below this is generally inconsequential because the field effectively vanishes well before the edge of our computational domain. For all our data runs, a spatial step size of ∆x = 0.02 was used. Our computations typically use Ω = 0.15 and a spatial grid running in the interval [−30, 30]. For smaller values of Ω, we use a larger interval as necessary to ensure the decay of the wavefunction within the computational domain.
The stationary ground state of the system uΩ(x) solves the ordinary differential uΩ is not known explicitly but in the case Ω ≪ 1, it is well-estimated by the Thomas-Fermi (TF) approximation which neglects the spatial variation of u (again motivated by similar considerations in atomic condensates [24,25]) defined by Eq. (4) solves Eq. (3) when the uxx term is neglected. We use this approximation as our initial guess for the Newton-Raphson iterations converging to the true numerical ground state of the system. The comparison of the TF approximation with the true numerical result is given in Fig. 1. At the level of the TF approximation, we note that the ground state is zero whenever |x| > xs, where uTF(xs) = 0, which gives As we can observe in the figure, the TF approximation is very good except for the vicinity of xs, where its disagreement with the true ground state has been quantified in a series of rigorous studies [29,30].
To study the stability of the ground state, we consider the linearization around the solution using the ansatz: Here uΩ is the numerically obtained ground state, and (λ, χ) represent the corresponding eigenvalues and eigenvectors of the linearization. The equation that these satisfy at O(ǫ) is of the form: achieve effective one-dimensionality of the system the condition Ω ≪ 1 needs to be enforced. This is what motivates our selection here too.
For the ground state we find, as is shown in the right panel of Fig. (1) that all the eigenvalues lie on the imaginary axis, hence the relevant state is spectrally stable; for instability to arise, at least one of the relevant eigenstates should correspond to growth, i.e., to Re(λ) > 0.

Single Kink Solutions
We now turn to excited states in the system in the form of kinks and anti-kinks, as well as bound states thereof. We start with the single static kink whose solution is not known explicitly but can be numerically approximated by once again in the TF approximation. For various values of Ω, the numerically calculated solutions are shown in Figure 2.
In the figure, also the spectrum of linearization around the single kink is shown (right panel). We see that in all the cases of nonzero Ω, there exists a pair of nonzero real eigenvalues. More concretely, the eigenvalue pair that used to be at λ = 0 due to translational invariance now exists as real and leads to the instability of the kink, which is progressively stronger, the larger the value of Ω.

Numerical Results (PDE)
We now examine the dynamical evolution of such unstable kinks. As a means of obtaining moving kink solutions uK (or antikink solutions uK ) of (1), we can employ the Lorentz transformation, where x0 is the initial position and v is the velocity of the moving kink. The initial conditions of the PDE are Of course the model in the presence of the potential is no longer Lorentz invariant (as it is in the case of Ω = 0). Nevertheless, our results have shown that from a numerical perspective producing a Lorentz-boosted kink and multiplying it by the (stationary) TF background is an efficient method to produce "moving" waveforms starting from standing ones.
The presence of the the trapping term causes the velocity of the kink not to remain constant (as it would under stationary initialization in the untrapped, translationally invariant problem). When the initial velocity (vin) is zero, the numerical results show that the motion of the kink solution depends sensitively on its starting position. When the kink starts to the left of the origin (x0 > 0), the kink slides to the left, accelerating as it does so. If the kink starts to the right of the origin (x0 < 0), then it slides to the right. In both cases, the kink eventually ends up being expelled from the system. Fig. 3 shows examples of initializing the kink both to the left and to the right of the fixed point at x = 0 with no speed and illustrates how the kink slides along the side where it starts. This corroborates in a definitive way the results of the stability analysis in that the kink now (in the presence of the parabolic trap) encounters a saddle point at x = 0 (the center of the trap). Hence, on each side of this unstable equilibrium, it will slide along the unstable manifold. In fact, it is interesting to note that in panel (c) of the figure the kink is initialized at x = 0 (up to numerical roundoff error). We can see that in that case, it stays at x 0 = 1, the kink slides away from zero to the right, because x 0 is positive and gradually accelerates. (c) With x 0 = 0, the kink preserves its steady form for a while but eventually as its instability manifests, it slides (in this particular case) to the right.
the origin for a long time, until eventually the projection of the roundoff error in the unstable eigendirection grows to O(1) and kicks the kink away from the unstable equilibrium (in the case of this particular example its center moves to x > 0, but it can just as well move to the x < 0, depending on the perturbation).
With nonzero initial velocity, the behavior of the kink depends on the magnitude of the velocity. For a kink starting away from the origin, the saddle point at the x0 = 0 represents a finite (potential) energy barrier which the kink must overcome to break through to the other side. Given sufficient kinetic energy, a kink moving towards the origin will overcome this potential energy barrier. On the other hand, if the velocity (and hence the kinetic energy of the kink is "subcritical", i.e., below the value needed to overcome the potential barrier) then the kink will bounce off the barrier and return to be expelled on the side that it started. Fig. 4   demonstrate these two scenarios in a convincing way numerically. In particular, between vin = 0.577 and vin = 0.578 of the left and right panel, respectively, there is clearly a critical velocity as the former results in a bounce-back of the kink from x = 0, while the latter leads to its transmission through the barrier and on to the other side. For any vin < 0.577, the kink stops at x = x1, called a turning point and then moves in the opposite direction. The Fig. 4 (c) shows the relation between the effective potential approximated by the associated kinetic energy (per unit mass) 1 2 v 2 in and the turning point x1. Of course, in a way consistent with the above picture, if the kink is initially moving away from the origin, then it will continue to move away and gradually accelerate over time.

Collective Coordinate Approach (ODE)
In this section we formulate a "collective coordinate" (CC) approach for the dynamics for the single kink. Our aim is to reduce the full PDE with infinitely many degrees of freedom to an ODE model with a single degree of freedom. This is partly driven by previous CC studies of the motion of a kink in the homogeneous case (Ω = 0), including the work of Sugiyama [3] (see also [10]). It is also importantly motivated by the nature of the motion in the vicinity of the saddle point discussed above, and by the desire to obtain a simplified description of this motion.
To begin, we note that the full PDE conserves the energy (Hamiltonian) where the kinetic and potential energies of the kink, respectively, are The corresponding Lagrangian is Assuming that we operate in the TF limit of Ω ≪ 1, we seek to characterize the kink motion by utilizing the ansatz where X(t) is the time-dependent displacement of the kink from the origin and uΩ is the ground state of the system in the presence of the trap. Then, the Lagrangian L becomes as follows: Since we have an extremely accurate numerical solution for uΩ, we compute the coefficients a0 and a1 and solve the Euler-Lagrange equation numerically. Fig. 5 shows the plots of a0 and a1 as a function of X for various values of Ω. By applying the Euler-Lagrange equation, we obtain the dynamical evolution: We solve these equations numerically by using the initial conditions X(0) = x0 and X ′ (0) = vin where x0 is the initial position of the kink and vin is the initial velocity of the kink. We use MATLAB's built-in fourth-order Runge-Kutta variable-step size solver ode45 with built-in error control. Note that here, overdots denote time derivatives.

Connection to Numerical Results
We now compare our full numerical results (PDE) with the findings of the CC ODE Method. The numerical CC computations show that when the initial velocity is zero (vin = 0), with nonzero displacement (x0 < 0 or x0 > 0), the ODE results agree with PDE results quite well up to the time/location when/where the kink is expelled from the system; from there on the kink cannot be clearly discerned anyway, and hence it is not particularly meaningful to seek to track it in the PDE or to match it to the ODE results. Fig. 6(a) and (c) show that when x0 = −3 and x0 = 3, with zero initial velocity, both results mostly agree until t ≈ 11 when the kink is expelled from the TF background.
In the case of zero displacement (x0 = 0) and zero velocity (vin = 0) of Fig. 6(b), both the ODE and the PDE show the kink as residing on the saddle point for a long time interval. In the ODE case, the absence of additional degrees of freedom and the conservation of the energy does not allow the kink to depart from the unstable equilibrium over the time horizons considered. In the infinite dimensional system, the instability is manifested at t ≈ 250. When the initial velocity is nonzero, we find similar results. The ODE results agree with the PDE results until the ODE kink is expelled from the system. Fig. 7 (a) exhibits the ODE and PDE agreement until t ≈ 15 when x0 = 0 and vin = 0.2. Fig. 7 (b) exhibits this type of agreement until t ≈ 23 when x0 = −5 and vin = 0.5. In summary, we can infer that while the kink remains within the region where the TF background is accurate, the CC method appears to yield an excellent qualitative and a very good quantitative characterization of the resulting motion.

Kink-Antikink Solutions
Armed with our understanding of the single kink energy landscape, we now turn to case of a kink-antikink pair. Such a state can in fact be stationary in our model, although it is never genuinely stationary in the case of Ω = 0. The reason is fairly similar to the corresponding existence of a stationary state in a nonlinear Schrödinger (NLS) model with a trap analyzed in detail, e.g., in [31]. However, there is a crucial difference. In the latter case, the trap has a restoring contribution, while the dark solitons (the kinks) of the NLS repel each other. Here, the situation is reversed: namely, as we saw in the previous section, the trap tends to expel the kink, while (as is known from the study of kink-antikink interactions in the homogeneous Ω = 0 case) the kink and antikink attract each other [10,2], i.e., the role of both contributions is reversed. As a result, the equilibrium point instead of being a center as it is in the case of NLS, here, again, it is a saddle. Indeed, if the kink pair is separated by a small distance, then the interaction prevails and forces the kinks to approach each other. On the other hand, a large separation leads to a dominant effect due to the trap, and the kink-antikink pair are pushed apart.
The static kink-antikink solution of (1) can be numerically approximated by where x0 is unique and dependent on Ω. We define that particular x0 as xcr, since it operates as a critical value separating between expulsion from and attraction of the kink and antikink to the center. For various values of Ω, the numerically calculated stationary kink-antikink solutions are shown in Figure 8. As usual, the right panel of the figure identifies the spectrum of the linearized operator about the static (in this case kink-antikink) solution. Here, we get two pairs of real eigenvalues, indicating instability. The pattern of the spectral problem is thus transparent: for one kink, there is a single unstable mode, for two kinks, two such modes and so on. Another way to state this is that all the negative energy eigenstates of [31] pertaining to the normal modes of vibration of the kinks in the real-field version of the model lead to unstable eigendirections.

Numerical Results (PDE)
We can obtain moving kink-antikink solutions by picking the initial conditions u(x, 0) = uΩ(x)(tanh(γ(x + x0)) − tanh(γ(x − x0)) − 1); ut(x, 0) = −vγuΩ(x)(sech 2 (γ(x + x0)) + sech 2 (γ(x − x0))); where γ = 1/ √ 1 − v 2 , ±x0 are the initial positions and ±v are the velocities of the moving kink and antikink. Once again, we utilize the (no longer representing an invariance) Lorentz boost and find it to be an efficient method of producing moving kinks-antikinks. Fig. 9 shows that when the velocity of the kink and antikink is initially zero, then for any x0 < xcr ≈ 1.6907 (for Ω = 0.15), the kink-antikink system forms a bound state. For any x0 > xcr, the system exhibits expulsion. This confirms our spectral understanding of the relevant equilibrium representing a saddle point. Now  the unstable manifold consists of two eigendirections, namely the in-phase and outof-phase motion of the pair, although here only the out-of-phase mode is excited.
When Ω = 0, the kink-antikink system is relatively insensitive to changes in the initial position x0 (provided that the kinks are not too close). This, however, changes in the presence of the parabolic trapping term. For x0 sufficiently small, the kink and antikink approach each other and collapse into a bound state for small input velocities, even with zero initial velocity. For x0 < xcr and initial velocities sufficiently large we observe the typical behaviors that are familiar from the homogeneous case [2,10]: reflection, n-bounce windows (i.e., windows where the kinks escape each other's attraction upon n collisions), and "bion/breather" [26,2] formation (i.e., emergence of a bound state that oscillates in time and gradually decays). In fact, the types of resonance windows observed in the homogeneous case [2,10,13] persist in the inhomogeneous case and present a similar type of pattern. A visualization of this can be seen in Fig. 10. In this figure and in subsequent figures, we opt not to plot velocities for which four or more bounces occur, as the windows become exceedingly narrow for this many bounces and it is not practical to display them together with one-and two-bounce windows. Just as when Ω = 0, the regions in between resonance windows correspond to "bion" formation (numerically these are treated as events with many bounces, so they are omitted from the figures). For x0 > xcr, the presence of the saddle point of the energy once again imposes a barrier that the kinks need to overcome in order to interact. Thus, in this case, there is a critical velocity associated to xcr, which we denote by vcr, such that whenever x > xcr and v < vcr the kinks do not have enough kinetic energy to reach the origin, so they stop (at a turning point) before they reach it. Then they recede towards spatial infinity without ever having collided.
If both x > xcr and v > vcr, then the kinks will have enough kinetic energy to overcome the barrier, reach the origin and collide, and then we observe the typical phenomena associated with the φ 4 model, including most notably the existence of multi-bounce windows. These can be visualized in Figs. 10 and 16 with both small and large separation x0, respectively. The velocity interval data that was used to create these figures is contained in Tables 1-4 in Appendix A.

Collective Coordinate Approach (ODE)
To complement our PDE simulations, we formulate a collective coordinate model for a kink-antikink configuration. As before, our aim is to reduce the full PDE with infinitely many degrees of freedom to a simple model with two degrees of freedom in a way similar to what has been previously done in the homogeneous φ 4 case [2,10,13,14]. To account for the internal modes of the kink and antikink, we introduce a new unknown A(t), representing the amplitude of the internal mode. We use only one function for both modes since the collision picture is symmetric between kink and antikink, so we expect the amplitudes of the kink mode and antikink mode to coincide if the initial configuration is symmetric. Our corresponding ansatz in this case will be: where u0(x) = tanh(x), and χ1 is an eigenfunction corresponding to the smallest positive eigenfrequency that satisfies the linearization equation for a single kink.
The particular eigenfunction we chose is given in Figure 11 for three different Ω values.
Surprisingly, these modes are related chiefly to the background state uΩ rather than to the internal vibrational mode of the kink. More specifically, they are localized near xs. In the homogenous model, the only nonzero discrete eigenfrequency corresponds to the internal vibrational mode, and it is this mode that is used for the collective coordinate ansatz. For our model, the mode with the smallest eigenfrequency results in a better match of collective coordinate and PDE results than when the (analogue of the homogeneous case kink) internal mode is used. This may be due to the former mode incorporating more adequately the role of the potentialinduced background, although, admittedly, this issue warrants further theoretical investigation. As proposed in [3], we work on a reduced effective Lagrangian that captures the fundamental features. The derivations of the formulas are in Appendix B.
Here, overdots denote time derivatives. We solve these equations numerically by using the initial conditions A(0) = 0,Ȧ(0) = 0, X(0) = x0 andẊ(0) = vin where x0 is the initial half-distance between the kink and the antikink, and vin is the initially prescribed speed of the kink (and anti-kink). Note that this is the reduced system where the higher order terms of A(t) andẊ(t) are ignored. We make this assumption based on the previous work [3], that higher terms of A andẊ can be dropped because they are small. The expressions involving integrals which are X-dependent in Eqs. (15) are computed by numerical integration.
For Ω values not much larger than the value considered so far (of 0.15), the phenomenology is virtually identical to that of the previous section. All of the interactions we have come across so far, namely bions, multi-bounce interactions, and expulsion, are observed for kink-antikink systems as the initial velocity is varied. Figure 16 shows the presence of two-and three-bounce windows and a threshold reflection velocity for Ω = 0.2 and Ω = 0.3. However, the widths of existing bounce windows shrink drastically as the trap strength is increased. Increasing Ω from 0.2 to 0.3 as in the same figure shrinks the resonance windows by approximately a factor of 1 10 . Our collective coordinate model also exhibits these features -it is capable of capturing one-and two-bounce interactions, and exhibits the same sensitive dependence on the initial speed we observe when Ω is small for configurations with more bounces. Examples can be seen in Figure 17.
Thus far we have confined ourselves to the case Ω ≪ 1. However, the behaviors for larger Ω are interesting in their own right and can differ significantly from those observed for weak trapping. Intuitively, we expect that a large value of Ω will constrain the dynamics severely (as it only allows the kinks to move within a restricted As Ω increases, existing multi-bounce windows become significantly narrower (note the change in the velocity increment in the rightmost figure), and many disappear completely. TF region). Our numerical simulations suggest that two-bounce and three-bounce still exist there, yet they are extremely narrow and shrink rapidly as Ω increases. Figure 18 shows some prototypical examples of multi-bounce cases for Ω = 0.4.

Conclusions and Future Challenges
In this work we have taken the well-known φ 4 model in a novel direction: motivated by considerations in atomic BECs where a complex variant of the model is examined in the presence of an external potential, we have explored how the standard φ 4 phenomenology is affected by a parabolic trapping. We have found that the results can be fundamentally different in comparison to the standard homogeneous (untrapped) case. In particular, even the single kink state is dynamically unstable, representing a saddle point in configuration space and leading to the expulsion of the kink from the trap center of symmetry. Moreover, in the presence of multiple kinks and anti-kinks, the number of coherent structures determines the number of unstable eigendirections (and the dimension of the corresponding unstable manifold). It is interesting to also observe the reverse role of the potential and of the kink-antikink interaction in comparison with the atomic BEC case. These features generally favor the expulsion of the kinks from the system unless the kinks are sufficiently close to each other that their interaction dominates. In that setting, n-bounce windows may still be present, although these windows were found to shrink radically as the potential parameter Ω is increased. Lastly, we have developed a collective coordinate low-dimensional ODE approach, under relevant approximations (such as the truncation of the higher power terms) which yields, especially for low numbers of bounces, reasonably consonant results between the PDE field theory and the ODEs. On the other hand, when presenting a more "complete" variant of the equations, issues similar to those encountered by [17,18] in the homogeneous variant of the model emerge. Indeed, the prospect of resolving these collective coordinate issues and improving the ODE models in this inhomogeneous φ 4 variant (and even for the homogeneous case within the model) is an important avenue for possible future work. From a mathematical standpoint, understanding better the stability properties of the model (e.g. proving rigorously the eigenvalue estimate in the case of N kinks, as likely possessing N unstable eigendirections) could be an interesting direction for further work. Also, extending inhomogeneous model considerations to other models of recent interest such as the φ 6 [15,17,19,20,32], and φ 8 models [23,33,34] would be valuable in understanding the genericity (or not) of the conclusions drawn in this study. Some of these directions are currently under consideration and will be reported in future publications.

A Tables for n− bounce windows
We list the velocity intervals [v1, v2] that result in multi-bounce windows when the initial velocity vin is picked within these intervals. In Tables 1-2, we list those intervals for a small separation (x0 = 1.4) and a big separation (x0 = 7) with Ω = 0.15. In Tables 3-4 Table 1: n-bounce windows for x 0 = 1.4. One-bounce windows occur for v in > 0.25845.  Table 2: n-bounce windows for x 0 = 7. One-bounce windows occur for v in > 0.74414.  Table 3: n-bounce windows for x 0 = 2, Ω = 0.2. One-bounce windows occur for v in > 0.28942.  Table 4: n-bounce windows for x 0 = 2, Ω = 0.3. One-bounce windows occur for v in > 0.36317.
Substituting (16) into Eq. (10) gives To handle the V (u) terms we first write u = ua + u b where ua = uΩ(φ+ + φ− − 1) and u b = A(χ+ + χ−). Then using the Taylor series expansion, we get The corresponding reduced Lagrangian (ignoring higher order terms) that is used in our simulations is given by Eq. (14). Applying Euler-Lagrange equations, we obtain Eq. (15).
We will list the formulae of the coefficients below. Note that they are all functions of X(t). Since χ1 is not known explicitly, the coefficients presented are in the integral form. These coefficients are calculated numerically.