Fermions on wobbling kinks: normal versus quasinormal modes

The system consisting of a fermion in the background of a wobbling kink is studied in this paper. To investigate the impact of the wobbling on the fermion-kink interaction, we employ the time-dependent perturbation theory formalism in quantum mechanics. To do so, we compute the transition probabilities between states given in terms of the Bogoliubov coefficients. We derive Fermi's golden rule for the model, which allows the transition to the continuum at a constant rate if the fermion-kink coupling constant is smaller than the wobbling frequency. Moreover, we study the system replacing the shape mode with a quasinormal mode. In this case, the transition rate to continuum decays in time due to the leakage of the mode, and the final transition probability decreases sharply for large coupling constants in a way that is analogous to Fermi's golden rule. Throughout the paper, we compare the perturbative results with numerical simulations and show that they are in good agreement.


Introduction
Solitons and solitary waves are important solutions of field theories that are topologically stable, propagate without losing the shape, and have localized energy density [1][2][3]. These solutions may couple to bosons and fermions with fascinating properties such as fractional fermionic charges [4,5]. In particular, kinks are solutions of scalar field theories when a discrete symmetry is broken and degenerate vacua appear. The prototypical models of kinks are the sine-Gordon, an integrable model, and the φ 4 model, which is non-integrable. The system composed of a fermion interacting with background φ 4 kink is analytically solvable and exhibits bound states as well as scattering states [6][7][8]. The properties of fermion-kink systems, such as the fermion energy spectrum and Casimir energy, have been studied for many different models [9][10][11][12][13][14]. In these systems, it is also possible to consider the back-reaction of the fermion on the soliton, which modifies the kink profile, especially for strong coupling constants [15][16][17][18]. This problem becomes even more fascinating when one considers the fermion spectrum in the presence of a static kink-antikink configuration. It was shown that the spectrum approaches two copies of the one of a single kink in the limit of large separation for both the φ 4 model [6] and the sine-Gordon [19], as expected. A more interesting problem is the time evolution of the fermion in the presence of a scalar field that is not static anymore and consists of a kink-antikink collision [20][21][22]. In this case, the fermions may be transferred from one soliton to the other [20,21] with decreasing probability as the system approaches the Bogomol'nyi-Prasad-Sommerfield limit, where the scalar field interaction vanishes [22]. An exciting application to this problem appears in cosmology, where domain walls can describe branes that live in higher-dimensional universes leading to fermions localization in the extra dimensions [23][24][25]. In this picture, we could live in a domain wall [26], and wall collisions could cause the big-bang [27].
The kink has a shape mode in φ 4 and many other scalar field models and can wobble. One of the earliest works to report this property was [28], which is now well understood in terms of perturbation theory [29][30][31]. In one of the earliest works of kink-antikink collisions, the authors considered an effective model to describe kink-antikink interactions with two collective coordinates, the position and the vibrational amplitude [32]. After that, in a triplet of seminal papers [33][34][35], the authors realized that the vibrational mode existence was essential to explain the resonance phenomenon where the kink and antikink collide multiple times before separating since the translational mode exchanges energy with the vibrational mode at each bounce. The structure of the resonance was shown to have a fractal structure [36] and was found in several models [37][38][39][40][41][42][43][44]. Surprisingly, the φ 6 , where the isolated kinks do not have a vibrational mode, also exhibits resonance windows. In this case, it was shown that the energy is stored in the vibrational mode of the compound kink-antikink configuration [45]. Interestingly, the resonance exchange mechanism only received a satisfactorily quantitative confirmation very recently [46,47].
To understand the role of the shape mode at each bounce, some authors considered collisions between initially wobbling kinks and showed that the initial vibrational energy could be converted into translational energy [48,49]. It is clear that the kink's vibrational mode is essential and that, even if the kink is not vibrating initially, this mode will get excited after a kink-antikink collision. Therefore, if a fermion interacts with a kink, it is important to ask how this interaction is modified when the kink starts wobbling.
An interesting modification of the usual scalar field models occurs when the potential is constructed to turn the vibrational mode into a quasinormal mode. In fact, many models already possess a tower of QNM for the stability equations around a kink solution. The quasinormal modes of the Schrödinger and Klein-Gordon equations can be found analytically for a few potentials, as described in [50]. See also [51] for a similar analysis. However, we are interested in the shape mode solution of the stability equation, which is the most relevant in the kink context. This mode can be turned into a quasinormal mode if one modifies the scalar potential, which was shown to suppress the resonance windows, as the decay rate of the quasinormal mode is increased [52,53].
In this work, we consider a fermion in the presence of wobbling kinks with both normal and quasinormal modes. In our case, the wobbling kink solution is considered as a background, that is, without a back-reaction. This is a common simplifying approximation that allows analytical treatment and is almost exact, in the limit where the coupling between the fermion and the kink is small [6,7]. We will show that the transition probabilities between the fermion bound and continuum states can be computed by utilizing the time-dependent perturbation theory formalism in quantum mechanics. We find that the fermions radiate to infinity if the coupling constant is small. This effect is quantified by Fermi's Golden Rule and can be reduced if the vibrational mode is turned into a quasinormal mode. In the next section, we discuss the framework we will use to study our system. In section 3 we apply our framework to the φ 4 model and in section 4 to a toy model that exhibits quasinormal besides normal modes. Finally, in section 5, we summarize our conclusions.

Model
Consider the following model representing a scalar field interacting with a fermion field in We ignore the back-reaction of the fermion field on the scalar field, which simplifies the system considerably and makes it possible to perform some analytical analysis. This is only a good approximation to the full equations of motion when g 1. However, we will also extrapolate our results for larger values of g. Besides that, it creates the possibility to explore the effect of fermion excited states, which only appear when g ≥ 1.
The potential is chosen such that the scalar field has a kink solution φ k . In this case, the fermion equation of motion in the presence of the kink is given by Taking the representation γ 0 = σ 1 , γ 1 = iσ 3 for the gamma matrices, where σ i are Pauli matrices and multiplying the above equation by γ 0 , it becomes with the Hamiltonian operator in the form This operator can be diagonalized by solving the Schrödinger-like equations − ∂ 2 x ψ n,± + V ± ψ n,± = E 2 n ψ n,± .
with the effective potentials given by V ± = g(gφ 2 k ∓ ∂ x φ k ) and the eigenfunctions by ψ n = (ψ n,+ , ψ n,− ) T . Because they are partner Hamiltonians [54] it is guaranteed that both have the same eigenvalues, except for the zero mode (E = 0), and that the eigenfunctions obey the first order equation Now, if we perturb the kink with the shape mode, the scalar field becomes φ = φ k + Aη S cos(ω S t) and the Dirac equation gives where the perturbation is given by H 1 (t) = gη S cos(ω S t)σ 1 and A is the amplitude of the perturbation. Then, assuming that the perturbation term is small compared with H 0 , one can write the fermion field as a linear combination of H 0 eigenfunctions, which in the Schrödinger picture evolves as ψ = n c n (t)e −iEnt ψ n .
We expand c n (t) in powers of A in the following form Considering the initial condition c n (0) = δ ni , first-order time-dependent perturbation theory gives [55] c (1) where ω ni = E n − E i and the matrix element of any operator O is defined as usual For a given i, the coefficients c n (t) are the Bogoliubov coefficients, which quantify the probability amplitude of the transitions between states [21]. Thus, the transition probability for n = i is given by This result is valid as long as A 2 |c (1) We can perform the time integration in eq. (10) to obtain For large t, one can compute the transition rate from the initial to the final state. It is defined as the large t limit of the ratio between the transition probability P i→n and time t. According to Fermi's golden rule, it is given by for n = i, which must be integrated over a set of final states.

φ 4 model
Let us start with the perturbative analysis on the fermion field interacting with the kink of the famous φ 4 model where the system is analytically solvable [6,7]. We need to set In this case, the effective fermionic potentials are in the Pöschl-Teller form They have analytical solutions, which are listed in Appendix A. In the following analysis, we will focus on the case where the initial fermion state is the zero mode with i = 0. Using the analytical expressions for the bound and scattering states, we can compute the matrix elements n|η S σ 1 |0 .
For the φ 4 model, a fermion bound state appears whenever g is an integer. The first excited state appears when g = 1.0, and the matrix element between this state and the zero mode is given by for example. Because ψ 0,+ and ψ 1,− are even and η S is odd, the first term vanishes and, because ψ 0,− = 0, the second term also vanishes. In general, the integral vanishes by parity arguments whenever the difference between the index of the two bound states is an odd Figure 1: Matrix elements for the operator η S σ 1 between the eigenstates of H 0 with fermion interacting with the φ 4 wobbling kink. We fix g = 1.6 in (b). Figure 2: Transition probabilities from the zero mode to continuum and bound states for a perturbation with amplitude A = 0.1. We compare simulation results with first-order perturbation theory. The rate from the zero mode to the continuum (eq. (17)) is also shown in (a) with the dotted line.
number. In contrast, the transition from the zero mode to the next excited state, appearing at g = 2.0, is allowed. In Fig. 1, some of the matrix elements n|η S σ 1 |i as a function of the coupling between the fermion and wobbling kink, g, are shown in the left panel. In the right panel, we depict the transition from the zero mode and also the first excited state to any continuum state 1 |k , varying k. This shows the amount of fermion radiating away for the fixed value of the coupling g = 1.6 and, in general, these transitions are allowed. We can also see that the radiation to the high continuum state, i.e., large k, is suppressed as the energy gap between the initial and final states becomes too large, making the transition impossible. We observe the same behavior in the other models studied here. The time evolution of the transition probability is shown in Fig. 2. In the figure, we compare a numerical simulation of the full equations of motion (eq. (3)) with the approximate result of first-order perturbation theory (eq. (13)). The numerical methods used in the simulations are described in Appendix B. In these calculations we take the analytical expression for the scalar field φ = φ k + Aη S cos(ω S t). We set A = 0.1 because it is of the same order of magnitude of typical wobbling amplitudes of a kink after a kink-antikink collision [48]. As one can see, the result of the simulation and the perturbative one match well. For g = 0.8, we compute the transition probability from the zero mode to the continuum. This is found semi-analytically by integrating eq. (13) over all the continuum states. We see that the probability increases at a constant rate, except for a small oscillation with frequency ω S , the kink wobbling frequency. This increase is expected for large times whenever Fermi's golden rule allows the transition to the continuum. For a final state in the continuum and the initial state in the zero mode, we have that ω ki = ± g 2 + k 2 . Thus, the transition to some state in the continuum is allowed by Fermi's golden rule whenever g < ω S = √ 3. More explicitly, we integrate eq. (14) over k from 0 to ∞ for the scattering states given in Appendix A, ψ k,L , with incident waves from the left. Then, including a factor of 4 to account for the scattering states with incident waves from the right as well as the ones with negative energy, it yields for g < ω S and zero otherwise. The transition rate is shown in Fig. 2a, where one can see that it gives the correct slope of the transition probability curve. This is an interesting result because it allows one to probe whether a kink is wobbling or at least to set an upper bound on the wobbling amplitude by measuring the fermions escape rate from the kink. For g = 2.1, we see that the transition to the continuum is not allowed anymore by Fermi's golden rule. Therefore, the transition probability to the continuum eventually stops increasing returns to zero, and this is repeated many times. Furthermore, the fermion now has a second excited state, and the transition probability to this state is finite and does not keep increasing for the same reason. In Fig. 2b we plot the sum of the transition probability for both second excited states with positive and negative energy. By symmetry, they both have equal transition probabilities and, therefore, the sum is just twice the transition probability to the second excited state with positive energy.

Toy model
Now we repeat the same procedure for a toy model, which allows kinks with either a normal mode or a quasinormal mode [53]. After rescaling, the toy model is described by the following potentials and The latter potential allows for the appearance of quasinormal modes, in contrast with the former one (18). All quantities in the definition of the potentials are constants, except for φ. The dependence on φ for negative values can be found by assuming that the potentials are even functions. The potential (18) is constructed such that the linearized potential of the stability equation has a square-well shape. Eq. (19) modifies this potential to the square-well shape with two barriers at the sides, allowing the bound states of the previous potential to tunnel and escape from the kink if γ is small enough. To ensure that the kinks behave correctly, we imposed continuity of the potentials and their derivatives. In this case, it is easy to show that, in the limit that goes to zero, eq. (19) becomes eq. (18). The parameter encodes the information about the energy leak when normal mode turns into the quasinormal mode, and as it increases, the energy leakage increases. Another consequence of the continuity relations is that we have only one free parameter in eq. (18), which we fix by setting γ = 3. This value is interesting because the kink has only a single vibrational mode and, thus, behaves similarly to the φ 4 model. For eq. (19), the continuity relations imply that there are three free parameters, which we fix at γ = 3.0, γ = 1.0 and = 0.05. Choosing these values, we guarantee that the only shape mode in the square-well stability potential turns into a quasinormal mode. The kink solution for the normal mode case is given by for positive x, where L is the point where φ 1 = A sin( √ λL). For the quasinormal mode case the kink is given by with the following definitions and Similarly, L and L are defined as the points where φ K = φ 1 and φ K = 1 − , respectively. The fermion spectra for both kinks, resulting from eq. (5), are shown in Fig. 3. They are very similar, except that the bound states appear earlier in the quasinormal mode case. This occurs because, in this case, the kink's tail decays more slowly, and the potentials V ± become wider. A typical evolution of the fermion field in the presence of a wobbling kink is shown in Fig. 4. These were obtained by solving Dirac's equation numerically and starting the fermion at the zero mode. The plots are similar for all the models that we considered. The kink stays bound to the kink at the origin but with noticeable oscillations. As g is increased, the fermion becomes more localized at the kink because the potentials V ± become deeper as g increases. However, the effect of the perturbations may increase with g in a certain range and have the opposite effect. If the transition to the continuum is allowed according to Fermi's golden rule, it is possible to see some radiation as shown in Fig. 4a and 4c. For the quasinormal case, the oscillations of the fermion field decrease slowly due to a decrease in the amplitude of quasinormal mode perturbation. This occurs because the quasinormal mode can tunnel the potential barrier created by the kink and radiate to infinity. However, the decay is slow and hence difficult to observe in Fig. 4c and 4d. We will discuss the perturbed system for both normal and quasinormal mode cases in detail and compare both numerical and perturbative solutions in the following subsections.

Normal mode
As we mentioned before, considering γ = 3.0, the system possesses one shape mode, which can be found analytically [53]. It has the following form with arbitrary normalization where we defined k S ≡ γ − ω 2 S and p S ≡ ω 2 S + 1 with ω S being the solution of the transcedental equation k S = −p S cot(p S L). For γ = 3.0 we find ω S 1.568. The matrix elements of the operator η S σ 1 between the eigenstates of H 0 are shown in Fig. 5. It is clear that the result is similar to the φ 4 case.
Starting with the fermion on the zero mode we can evolve the system integrating the equations of motion subject to the perturbation numerically. These results are shown in Fig. 6, where we compare them with the analytical expressions. Again we find results very similar to the φ 4 model, including the increase of the transition probability for g < ω S with transition rate given by the Fermi's Golden Rule expression in eq. (17). In fact, we expect similar results for a large class of models with kink solutions containing a shape mode.
The main difference between the φ 4 and the toy models is that H 0 has to be diagonalized numerically in the latter. We find the bound states using the NDEigensystem method in Mathematica and the solve_bvp method of the Python library SciPy. The continuum states are found as described in [12]. The method deals with the scattering states instead of the parity eigenstates. However, a simple linear combination of the scattering states gives the parity eigenstates. We consider, for instance, the scattering state with an incoming wave from the left with wavenumber k and energy E k = g 2 + k 2 . In short, the method consists of factoring out the highly oscillatory transmitted wave dependence as follows Then, we find the initial condition for χ k by taking the limit x → ∞ in the eigenvalue equation (6). Using φ k (∞) = 1.0 and χ k (∞) = 0, this results in (a) g = 0.8 Figure 6: Transition probabilities from the zero mode to continuum and bound states for the toy model. We consider the normal mode case with parameters A = 0.1, γ = 3.0. We compare simulation results with first order perturbation theory. The rate from the zero mode to the continuum (eq. (17)) is also shown in (a) with dotted line.
Numerically, we set infinity at x = 20.0, and we integrate χ k up to x = −20.0 according to the eigenvalue equation (6). After that, we multiply all eigenfunctions by appropriate normalization constants.

Quasinormal mode
For the toy model with potential (18), there is already a tower of quasinormal modes of the stability equation. This is also true for the φ 4 model with non-integer g. However, physically, the most important solution of the stability equation is the shape mode for the kink systems. We want this mode to be turned into a quasinormal mode, and this can be achieved by considering potential (19). In this case, it is possible to find this mode analytically, and it could be used as our perturbation to the scalar field if carefully truncated. This method can work well for large times. However, as we are interested in the whole transition process and not only the large-time behavior, a more sensible approach is to use the shape mode of the model with = 0.0 as the initial condition of the perturbation for the = 0 case. In this work, we will adopt this approach. The kink's shape mode becomes a quasinormal mode because the linearized potential of the normal mode case has been modified such that this mode can tunnel through the potential barrier. Therefore, this perturbation describes precisely the quasi-bound state localized at the kink, which decays in time according to the decay rate of the quasinormal mode. Writing φ = φ k +Aη(x, t), the evolution of this perturbation can then be found numerically using the time-dependent version of the stability equation where U (x) is the linearized potential. In this case the perturbed Hamiltonian becomes H 1 (t) = gη(x, t)σ 1 . If we neglect the radiation generated by the evolution, we can approximate the perturbation by η(x, t) η S cos(Ωt)e −Γt , where η S is given by eq. (24) and Ω and Γ are the analytical values of the frequency and decay rate of the quasinormal mode. For γ = 3.0, γ = 1.0 and = 0.05, we find Γ 0.0564 and Ω 1.584 according to the analytical expression in [53]. In this case, eq. (13) is modified to  If we wish to include the effect of the radiation, we should take eq. (10) to compute c (1) n with the full numerical solution of the scalar field perturbation as described above.
In Fig. 7 we show the time evolution of the transition probability. It behaves similarly to the previous cases, but the effect of the perturbation decays gradually in time due to the leak, which is expected from eq. (28). In the figure, we compare the numerical simulations with both perturbative results using eq. 28 and eq. (10) with the inclusion of radiation. The former does not have a good agreement for small values of g, where the effect of the radiation is more relevant but improves as g is increased.
Due to the decrease of the quasinormal mode's amplitude, the transition probabilities will reach an asymptotic value if we wait long enough. This can be used to make a more detailed analysis of the system's dependence on the coupling parameter g. The result is shown in Fig. 8, where we plot the asymptotic value of the transition probability as a function of g. The transition probability from the zero mode to the continuum, P 0→cont. , initially increases with g due to the dependence of H 1 on this parameter. After a critical value of g, the transition becomes forbidden because the energy gap between the zero mode and the continuum becomes too large, and consequently, the transition probability decreases considerably. This occurs for g Ω and is an effect analogous to Fermi's golden rule for normal modes. Near g = 1.54, the second excited state appears, and there is a finite transition probability to this state, which has a small peak and then also decreases considerably as g increases due to the large energy gap between the states.

Conclusion
In this paper, we considered a system composed of a fermion in the presence of a wobbling kink and treated the wobbling as a perturbation over the background kink configuration. We were able to write the Dirac equation for the fermion field interacting with the static kink background in the Hamiltonian formalism and diagonalize the unperturbed Hamiltonian. Moreover, we studied the fermion interaction with the wobbling φ 4 kink and also kinks of a toy model with excited normal and quasinormal modes, which have analytical expressions. The toy model potential possesses a parameter , which encodes the transition from the normal mode of the kink to a quasinormal mode as it becomes nonzero. As approaches zero, the toy model behaves similarly to the φ 4 model [53]. For the φ 4 model, we used the analytical results for the energy and eigenstates, while for the toy model, we diagonalized the Hamiltonian numerically. In both cases, the threshold energy increases linearly with g, and more bound states with positive and negative energy appear when g is increased. After diagonalizing, we include the perturbation and use the first-order perturbation theory to find the transition probabilities.
The perturbation theory results depend on the matrix elements of the Hamiltonian perturbation, which can be easily computed in all cases considered. The matrix elements from the bound states to the continuum are suppressed when the energy gap becomes too large. Moreover, we found that these matrix elements obey parity selection rules. The time dependence of the transition probabilities can be computed analytically, except for the quasinormal mode case, where the evolution is computed numerically. From the first-order perturbation theory results, we derived Fermi's golden rule for our system, which appears because the perturbation is periodic. It states that the probability rate to the continuum vanishes if the fermion-kink coupling g is greater than the wobbling frequency. The transition to the continuum was obtained by integrating the transition probability and transition rate over all the final momenta k for both positive and negative energy. For the transition rate, the integral was performed explicitly.
The perturbative results were compared with numerical simulations of the Dirac equation, which showed a satisfactory agreement. In the models containing a normal mode, the transition probability to the continuum shows a steady increase when g < ω S , in agreement with Fermi's golden rule. For larger g, the transition probabilities oscillate around an average value, leading to a vanishing transition rate over a large time interval. The simulations show the evolution of the fermion density with time. The fermion is mostly bound to the kink and, ignoring the effect of perturbation, it becomes more localized around the kink center as g is increased. On the other hand, the effect of the perturbation can be opposite to that in a certain range of g. Furthermore, for small coupling constants, it is possible to see more pronounced radiation in the evolution.
The scalar potential was modified to turn the normal mode of the kink into a quasinormal mode. In this case, the vibrational mode of the toy model with = 0 was considered as the initial condition of the perturbation. Then, we evolved the perturbation numerically according to the time-dependent stability equation. When the normal mode becomes a quasinormal mode, the perturbation amplitude decays over time because it leaks energy. This causes the transition rate to also decrease over time and, consequently, the probability reaches an asymptotic value for large times. Interestingly, the transition to the continuum is only significant for coupling constants smaller than the oscillation frequency, which is a reminiscent behavior of Fermi's golden rule for the normal mode case. These effects were observed in the numerical simulations and compared to results from first-order perturbation theory. We included the full evolution of the scalar field perturbation in the perturbation theory to have a good agreement between the two. However, we have also considered the approximation, neglecting the radiation in the scalar field because, in this case, it is possible to find the time dependence of the transition probabilities explicitly. The latter approximation gives a reasonable estimation of the asymptotic value of the transition probability.
Our work clarifies that the fermion decouples from the wobbling kink if the coupling constant is small and that this effect is attenuated if we have quasinormal instead of normal modes. This occurs because the perturbation decreases over time in the former case. Moreover, if it is possible to measure the fermion escape rate, one could, in principle, probe whether the kink is wobbling or set a bound on the wobbling parameters, such as the wobbling amplitude, quasinormal mode decay rate, and the coupling constant g. This could be relevant for electrons in polyacetylene [56,57] and in the scenario where we live in a domain wall [26], to name a few examples.

A Fermion Bound States of the φ 4 model
The bound and scattering states of a fermion coupled to a φ 4 kink are listed in references [6,7]. The effective fermionic potentials, in this case, are given by eq. (15), which are Pöschl-Teller potentials. In this case, the number of bound states is given by the largest integer smaller than g. The bound states are as follows ψ n,+ = N n + sech g−n (x) F −n, 2g − n + 1, g − n + 1; 1 2 (1 − tanh(x)) , ψ n,− = N n − sech g−n (x) F −n + 1, 2g − n, g − n + 1; where N n ± are normalization constants obeying N 0 − = 0 and, for n ≥ 1, N n − /N n + = n/E n . Also E n = n(2g − n) is the energy of the state and F is the hypergeometric function. On the other hand, the scattering states of an incoming incident wave from the left and from the right are respectively, where ζ ± = g ± 1 2 , and N k,L ± and N k,R ± are normalization factors obeying N k,L − /N k,L + = −N k,R − /N k,R + = (ik + g)/ k 2 + g 2 . Both scattering states have energy E k = g 2 + k 2 .

B Numerical technique
To solve the equations of motion numerically for the fermion field, we discretize space on the interval −L < x < L with spacing h = 0.1 and L = 100.0. We substitute the partial derivative with respect to x by a five-point stencil approximation which is fourth order accurate in h. At x = −L + h, the finite difference expression is as follows where we impose ψ(x = −L) = 0. Similar expressions are used at x = L. After discretizing the space, we integrate the resulting ordinary differential equation using a fifth-order Runge-Kutta method [58] implemented by the solve_ivp method of Python library SciPy. For the scalar field, we use the analytical expression for the kink and, whenever possible, for the perturbation. This way, we ignore the higher-order coupling of the shape mode with the radiation mode [29]. We cannot use the analytical expression for the quasinormal modes, and we must integrate the time-dependent stability eq. (27). In this case, we follow the same procedure to integrate the fermion field, i.e., the spacetime is discretized the same way, partial derivatives over x are obtained using a five-point-stencil approximation, and integration over time is done using a fifth-order Runge-Kutta method.