Characteristic approach to the soliton resolution

As a toy model for understanding the soliton resolution phenomenon we consider a characteristic initial boundary value problem for the 4$d$ equivariant Yang-Mills equation outside a ball. Our main objective is to illustrate the advantages of employing outgoing null (or asymptotically null) foliations in analyzing the relaxation processes due to the dispersal of energy by radiation. In particular, within this approach it is evident that the endstate of evolution must be non-radiative (meaning vanishing flux of energy at future null infinity). In our toy model such non-radiative configurations are given by a static solution (called the half-kink) plus an alternating chain of $N$ decoupled kinks and antikinks. We show numerically that the configurations $N=0$ (static half-kink) and $N=1$ (superposition of the static half-kink and the antikink which recedes to infinity) appear as generic attractors and we determine a codimension-one borderline between their basins of attraction. The rates of convergence to these attractors are analyzed in detail.


Introduction
According to the soliton resolution conjecture, global-in-time generic solutions of nonlinear dispersive wave equations resolve for t → ∞ into a superposition of decoupled nonlinear bound states (solitons) and radiation [1]. There are numerous physical manifestations of this phenomenon, ranging from the formation of solitons in optical fibers (modelled by a nonlinear Schrödinger equation) to the formation of stationary black holes in binary black hole mergers (modelled by the Einstein equation).
The past decade has seen a significant progress in mathematical understanding of the soliton resolution, especially for radial solutions of the energy critical nonlinear wave equations (see [2] for a survey and references therein). Notably, the soliton resolution was recently proved for the equivariant wave maps R 2+1 → S 2 [3,4] and the equivariant Yang-Mills (YM) equation in 4 + 1 dimensions [4] (both for the global and blowup solutions). These remarkable results are abstract in the sense that they enumerate all possible asymptotic scenarios but do not settle which scenarios are actually realized. As a toy model for more quantitative description of the soliton resolution phenomenon, in this paper we consider the equivariant YM equation in 4 + 1 dimensions where W (t, r) is the YM potential. As the spatial domain we take the exterior of the unit ball, i.e. r ≥ 1, and impose the Dirichlet condition on the boundary W (t, r = 1) = 0. The associated conserved energy is Finiteness of energy requires that |W (t, ∞)| = 1. Since the singular point r = 0 is outside the domain, it is easy to see that solutions starting from smooth, finite-energy initial data (W (0, r), W t (0, r)), which are compatible with the boundary condition, remain smooth for all times. Our goal is to describe their asymptotic behavior for t → ∞.
The key role in our analysis will be played by the half-kink which is the unique (modulo sign) static solution of equation (1) satisfying the boundary condition Q(1) = 0. The half-kink is a global minimizer of energy (see the Bogomolnyi inequality (20) below) and thereby a natural candidate for an attractor. Indeed, we will see that on any compact interval [1, R) every smooth solution W (t, r) converges to Q(r) or −Q(r) as t → ∞. However, for sufficiently large energies a nontrivial coherent structure can simultaneously develop in the asymptotic region (R, ∞). This behavior is intimately related to the energy criticality of the model and is absent in supercritical dimensions; cf. the soliton resolution for the equivariant wave maps exterior to a ball in 3 + 1 dimensions [5,6].
The paper is organized as follows. In section 2 we first recall from [4] the formulation of the soliton resolution conjecture for equation (1) in the whole space. Then we present an analogous conjecture in our model and support one special case by a heuristic argument based on the method of collective coordinates. In section 3 we reformulate the initial-boundary problem in terms of null foliations of constant retarded time and compactified spatial domain. Using this formulation, in section 4 we consider the late-time behavior (the quasinormal ringdown and the polynomial tail) for the linearized problem. Finally, in section 5 we present numerical evidence supporting the soliton resolution conjecture.

Soliton resolution
Let us first recall what is known about equation (1) posed on the whole space r ≥ 0. In this case, equation (1) is invariant under scaling, i.e. if W (t, r) is a solution, so is W λ (t, r) = W (t/λ, r/λ) for any positive number λ. The conserved energy , which is a distinctive feature of the critical dimension d = 4 allowing for the existence of nontrivial static solutions in the presence of scaling symmetry. These static solutions, hereafter called kinks (also referred to in the literature as instantons, solitons, or bubbles), form a one-parameter family with energy E 0 [Q λ ] = 4 3 which is the minimum energy for solutions interpolating between different vacuum states W = ±1 at the origin and at infinity. Obviously, the antikink −Q λ (r) is also the solution with the same energy.
Jendrej and Lawrie proved (see Theorem 1 in [4]) that any finite-energy solution of equation (1) posed on the whole space tends (modulo sign) either to the vacuum W = 1 or to an alternating chain of N rescaled kinks and antikinks 1

+
Here λ j (t) are continuous positive functions such that for each j = 1, ..., N where by convention λ N +1 (t) = t (in the global case) or λ N +1 (t) = T − t (in the blowup case), corresponding to the non-existing self-similar expansion and collapse.
We return now to our toy model and make the soliton resolution conjecture: 1 Strictly speaking [4] deals with 2d equivariant wave maps which split into equivariance classes indexed by a positive integer k. The case k = 2 is essentially equivalent to the 4d equivariant YM. The only qualitative difference is that for wave maps there are infinitely many topological sectors, while for YM there are only two sectors (modulo sign). For this reason, in the case of wave maps the chain of kinks and antikinks in Theorem 1 in [4] need not be alternating.

Conjecture 1.
Any smooth, finite-energy solution W (t, r) of equation (1) subject to the boundary condition W (t, 1) = 0 tends for t → ∞ (modulo sign) either to the half-kink or to the rescaled half-kink plus an alternating chain of N rescaled kinks and antikinks: Here λ j (t) are continuous positive functions such that for each j = 1, ..., N In the rest of the paper we confirm this conjecture for N = 0 and N = 1 and determine the rate of convergence to the attractors. In addition, we find a borderline between the basins of attraction using bisection along an interpolating one-parameter family of initial data.
Before presenting the results of numerical simulations, we wish to put forward a heuristic argument based on the method of collective coordinates that helps to understand some aspects of asymptotic dynamics. According to (7), in the case N = 1 the attractor has the following form: where the formula for µ(t) follows from the boundary condition W (t, 1) = 0. Inserting this ansatz into the lagrangian and integrating over r, in the limit of large λ we get the effective lagrangian (we retain only the first two leading terms) Thus, the λ-particle starting at some large λ(0) with velocityλ(0) > 0 escapes to infinity if E eff ≥ 2, while if E eff < 2 it reaches a turning point in a finite time. The ODE (12) provides a qualitative picture of the attractive interaction between the anti-half-kink and the expanding outer kink. This approximation ceases to work when the outer kink starts shrinking because the PDE solution is no longer close to the ansatz (9) (see Fig. 4 below). At the quantitative level the predictions of the effective model should be taken with caution because the ansatz (9) neglects radiation. In particular, according to Conjecture 1 the expansion rateλ(t) in (7) must go to zero as t → ∞, whereas the ODE yieldṡ

Characteristic formulation
We now reformulate our problem as the characteristic initial boundary value problem. To this end we define new coordinates where the function g(x) is assumed to be smooth and satisfying the finite energy condition g(0) = 1. For such data, the solution w(u, x) remains smooth for all future times u > 0. Moreover, the results by Chruściel and collaborators [7,8] imply 2 that the following asymptotic expansion holds near x = 0: Inserting this expansion into equation (13) and equating the coefficients of the same powers of x, we obtain an infinite system of ordinary differential equations for the coefficients c n (u). This system can be solved recursively oneby-one starting from the radiation coefficient c 1 (u) which is free. For n = 2 we getċ 2 (u) = 0, hence the coefficient c 2 is constant (so called Newman-Penrose constant). For large n the nonzero coefficients c n (u) grow polynomially for u → ∞ which is a reflection of the well-known fact that the expansion (15) is not uniform; see e.g. [9].
Multiplying equation (13) by x −3 w u we get the local conservation law Integrating this over x and using (15), we obtain the energy loss formula where is the Bondi-type energy (hereafter just called energy). Note that E[w] is equal to the potential part of the total conserved energy E[W ]. In terms of x the half-kink reads It is the global minimizer of E[w] as follows from the Bogomolnyi inequality which is saturated on is non-increasing and bounded below, there exists a limit According to Conjecture 1 the endstate (which clearly must be non-radiative) has the form (7), hence the final energy is a sum of the energies of the half-kink and N kinks/antikinks In section 5 we describe the relaxation to the N = 0 and N = 1 attractors. In our numerical simulations we have not observed N ≥ 2 attractors which suggests that, if they exist, they are nongeneric.

Linearized dynamics near the half-kink
Let w = q(x) + xf (u, x). Substituting this into equation (13) we obtain where Dropping the nonlinear terms and the O(x 5 ) term in the potential, we get the linear equation (corresponding to linearization around w = 1 rather than q) This equation has an explicit solution General solutions of equation (24) for initial data with the vanishing NP constant behave similarly to f 0 (u, x) for u → ∞, i.e. they decay as u −5 in the interior (x > 0) and u −5/2 at future null infinity (x = 0). This can be shown directly, or by defining F (t, r) = x 5 f (u, x) and rewriting equation (24) in terms of the original coordinates (t, r). Then, F (t, r) satisfies the free radial wave equation in 6 + 1 dimensions F tt − F rr − 5 r F r = 0, for which the late-time pointwise decay F (t, r) ∼ (t − r) −5/2 (t + r) −5/2 is well known [10]. The rate of decay of linear perturbations about q(x) is the same because the term O(x 5 ) in the potential (23) is asymptotically negligible. The numerical verification of this claim is depicted in Fig. 1 where we also show that solutions with nonzero NP constant exhibit slower decay. Having determined the late-time linear tail, now we turn to the computation of quasinormal modes. They govern the relaxation to the half-kink for intermediate times before the tail is uncovered. It is convenient to rewrite the linear part of equation (22) in terms of y = x 2 . Substituting we get the eigenvalue problem with v(1) = 0. Following Leaver [11] we seek solutions of (27) in terms of the power series v(y) = n≥1 a n (1 − y) n , a 1 = 1.
Since the nearest singularity from y = 1 is located at y = 0, this power series is absolutely convergent for y ∈ (0, 1]. The eigenvalues s n , called quasinormal frequencies, are selected by the condition that the power series is absolutely convergent at y = 0; as follows from (26), the corresponding solutions f n (u, y), called the quasinormal modes, are purely outgoing 3 . Inserting (28) into equation (27) we get a 7-term recurrence relation. Among its six linearly independent solutions, four solutions decay as a n ∼ 2 − n 2 for n → ∞, hence they do not affect the convergence properties of the series at y = 0. Using the method of successive approximations [12] one can show that the remaining two solutions have the following asymptotic expansions where the coefficients c (±) k can be determined successively by plugging the expansions (29) into the recurrence relation. We conclude that for n → ∞ a n = C + (s)a (+) For | arg(s)| < π, the series |a (+) n | diverges while the series |a (−) n | converges, therefore the quasinormal frequencies are given by the roots of the coefficient C + (s). There are several alternative ways to find these roots. The most frequently used is Leaver's method of continued fractions [11]. It is stable and accurate but tedious in the case at hand because the recurrence relation must first be reduced to three terms by Gaussian elimination. Employing this method we found exactly one pair of complex conjugate frequencies s ≈ −0.364322±0.476858i. To verify this result, we have reproduced it by two different methods: a brute force evaluation of a dominant solution by forward recurrence and an algebraic method introduced in [17]. We skip the details of these straightforward but dull computations. We confirmed the above perturbative analysis by the direct numerical integration of the linearized equation; see Fig. 2. Remark. It is instructive to compare the above computation of the quasinormal modes for the half-kink with an analogous computation for the vacuum solution w = 1 of equation (13) with the boundary condition w(u, 1) = 1 (as mentioned above, this is equivalent to the free wave equation in 6 + 1 dimensions). For the ansatz w = 1 + √ y e su v(y) (where y = x 2 ), we obtain the same eigenvalue equation as (27) but with the potential V = 15/4. Repeating the above analysis, we get a three-term recurrence relation having two linearly independent solutions a (±) n with the asymptotic behavior (29), hence as before the quasinormal frequencies are given by the roots of the coefficient C + (s) of the dominant solution. We remark that in this case the analysis based on the recurrence relation is purely academic because the eigenvalue equation can be solved exactly and the outgoing solution is given by v out = y − 1 2 e s/y K 2 (s/y), where K 2 (z) is the modified Bessel function of the second kind. Thus, the quantization condition for the quasinormal modes is K 2 (s) = 0, which has exactly one pair of complex conjugate zeros on the principal branch s = −1.281373 ± 0.4294849i [18]. We verified that the roots of C + (s) are the same, which provides a reassuring benchmark test for Leaver's method.

Numerical results
In this section we corroborate Conjecture 1 with direct numerical simulations of the initial boundary value problem (13)- (14). As in section 4, we write w = q(x) + xf (u, x) and then solve equation (22) numerically using the method of lines. To this end, we first discretize equation (22) in space using the pseudospectral approach. For numerical convenience, we rescale the spatial domain to the interval [−1, 1] using z = 2x − 1 and work with function values {f j (u) ≡ f (u, z j )} evaluated at K Chebyshev points of the second kind {z j = cos (j−1)π K−1 }, 1 ≤ j ≤ K. Spatial derivatives ∂ z and ∂ 2 z are replaced by the corresponding spectral differentiation matrices D (1) K and D (2) K [21] and then both the derivatives and nonlinear terms are evaluated using the grid function {f i }. The resulting semi-discrete system takes the following schematic form where the boundary condition f 1 = 0 replaces the equation at the grid point z 1 = 1. We bring this system to an explicit form by solving the linear system for u derivatives of {f j }. This requires inverting the operator D N with the first row replaced by a condition ∂ u f 1 = 0. Note that this linear operator is invertible and non-degenerate. The solution uses the LU decomposition of the resulting matrix. For efficient time integration we use an implicit scheme. We employ the BDF method (variable-order, variable-coefficient, in fixed-leadingcoefficient form) which, for the best performance, we limit to the second order (higher-order methods struggled to find the optimal step size/order, probably due to the stiffness of the equation). We used the IDA code [22], as available in Wolfram Mathematica [23], in which we set the error tolerances to very conservative values (typically 10 −13 ) so that the spatial resolution determines the errors in the numerical solution. Tests of the final algorithm show the spectral (exponential) convergence with increasing K.
Using the above method we have simulated the evolution of various initial data. Here we illustrate the results for a sample one-parameter family where b is a free parameter. For this data the NP constant is equal to zero. To see how a nonzero NP constant affects the dynamics we look in parallel at the evolution of initial data (32) with the additional term x(1 − x). The energy of initial data (32) attains the minimum value 0.6672 at b ≈ −2.1022. For −7.7295 b 2.5933 we have E < 2, hence according to (21) the halfkink q is the only possible attractor (the N = 0 case in our terminology).
In agreement with this, we observe rapid convergence to q through a short ringdown and then the late-time tail; see Fig. 3. Notice that the linear decay rates determined in the previous section, namely u −5 in the interior (x > 0) and u −5/2 along the future null infinity (x = 0), are not propagated by the nonlinear flow 4 for which the decay rates are slower by one power of u. Most important for us is the decay rate for the radiation coefficient c 1 (u) ∼ u −3/2 , however establishing this fact rigorously is a task that goes beyond the scope of this paper. We remark that similar nonlinear tails (but only in the interior) have been studied in the literature for semilinear wave equations in high even spatial dimensions (in particular, for the quadratic wave equation in 6 + 1 dimensions which is relevant in our context); see [19] and [20]. Next, we consider initial data with energy greater than 2. For moderate positive values of b the solution again quickly converges to the half-kink, however for larger values of b we observe formation of the superposition of the anti-half-kink and the expanding kink (the N = 1 attractor in our terminology). We find that the transition between these two scenarios occurs at b 0 ≈ 12.458288341217909. For marginally subcritical solutions (i.e. for b = b 0 − ε with small positive ε) a superposition of the anti-half-kink and the expanding kink appears for intermediate times but at a later time the expansion stops, the kink starts shrinking and is quickly annihilated. In this process the energy of the kink is rapidly radiated away and the solution settles down to the half-kink. This behavior is shown in Figs. 4 and 5 where we plot the snapshots of marginally subcritical and supercritical solutions and the corresponding energies, respectively. For a more quantitative description of the expanding phase of subcritical solutions, let x 0 (u) be the zero of the solution w(u, x) and u R be the return time when the expansion stops. We find that u R ∼ ε −1 and x 0 (u R ) ∼ ε 1/4 ; see Fig. 6. Translating these scaling relations to the original variables and comparing with the ansatz (9) we get t R ∼ ε −1 and λ(t R ) ∼ ε −1/2 . This is in agreement with the ODE approximation (12) which predicts that for marginally subthreshold effective energies we have λ(t R ) ∼ (2 − E eff ) −1/2 . Increasing b we find that above b 1 ≈ 47.90418049175238 the solution again settles down on the half-kink after a few rapid nonlinear oscillations. A similar transition occurs below b −1 ≈ −53.94479194728988. As |b| grows further the endstate keeps flipping back and forth between the N = 0 and N = 1 attractors. We conjecture that there are infinitely many critical values b n (n ∈ Z) at which the curve of initial data (32) intersects the N = 0 and N = 1 basins of attraction.
In the case N = 1 of Conjecture 1, it remains to verify that the speed of expansion of the outer kink goes asymptotically to zero, i.e. λ(t) t → 0 for t → ∞. This is shown in Fig. 7. Unfortunately, we are not in position to say more about the dynamics of λ(t). The computation of the precise asymptotic behavior of λ(t), which takes into account the loss of energy by radiation, is a challenging open problem that we leave to future work 5 .