Long-range interactions of kinks

We present a computational analysis of the long-range interactions of solitary waves in higher-order field theories. Our vehicle of choice is the $\varphi^8$ field theory, although we explore similar issues in example $\varphi^{10}$ and $\varphi^{12}$ models. In particular, we discuss the fundamental differences between the latter higher-order models and the standard $\varphi^4$ model. Upon establishing the power-law asymptotics of the model's solutions' approach towards one of the steady states, we make the case that such asymptotics require particular care in setting up multi-soliton initial conditions. A naive implementation of additive or multiplicative ans\"atze gives rise to highly pronounced radiation effects and eventually leads to the illusion of a repulsive interaction between a kink and an antikink in such higher-order field theories. We propose and compare several methods for how to"distill"the initial data into suitable ans\"atze, and we show how these approaches capture the attractive nature of interactions between the topological solitons in the presence of power-law tails (long-range interactions). This development paves the way for a systematic examination of solitary wave interactions in higher-order field theories and raises some intriguing questions regarding potential experimental observations of such interactions. As an Appendix, we present an analysis of kink-antikink interactions in the example models via the method of collective coordinates.


I. INTRODUCTION
Field-theoretic models with polynomial potentials are of great interest in many areas of modern theoretical physics: from cosmology [1, 2] to condensed matter [3,4]. For example, the scalar ϕ 4 model with two minima of the potential is widely used to model spontaneous symmetry breaking. Besides that, the quartic potential arises in the phenomenological Ginzburg-Landau model of superconductivity [5,6]; see also Refs. [4,7] for an overview of different areas of application. In this setting, the dynamics of the complex scalar field of Cooper pairs is described by a polynomial self-interaction of the fourth degree. Models with polynomial potentials of higher degrees are commonly used, e.g. to model consecutive phase transitions [8], which arise in material science [9], as recently summarized in [10].
It has also been shown that scalar field theories can describe distinct quantum mechanical problems (including supersymmetric ones) [11], leading to new applications of field theories of the type ϕ 2n . Another possibility is to consider scalar ϕ 2n field theories as Lane-Emden truncations of a periodic potential, which then leads to applications of these theories as toy models for dark matter halos [12]. Field-theoretic models with polynomial potentials that can exhibit topological solutions (kinks) are also important in cosmological applications of been studied systematically. Nevertheless, the work has been started [8,10,[32][33][34][35][36][37][38][39]. In particular, exact (but implicit) solutions for the kink shapes of various ϕ 8 , ϕ 10 and ϕ 12 field theories have been obtained [8], the excitation spectra of the ϕ 8 kinks have been studied, resonance phenomena in the kink-antikink scattering have been found, and their relation to the kinks' vibrational excitations has been discussed [35,36,38,39]. As described below, we use the purely theoretical foundation established in [8] to develop a novel understanding of long-range kink interactions, re-assessing in significant detail the interaction picture first put forth in [38,39]. However, issues of vibrational modes in higher-order field theory [35,36] remain beyond the distinct scope of the present work.
Until recently, the dynamics of kink-(anti)kink interactions had been studied only for kinks with exponential tail asymptotics. At the same time, it is easy to show that, for certain potentials of sixth or higher degree, there exist kinks that exhibit power-law asymptotics of either or both tails connecting two distinct equilibria. Although conditions for algebraic soliton solutions of the nonlinear Schrödinger, Korteweg-de Vries and related integrable models have long been known [40,41], the case of kinks with algebraic tails in non-integrable ϕ 2n field theories remain less explored in comparison.
In this paper, building upon the preliminary report of some of the present authors [37], we systematize restrictions on the potential and obtain general formulae for the kinks' tail asymptotics (see also Refs. [8,32,34,42]), including the conditions for power-law tail asymptotics in a sextic potential. The existence of kinks with power-law tails is of particular interest because such tails lead to long-range interaction between a kink and an antikink. Studying such long-range interactions at the effective "particle" level is a topic of significant current interest. Thus, we undertake and exploration of the interaction chiefly via direct numerical simulations of the relevant field equation. An Appendix complements our computations with analytical considerations based on the variational technique known as the collective coordinate approximation, which is widely used in various field-theoretic problems, see, e.g. [7, 25-27, 30, 43-50]. We note, in passing, that other approaches have been used to interrogate long-range interactions, such as evaluating the field's potential energy at the center of mass of two superimposed (anti)kinks [51,52]. Also, the effect of minima of the potential and their relative depth (including inflection points in-between) on the kink asymptotics and their mutual forces was considered in [53]. An alternative method for identifying the interactions of solitary waves is the so-called Manton's method that has been widely used in solitonic equations [2, [54][55][56][57]. A very recent attempt to utilize this and other methods to infer a power-law dependence of the force on the kink-antikink separation, in the case of the ϕ 8 model, has just been posted in [58,59]. The result in [58] corroborates our previous observation in [37] that, for an example eight-order potential with three degenerate minima and one-sided power-law tail asymptotics of the kink, a kink and an antikink attract each other with a force proportional to their separation to the −4 power. In this work, we provide numerical evidence for this type of power-law attractive force. However, a more conclusive theoretical investigation (including discussion of the pre-factor on this power law) is deferred to future work; see also the relevant discussion in Sec. V.
Our presentation is structured as follows. In Sec. II, we begin by providing the kink asymptotics in higher-order ϕ 2n field theories (and specifically for our principal ϕ 8 example). Subsequently, in Sec. III, we delve into our numerical considerations, starting with how numerical experiments of kink-antikink collisions are set up, explaining the difficulties of such a setup for higher-order field theories, and proposing a corresponding methodology for handling such difficulties in the ϕ 8 model. The latter are complemented by parallel considerations of example ϕ 10 and ϕ 12 field theories. A discussion of the force of interactions between a kink and an antikink via power-law tails is presented in Sec. IV. We then summarize our work and propose some (among the many intriguing) questions for future work in Sec. V. In the Appendix, we explain how to perform a calculation of long-range interactions based on the method of collective coordinates.

II. ANALYTICAL CONSIDERATIONS
A. Power-law asymptotics of kinks Consider a real scalar field ϕ(x, t) in (1 + 1) dimensional space-time, with its dynamics determined by the Lagrangian density where V (ϕ) is a potential that defines the self-interaction of the field ϕ. The energy functional corresponding to the Lagrangian (1) is Assume that the potential V (ϕ) is a non-negative polynomial of even degree (denoted in short-hand as "ϕ 2n ") having two or more minimaφ 1 ,φ 2 , . . .,φ n of equal depth V (φ 1 ) = V (φ 2 ) = · · · = V (φ n ) = 0. Consider two adjacent minimaφ i andφ i+1 . Let ϕ K (x) be a kink (see, e.g. [2, Ch. 5]) interpolating between these minima, i.e., According to conventional notation, we can say that this kink belongs to the topological The Euler-Lagrange equation of motion, which is the condition for extremizing the action generated by the Lagrangian density (1), is The kink shape function ϕ = ϕ K (x) is a time-independent solution of Eq. (4), and it can be shown that it satisfies a first-order ordinary differential equation: Taking into account that 2V (ϕ) dx = dϕ, the energy (rest mass) of a static BPS-saturated field configuration is Now, we formulate general conditions that must be satisfied in order for the model to have kinks with power-law tail asymptotics. We also give general formulae for such asymptotics.
Let us turn our attention to the potential V (ϕ). Letφ i andφ i+1 be zeros of the function V (ϕ) of orders k i and k i+1 , respectively. Notice that k i and k i+1 must be even. We assumē ϕ i <φ i+1 , for definiteness. Then, the potential V (ϕ) can be written as Inserting Eq. (8) into Eq. (5) and recalling thatφ i < ϕ < ϕ i+1 , we obtain after integrating: To find asymptotics of the kink ϕ (φ i ,φ i+1 ) (x) as x → −∞, we use the fact that ϕ →φ i in this limit. Then, there are slowly varying factors within the integrand in the right-hand side of Eq. (9) at ϕ →φ i . These factors can be (approximately) taken out from the integral. As a result, we obtain an asymptotic equality: Integration of the right-hand side gives power-law dependence, if k i > 2. Taking this observation into account, we obtain the asymptotics of the kink as x → −∞: Similarly, for case of k i+1 > 2, we obtain asymptotics of the kink as x → +∞: Below, we will use Eqs. (11) and (12) to find power-law asymptotics of kinks of a ϕ 8 model.
Note that exponential asymptotics could be found by integrating Eq. (10) for the case of k i = 2, but a further refinement of our approximations is needed to capture the prefactor of the exponential.
Inserting the potential (13) into Eq. (5) and requiring that 0 < |ϕ| < 1, we have Integrating this equation, we obtain (implicitly) the two kinks belonging to the topological sectors (−1, 0) and (0, 1): The corresponding antikinks can be obtained from (15) by the transformation x → −x: It can be directly inferred from Eq. (16) that the asymptotics of the kink ϕ (0,1) (x) as x → −∞ is of power-law type; likewise, for the asymptotics of the kink ϕ (−1,0) (x) as x → +∞. Using (7) we obtain the energy of the static kink (15): Of course, this energy is the same for all possible kinks of the model with the potential (13), see also Ref. [8].
Now, we derive the asymptotic behavior of the kinks ϕ (−1,0) (x) and ϕ (0,1) (x). We can do this in two ways: firstly, by using the approach from Section II A, and secondly, by expanding Eq. (15) in a Taylor series around the appropriate limiting point(s). Consider the kink ϕ (−1,0) (x). According to the notation of Section II A, for this kink: Equations (11) and (12) are approximately applicable only if k i or k i+1 is greater than 2, respectively. So, Eq. (12) gives power-law asymptotics of the kink: This asymptotic expression can also be obtained from Eq. (15) as shown in [8]. Indeed, the tail asymptotics emerge rather straightforwardly from the implicit kink solution as the logarithmic term becomes a small correction in the limit. At the same time, the asymptotics at x → −∞ is exponential, and it can be obtained from Eq. (15): For the kink ϕ (0,1) (x), we have:φ Similarly to the previous case, from Eq. (11) we obtain power-law asymptotics at x → −∞: and the exponential asymptotics at x → +∞ can be obtained from Eq. (15): To summarize: in this Section, we discussed the kinks in topological sectors (−1, 0) and (0, 1) for the featured ϕ 8 model with the potential (13). In particular, we highlighted that both of these kinks have one power-law and one exponential asymptotic decay to their equilibrium background states (vacua) at |x| → ∞. In what follows, we use the kink ϕ (−1,0) (x) and its corresponding antikink to study their their interaction via their power-law tail asymptotics. We will argue that power-law asymptotics lead to long-range interaction: a kink and an antikink "feel" each other at very large separations. This situation is quite different from the case of exponential tail asymptotics (for example, if ϕ (−1,0) (x) and its antikink were reversed in their initial configuration, or as in the classical ϕ 4 field theory), in which case the kink-antikink interaction is exponentially decaying.

III. DIRECT NUMERICAL SIMULATION OF COLLISIONS
The non-integrability of models such as the above-mentioned ϕ 2n field theories suggests that exact multi-soliton solutions are not available in these systems. Nevertheless, as it is well-known from numerous prior works [7, 21-23, 26, 27, 44, 60-63], the study of kinkantikink collisions is both particularly interesting in its own right and potentially presents a very rich phenomenology. Among the many phenomena observed in such collisions are multibounce windows, the fractal structure thereof, the role of the presence (or even absence [23]) of internal vibration modes, the ability to describe such phenomena via collective coordinate methods (or complications [26,27] thereof), and many others. Of critical importance to all of the above is the computational feasibility of interrogating collision phenomena via direct numerical simulations of the (1 + 1)-dimensional field theory. Indeed, the main message of the present work is that the standard approaches for setting up such simulations do not work for higher-order field theories, such as the ones considered herein. In light of this observation, we begin by discussion the well-known ϕ 4 field theory. Subsequently, we compare/contrast it with the specific ϕ 8 model from Sec. II B, and finally we present some possibilities for handling the complications due to long-range interactions of kinks.
Prior to discussing the numerical results, it is relevant to briefly describe the methods used to obtain them. We discretize the governing equation of motion (4) on the spatial 200,200] with the increment ∆x = 0.2 (which fixes the number of Fourier modes used). We use a Fourier-based pseudospectral differentiation matrix D 2 as in [65] to approximate ∂ 2 ϕ/∂x 2 as D 2 ϕ subject to periodic boundary conditions. This step turns the PDE into a semi-discrete system of second-order-in-time ordinary differential equations (ODEs). These are trivially rewritten as a first-order system of ODEs and integrated forward in time using MATLAB's ode45 differential equations solver with adaptive time stepping and error control.
A. The standard example: ϕ 4 field theory The classical ϕ 4 field theory is determined by the potential V (ϕ) = 1 2 (ϕ 2 − 1) 2 (see, e.g., [2, 19,44]). The stationary kink solution of this model in the (−1, 1) topological sector, i.e., the solution of the BPS equation (5), is ϕ (−1,1) (x) = tanh(x). By Lorentz boosting the stationary solution, we obtain a traveling kink solution: for any kink velocity v such that −1 < v < 1. A traveling antikink solution is given by −ϕ v (x, t) for this model. A kink moving to the right with velocity v, shifted to the left by the amount x 0 , is then given by , and an antikink moving to the left with opposite velocity, shifted to the right by the amount x 0 , is likewise given by Therefore, the function represents a waveform consisting of a kink and an antikink with initial separation 2x 0 . The kink and antikink in this pair have equal and opposite velocities ±v. Furthermore, ϕ(x, t) given by Eq. (25) (4) to approximately standard machine precision, i.e., on the order of 10 −16 . Specifically by "satisfies the PDE" we mean that the residual, as measured by the value of is suitably small. Therefore, it is reasonable to use Eq. (25) to generate initial conditions for kink-antikink collisions with and where primes denote differentiation with respect to the function's argument. However, for separations 2x 0 20 units, the value of maxAbsP de (evaluated from the numerical solution of the PDE starting from the initial conditions in Eqs. (27) and (28) with v = 0) decreases exponentially with x 0 ; for stationary solutions, its magnitude is on the order of 10 −7 at a separation of 2x 0 = 10, and on the order of 10 −3 at a separation of 2x 0 = 5.
It is relevant to mention here that the ansatz in Eqs. (27) and (28) is suitable not only for direct numerical simulations, but also for collective coordinate approximations of the PDE dynamics [7, 25-27, 30, 43-50]. In particular, one can use Eqs. (27) and (28) with x 0 − vt → X(t) as a new variable (the collective coordinate) that determines the dynamic location of the kink's center. In some of the latter references, more elaborate ansätze involving also a coordinate characterizing the kink's internal (vibration) mode were considered. However, these considerations are beyond the scope of the present study. The principal features of a collective coordinate approach to kink-antikink interactions in the higher-order field theories of interest herein are presented in the Appendix, for completeness.
B. The present case: ϕ 8 field theory Now, consider the governing PDE (4) with the potential V (ϕ) = ϕ 4 (ϕ 2 − 1) 2 as in Sec. II B above. We can study a kink-antikink collision interaction by adapting the ansatz from Eq. (25) as follows: where ϕ (−1,0) and ϕ (0,−1) are given implicitly by Eq. (15). As in the ϕ 4 example above, we can use Eq. (29) to generate initial conditions for a collision simulation. However, because of the power-law tails of the kink and antikink, we are presented with a problem. In To this end, consider the PDE residual maxAbsP de defined in Eq. (26). Now, we substitute Eq. (29) into Eq. (26) to determine whether this ansatz provides a suitable approximate kink-antikink solution to the PDE (4). Evaluating maxAbsP de numerically (keeping in mind that ∂ 2 ϕ/∂t 2 = 0 for the stationary solution with v = 0), using the pseudospectral differentiation matrix D 2 to approximate ∂ 2 /∂x 2 as mentioned above, for x 0 = 5, 10 and 20, we obtain the values 0.72, 0.32 and 0.15 respectively, all on the order of 10 −1 (see Table I  kink-antikink pair at a separation of 10) we found that maxAbsP de was on the order of 10 −7 .
Thus, in contrast to the ϕ 4 case, even for a separation as large as 40, we find that the linear superposition (i.e., sum) ansatz in Eq. (29) does not provide an approximate solution to the ϕ 8 equation of motion in a quantitative sense. We thus warn the numerous practitioners of such numerical computations regarding the substantial obstacles to using the classical sum ansätze to study collisional dynamics of kinks and their interactions numerically.
Next, consider what happens when we use Eq. (29) to create initial conditions for a prototypical kink-antikink collision simulation. We restrict ourselves to the case in which the kink and antikink are initially stationary, i.e., v = 0; we do this to avoid the complication(s) of how an initial kinetic energy may affect the dynamics. In Fig. 2   Earlier work [39] has suggested that a repulsive force might exists between the example kink and antikink considered above. Yet, we argue that this apparent repulsion is a result of the ansatz selected in Eq. (29) being a poor quantitative approximation of a kink-antikink solution. Recall the undershoot below ϕ = −1 near x = ±20 in Fig. 1. We can think of this undershoot as providing a kind of initial "spring board," which pushes the "points" just above the spring board upward, and consequently causes the kink to move to the left and the antikink to move to the right. Furthermore, this upward motion creates disturbances at x ≈ ±20 that move upwards and then along the top of the kink-antikink combination until they meet at x = 0. One can see these effects clearly in the contour plot in Fig. 2(a). One can also see that, after meeting at x = 0, these disturbances are trapped between the centers of the kink and antikink, which they reach again just before t = 50. From the plot of the kink velocity in Fig. 2(b), it is clear that just at this time the arrival of the disturbance gives another boost to the velocity of the kink in the negative direction. Again, the disturbances are reflected back towards x = 0 and out again to the centers of the kink and antikink for another boost to the velocities, pushing them apart faster. This phenomenology holds for other values of the initial half-separation x 0 , down to x 0 = 2.5. Thus, we have accounted for the apparent repulsive force, in the long-range interaction between a ϕ 8 kink and antikink, by showing that this "force" is simply the result of initial conditions that have been derived from an inaccurate sum ansatz. More specifically, we have quantified how this ansatz does not lead to a sufficiently accurate description of the motion of a kink-antikink pair. An additional (albeit weaker) effect along this vein consists of the radiative wavepackets, emitted from each kink, that affect both of them in the process. As a remedy, we propose to improve the initial conditions introduced in the previous section by employing Eq. (29) as the initializer for a weighted nonlinear least-squares minimization of the objective function where · 2 is the usual Euclidean norm, and C is an empirical constant. Then, we can use the minimizer ϕ min (x) of I as the initial condition for a direct numerical simulation of kink-antikink collisions, ensuring that our initial condition quantitatively satisfies the PDE to some preset accuracy. We take C = 50, which is sufficient to keep the initial kink and antikink locations nearly fixed at ±x 0 during the minimization process. The optimization problem is solved using MATLAB's optimization toolkit, specifically via the lsqnonlin subroutine.
In Fig. 3, for x 0 = 20, we compare the sum ansatz from Eq.  maxAbsP de for the minimized and (non-minimized) sum ansätze. Specifically, the maximum value of maxAbsP de when taking ϕ = ϕ min (x) is several orders of magnitude smaller than when using the sum ansatz from Eq. (29). Thus, we conjecture that the initial conditions generated from ϕ min (x) will more accurately reflect the actual kink-antikink solution of this non-integrable field theory, at least considerably better than the non-minimized sum ansatz from Eq. (29).
To support our conjecture, Fig. 4 shows the result of a direct numerical simulation of the PDE (4), for our ϕ 8 model, using the minimized function ϕ min (x) to generate the initial conditions. As before, this simulation corresponds to a kink-antikink interaction with v = 0 because ϕ min (x) approximates a stationary solution of the PDE. As in Fig. 2, we show both a contour space-time plot and a velocity plot in Fig. 4. Now, we observe an attractive force between the kink and antikink. Also, there are no visibly discernible small disturbances moving back and forth between x = ±20. This example simulation, along with numerous other similar simulations for different x 0 values, suggest that we have eliminated the significant detrimental effects of the algebraic kink and antikink tails (and of radiation), and we can thus now observe the proper (effective) inter-particle interaction between the kink and antikink. Also, note that this interaction is more in line with what one would naturally expect from a PDE of the type in Eq. (4), in that the ∂ 2 ϕ/∂x 2 term tends to "pull points down" (pulling the kink and antikink together) when the function is concave down. Finally, the minimized ansatz ϕ min (x) serves not only as an initial condition for the PDE simulations, but also as an appropriate ansatz for the collective coordinates approach described in the Appendix.

D. Other possible ansätze
Besides finding ϕ such that I in Eq. (30) is minimized, there are other possible setups that could be used to generate initial conditions for direct numerical simulation of the PDE (4) (and possibly for collective coordinate approaches as well). For example, one can use a product (rather than a sum) ansatz. Such an ansatz, customized for our ϕ 8 model, is which we term the product ansatz.
Numerical simulations starting from Eq. (31) as an initial condition (again with v = 0), once again exhibit a repulsion initially (though, weaker than for the sum ansatz) for x 0 > x c where 6.2 < x c < 6.3. For x 0 < x c , attraction is found in the direct numerical simulations.  the product ansatz also leads to generic attraction, as we have come to expect, at this point, from the ϕ 8 field theory. All of these observations are illustrated in Fig. 7.
A third option is to treat the kink and antikink "completely separately," meaning to use the kink formula for x < 0 and the antikink formula for x ≥ 0. To accomplish such a feat, we define which we term the split-domain ansatz. Here, H(x) is the Heaviside unit-step function. Using this ansatz to generate the initial conditions for a PDE simulation, we plot the contours of ϕ and the kink velocity for x 0 = 6.2 and x 0 = 6.3 in Fig. 8. Contrary to what was the case for the product ansatz, we observe attraction for both x 0 values. It is perhaps natural to expect that the split-domain ansatz is the most accurate unminimized one (i.e., among the more standard ones that have not been "optimized" via our proposed minimization procedure), but it is expected to be limited in accuracy in the vicinity of x = 0 due to the derivative discontinuity introduced in Eq. (32). This observation is substantiated by the kink-antikink dynamics shown in Fig. 8.
In this case, the ansatz is continuous at x = 0 but its first derivative is not. The minimized version for the split-domain ansatz is quite similar to the non-minimized version; the "point" created where the kink and antikink meet at x = 0 (due to the discontinuity in the derivative) is smoothed by the minimization procedure, but the two look rather similar otherwise. Figure 9 shows how the minimized version of the split-domain ansatz differs from its non-minimized counterpart. The dynamics of the split-domain ansatz comes closest to the minimized case in that it generically leads to attraction of the kink and antikink.
One other property of the various ansätze that is worth mentioning is the relative smoothness of the velocity graphs. One sees significant oscillations in the velocity of the center of the kink (or antikink) for all of the non-minimized ansätze (see, e.g., Fig. 2). The minimized versions of all three ansätze, on the other hand, show a steadily increasing velocity function (see, e.g., Fig. 4).
In Table II, we show a comparison of the minimized versus non-minimized values of the PDE residual, maxAbsP de, for the product and split-domain ansätze. Note that for the non-minimized split domain case, ∂ 2 ϕ/∂x 2 is not defined due to the discontinuity in the derivative at x = 0, and hence the value of maxAbsP de is listed as NA ("not available").
Another reason as to why we have elected not to provide this value is that the ansatz was constructed from the (numerically evaluated) exact solution of the BPS equation (no minimization) for each x > 0 and x < 0, which already satisfy maxAbsP de = 0 numerically. Figure 10 shows the equivalent plots of those in Fig. 2 (non-minimized sum ansatz) and A plot giving a sense of how the initial conditions for the ϕ 8 model compare for the different ansätze is shown in Fig. 11 for x 0 = 4.5 and v = 0. We observe that all minimized ansätze lie between the sum/product and the split-domain ansatz. Interestingly, the results of the different minimization procedures are close to each other functionally, and the differences between them are difficult to detect without zooming in.
To summarize: our results indicate that the correct interpretation of the nature of the pairwise kink-antikink interaction is that the kink and antikink attract each other. We  II. maxAbsP de for minimized ("min") and non-minimized ("non") product ("prod") and split-domain ("split") ansätze for the example ϕ 8 model.
proposed the minimization procedure in Sec. III C as a way to "distill" the initial data and, thus, observe the genuine interaction dynamics of the kink and antikink without the detrimental side effects of the undershoot caused by their tails, as well as the radiation caused by the inexact initial conditions. While it is impossible to push the objective functional I from Eq. (30) to zero exactly (due to the absence of multi-soliton solutions for such a nonintegrable model), the minimization of I brings the initial ϕ field as close as possible to a distilled configuration involving the superposition of a kink and an antikink. On the other hand, in the absence of access to such a minimization procedure, our recommendation is to use the split-domain ansatz from Eq. (32) directly, as it is the one that bears the least spurious byproducts among the "standard" multi-soliton ansätze, even though it introduces a derivative discontinuity at x = 0. We find similar behaviors when considering ansätze for the corresponding ϕ 10 and ϕ 12 field theories with three degenerate minima. In particular, we considered the models represented x → −∞, but approaches 0 as k/x 1/2 (for some k constant) when x → +∞. Similarly, a kink of the model ϕ 12 potential above approaches −1 exponentially as x → −∞, but approaches 0 as k/x 1/3 (for some k constant) when x → +∞. Thus, these higher-order field theory models possess solutions with "fatter" tails. Table III shows the maxAbsP de residuals for the ϕ 10 model in a way that parallel Tables   I and II for the ϕ 8 case. Table IV shows the equivalent results for the ϕ 12 model. We observe  III. maxAbsP de for minimized ("min") and non-minimized ("non") sum, product ("prod") and split-domain ("split") ansätze applied to the example ϕ 10 model. similar trends for these models to what we saw in the ϕ 8 case; once again, the minimization procedure significantly improves the quantitative agreement between an initial condition ansatz profile and a hypothetical one that exactly satisfies the PDE (4).
for these cases, the ansatz leads to solutions that show oscillations in the range ϕ = −1 to ϕ = 1 rather than topological solitons connecting ϕ = −1 to ϕ = 0.
A way to explain this strange result is to consider the fact that for such fatter-tail cases, as the ones arising from the example ϕ 10 and ϕ 12 field theories considered herein, the undershoot of the kinks is so substantial that not only is the ϕ = 0 fixed point not reached between the kinks, but also neither is the asymptotic value of ϕ = −1 as |x| ≫ 1. In other words, the (non-minimized) sum ansatz provides an extremely poor initial conditions for the fatter-tail cases.
We can obtain some further insight into this quantitative disagreement by considering the graphs of the sum ansatz for the three cases; see  Table IV), and to a lesser extent for x 0 = 5, 10 and 20 for the ϕ 10 model (recall Table III). In all of these cases, the notions of "attracting kinks" and "repelling kinks" are no longer meaningful.

OF SEPARATION DISTANCE
Lastly, we can illustrate the attractive nature of the force between a kink and antikink in two additional ways. In Fig. 13, we show the potential energy E[ϕ] as defined in Eq.
Also, we have calculated the acceleration of the left kink (a proxy for the kink-antikink force of interaction), from the previously computed kink velocity v from the PDE evolution, as function of x 0 for a minimized split-domain initial condition ansatz and zero initial velocity. In this case, the acceleration is quite steady for a short period of time. Six data points were collected for the ϕ 10 and ϕ 12 models, and seven for the ϕ 8 model, with x 0 ∈ [20, 300], and a power-law model ax −b 0 was fit to the data. In all cases, an excellent fit was obtained. Specifically, we find b = 3.998 ± 0.002 for the chosen ϕ 8 model, b = 3.067 ± 0.019 for the chosen ϕ 10 model, and b = 2.764 ± 0.025 for the chosen ϕ 12 model, all within the 95% confidence interval for the fit. Figure 14 shows the simulation data and fits on a log-log plot, for all three cases.
We note that this scaling (i.e., b ≈ 4 for the example ϕ 8 model) of the acceleration (thus, the force of interaction) with the half-separation is in line with the theoretical prediction for the −4 power-law decay of the force for the same ϕ 8 model in [37] and, more recently, in [58,59]. A systematic study of this interaction force and its dependence on the kink-antikink separation, for arbitrary power-law tails, is the subject of future work [66], as it is a topic of interest in its own right; we would digress from the main theme of the present study if  we were to pursue it here.

V. CONCLUSIONS AND FUTURE WORK
In the present work, we have systematically interrogated the dynamics of kink-antikink interactions in higher-order polynomial field theoretic models (of eight degree and higher) with degenerate minima. The specific feature of these models that we have sought to capture is the presence of long-range interactions via kink tail asymptotics that do not decay exponentially but rather decay algebraically. Although a ϕ 8 model was used as our featured example, we also demonstrated that our discussion of how to properly generate initial conditions for direct numerical simulations of kink-antikink interactions applies to ϕ 10 and ϕ 12 models with "fatter tails." Our main finding was that, for all of these higher-order field theories, the standard sum ansätze (of a kink plus an antikink, spaced some distance apart) for direct numerical simulation of collision are problematic. These ansätze exhibit a significant undershoot around the two kinks in the combined profile, which leads to considerable radiation in the numerical solution for t > 0. These unwanted effects, in turn, are responsible for the apparent observation of unwarranted features such as repulsive dynamics (for a sum ansatz) or transitions between repulsive and attractive dynamics (for a product ansatz). We have argued that, among the simpler ansätze available, the one leading to the most realistic kink-antikink interaction dynamics (i.e., the final results are not contaminated by the details of the initial conditions for a simulation) is the one we have termed the split-domain ansatz.
Moreover, we have argued towards the usefulness of a suitable minimization procedure that "distills" a given ansatz further by making a closer match to a stationary kink-antikink solution of the problem. The minimization procedure reduces the undershoots in the combined profile and, thus, reduces radiation wavepackets (as well as their side effects) in simulations.
Once an initial condition was thus suitably prepared, we observed attraction between a kink and an antikink in all of the higher-order field theories (i.e., our prototypical ϕ 8 example, as well as ϕ 10 and ϕ 12 models), much like in the classical ϕ 4 field theory. This type of improvement in the kink-antikink state construction allowed us also to unambiguously obtain the power law nature of the kink-antikink interaction force and how its exponent varies among the different higher-order field theories examines.
Naturally, this study opens up numerous avenues for future work on the interactions of topological solitons in higher-order field theories. The most canonical extension of this work concerns the outcome of collisional events for different initial speeds of the kink and antikink (in this work, we took v = 0 in all of our examples), and across our proposed variety of initial condition ansätze. Such an exploration and a corresponding systematic study will be reported elsewhere in the future. Another open question is: how much of the abovedescribed interaction picture can be captured through a semi-analytical approximation such as the method of collective coordinates (CC)? A first attempt is given in the Appendix that follows, yet as can be seen there, it is somewhat limited in its ability to capture in detail the kink-antikink interactions. Moreover, at the present stage, the CC model is lacking the inclusion of the internal vibration mode of the kink; incorporating the latter appears to be extremely cumbersome in the present setup. Lastly, another important question is: how much of the above-described phenomenology can be captured in an experiment? We are not immediately aware of experiments involving higher-order field theories. However, for complex variants of the ϕ 4 field theory, such as nonlinear Schrödinger models, kinks can be  [2] N. Manton In this Appendix, we apply the method of collective coordinates (CC) using the minimized initial condition ansätze, which we introduced in the main text above to analyze the kinkantikink interactions numerically. To this end, recall that the Lagrangian for our neutral scalar field theories is where the Lagrangian density is given in Eq.
(1) and the potential is given in Eq. (13) for the chosen ϕ 8 model. Now, for all minimized ansätze, we can reduce the PDE (4) to a Hamiltonian dynamical system with one degree of freedom as follows. First, we obtain an effective Lagrangian by evaluating Eq. (33) using the ansatz for ϕ, having identified "x − vt" as the collective coordinate X(t). The manipulation is formally denoted as where different ansätze yield different functions b 0 (X) and b 1 (X). In the following subsections, we will present the formulae for these coefficients for each ansatz. The Euler-Lagrange equation rendering the functional L eff in Eq. (34) stationary is The resulting dynamical evolution equation, written as a first-order system, iṡ We solve this first-order ODE system (36), subject to the initial conditions X(0) = x 0 , Y (0) = 0 (corresponding to v = 0). As before, x 0 is the initial half-separation between the kink and the antikink, and it is assumed that the initial speed of the kink and antikink is zero. For integration of the system, we use MATLAB's ode45 differential equations solver with adaptive time stepping and error control. In (b), solid curves represent K a 1 (x), and dashed curves represent K a 2 (x). In both panels, x 0 is varied from 5 to 10 to 20 (darker color curves to lighter color curves, respectively).

CC method for the improved (minimized) sum ansatz
Let f a (x) be the function obtained by using the minimization of sum ansatz (corresponding to light blue curve in Fig. 11 for x 0 = 4.5) for the ϕ 8 model, i.e., Eq. (29). Then, assume a colliding kink-antikink scenario with the following field configuration ϕ(x, t) = K a 1 (x + X(t) − x 0 ) + K a 2 (x − X(t) + x 0 ) − f a (0), where X(t) is the half-distance between the kink and antikink and Observe that when X(0) = x 0 , we have ϕ(x, 0) = f a (x), which implies that the initial conditions for the ϕ 8 equation of motion (4) (i.e., "PDE model") and the CC approach (i.e., "ODE model") match. Figure 15 presents such functions for x 0 = 5, 10, 20.
When the initial velocity is zero, both the ODE model and the PDE model predict attraction for x 0 < x c , where 6 < x c < 7. The ODE model results agree better with the PDE model results as x 0 becomes smaller, as shwon in Fig. 16.