# Automation for Interpretable Machine Learning Through a Comparison of Loss
Functions to Regularisers

A. I. Parkes, J. Camilleri†, D. A. Hudson and A. J. Sobey∗
Maritime Engineering, University of Southampton, Southampton, SO17 1BJ
{A.I.Parkes, ajs502, J.Camilleri, dominic} @soton.ac.uk
∗Marine and Maritime Group, Data-Centric Engineering, The Alan Turing
Institute,
The British Library, London, NW1 2DB
†Silverstream Technologies Ltd, London, UK

###### Abstract

To increase the ubiquity of machine learning it needs to be automated.
Automation is cost-effective as it allows experts to spend less time tuning
the approach, which leads to shorter development times. However, while this
automation produces highly accurate architectures, they can be
uninterpretable, acting as ‘black-boxes’ which produce low conventional errors
but fail to model the underlying input-output relationships—the ground truth.
This paper explores the use of the Fit to Median Error measure in machine
learning regression automation, using evolutionary computation in order to
improve the approximation of the ground truth. When used alongside
conventional error measures it improves interpretability by regularising
learnt input-output relationships to the conditional median. It is compared to
traditional regularisers to illustrate that the use of the Fit to Median Error
produces regression neural networks which model more consistent input-output
relationships. The problem considered is ship power prediction using a fuel-
saving air lubrication system, which is highly stochastic in nature. The
networks optimised for their Fit to Median Error are shown to approximate the
ground truth more consistently, without sacrificing conventional Minkowski-r
error values.

##  1 Development of Interpretable Machine Learning

Machine learning regression models are increasingly being used in industrial
and engineering contexts for high-stakes decision making, automation and
control. These methods often produce low conventional error values, yet are
known to produce physically inconsistent results which cannot generalise off
test set and so cannot be relied upon to model the ground truth of the system.
The models are designed to be accurate, not interpretable, and so a human
cannot understand how changes in the inputs change the prediction. In real-
world applications, where the output can have a direct effect on human life or
the environment, model accuracy alone is not sufficient.

Trust is increased if a trained model approximates the true input-output
relationships, performing accurately within the bounds of the training data
set and beyond it. It has been demonstrated that for many applications
minimising traditional error measures cannot guarantee an accurate
approximation of the ground truth (Willard et al., 2020). This is due to a
poor inductive bias, the inherent prioritisation of one solution over another
(Battaglia et al., 2018), produced by conventional error measures which are
based on Minkowski-r metrics (Hanson and Burr, 1987).

This trust can be increased by manually tuning to remove overfitting or to
provide a solution that makes more sense to the user. However, the expert
knowledge and domain experience required to properly tune a machine learning
method manually are not always available in industry. Genetic algorithms are
therefore increasingly used to search a method’s hyperparameter space more
efficiently (Yang et al., 2021) (Kumar et al., 2021); which minimise
conventional error measures on a test set, often combined with lowering the
complexity of the network. This automation exacerbates the lack of
interpretability, as models have a large flexibility, and prediction accuracy
is prioritised, a low conventional error is achieved without certainty that
the method has modelled the correct internal functions. Regularisation
hyperparameters can be optimised alongside other neural network parameters
(Tani et al., 2021) (Luketina et al., 2016), which increases the search space
and creates more flexibility for methods to produce ‘accurate’ predictions and
avoid overfitting.

Common regression regularisation methods are l1 and l2 regularisation and
dropout. For l1 and l2 regularisation, large network weights are penalised in
the loss function (Nowlan and Hinton, 1992). The absolute value of the weights
is penalised in l1 regularisation and the squared value in l2, meaning l1
encourages weights towards zero and l2 encourages weights to be small but non-
zero. The l1, l2 and elastic net (l1+l2) regularisers improve a networks
generality, increasing the applications where the trained methods can be
applied, by penalising complexity. Dropout, where a randomly selected subset
of weights are optimised at each epoch rather than the full set, improve the
generality of the trained models by preventing co-adaption of weight values
(Srivastava et al., 2014). Dropout has been shown to be equivalent to l2
regularisation after scaling by Fisher information (Wager et al., 2013),
suggesting that the two should not be used in unison. The neural network
regularisation methods discussed above aim to improve generality, reducing
overfitting by simplifying the relationships modelled by the networks.

Regularisers improve the modelling of the ground truth in scenarios adhering
to the assumptions in the proof in Bishop (1995), under which minimum
Minkowski-r error values approximate the conditional average of the dataset.
This is because the inductive bias from the loss function guides the input-
output relationships towards the conditional average, while the regularisation
stops overfitting by simplfying the input-output relationships being modelled.
However, these assumptions are restrictive and it is noted that few regression
applications adhere to them. For example, one assumption is that the dataset
is homoscedastic. In scenarios not adhering to these assumptions, network
regularisation simplifies the relationships being modelled but this does not
necessarily improve the generality, or model the ground truth.

The Fit to Median Error measure (Parkes et al., 2021) produces more
interpretable regression, when used in conjunction with conventional error
measures. This is achieved by regularising the learnt input-output
relationships to the conditional median of the training dataset: the median
output value, conditioned on each isolated input variable in turn (Bishop,
1995). For many regression applications the conditional medians are a good
approximation of the ground truth input-output relationships but as yet it has
not been explored as part of an automated approach.

A challenging regression problem is ship power prediction for a vessel using
air lubrication to reduce fuel consumption. It is chosen to be used in this
study as it violates the assumptions in Bishop (1995), where the noise in the
output space is non-Gaussian and heteroscedastic. In this situation, correctly
modelling the ground truth and accurate prediction is required but there is
limited understanding of that ground truth (Parkes et al., 2018). The
literature shows that shaft powering of a vessel can be predicted with average
accuracies of between 1.5-5% error with the use of a regression neural network
trained with high frequency data from the vessel (Pedersen and Larsen, 2009),
(Petersen et al., 2012), (Le et al., 2020), (Jeon et al., 2018), (Liang et
al., 2019). All neural network applications to ship power prediction in the
literature use a combination of local searches and domain knowledge to
identify hyperparameter values. The addition of an air lubrication device
increases the complexity of the regression problem, as the system interacts
with a number of interrelated input variables.

This paper explores the automation of neural network training to a new
problem, with a focus on producing a network which accurately models the
ground truth. It compares the ground truth representation of a neural network
when a genetic algorithm optimises the network’s hyperparameters to reduce the
Mean Fit to Median Error measure and compares it to standard regularization
using l1, l2 and dropout, and to a network optimised to minimise the Maximum
Absolute Error. It is illustrated that neural network regularisation methods
(l1, l2 and dropout) can be replaced by the use of the Mean Fit to Median
performance measure as an objective in the genetic algorithm, reducing the
complexity of the search space and producing networks which more consistently
model the ground truth.

##  2 Neural Networks Parameters

Previous applications of neural networks to ship power prediction use between
1 and 3 hidden layers (Leifsson et al., 2008) (Parkes et al., 2019), and
between 5 and 300 neurons in each hidden layer (Jeon et al., 2018). To provide
a sufficiently large search space to allow verification, or otherwise, of
these parameters a maximum of 4 hidden layers and 1000 neurons in each layer
are used. The majority of the literature treats the problem as time-invariant
and use feed-forward networks, so no recurrent parameters are optimised. As
the optimiser or activation functions are rarely documented in the literature,
the state-of-the-art optimisers and activation functions available in the
Keras framework (Chollet et al., 2015) are used in the optimisation, Table 1.

Table 1: Selected Neural Network Hyperparameters Hyperparameter |  Value or set
---|---
Layers |  [1,4]
Neurons in each layer |  [1,1000]
Epochs |  Increasing from 1-20 for increasing generations
Early stopping patience |  5
Loss function |  Mean Absolute Error
Performance measures |  Mean Absolute Relative Error, Maximum Absolute Relative Error, Mean Fit to Median Error
Optimiser |  SGD, Adam (Kingma and Ba, 2014), Nadam (Dozat, 2016),
|  RMSprop (Hinton et al., 2012), Adagrad (Duchi et al., 2011),
|  Adadelta (Zeiler, 2012), Adamax (Kingma and Ba, 2014)
Activation function |  ReLU, sigmoid, softmax, softplus, softsign, tanh, selu, elu
l1 & l2 Rates |  0, 0.01,0.001,0.0001,0.00001
Dropout |  [0,0.9)
Initialiser |  Random Normal ($\mu=0,\sigma=0.1$)

The number of epochs and early stopping procedure are not optimised, as there
was a need for predictable compute requirements and allowing the optimisation
of these parameters leads to unpredictable run times. The number of epochs to
train each network increases for increasing generation number in the genetic
algorithm, from 1 epoch in the first 15 generations to 20 in the final 15.
This was also implemented to reduce compute and it was validated that when
more than 20 epochs were allowed, that the early stopping, with a patience of
5, stopped the training within 20 epochs for the majority of networks. The
loss function is similarly not optimised, the Mean Absolute Error is used, as
the conditional medians are closer to the ground truth input-output
relationships in these datasets than the conditional means.

The performance measures, or the genetic algorithm’s fitness functions, are
the Mean Absolute Relative Error, the Maximum Absolute Relative Error and the
Mean Fit to Median Error. Different combinations of these, alongside the use
of regularisation parameters in the search space are compared to illustrate
the effect of different types of regularisation.

##  3 cMLSGA Parameters

Table 2: Selected cMLSGA Hyperparameters Hyperparameter |  Value or set
---|---
Algorithm at Individual Level |  HEIA, IBEA
Crossover Type & Rate |  SBX & DE, 1
Mutation Type & Rate |  Polynomial, 0.08
Number of eliminated collectives |  1
Generations between elimination |  10
Population size |  1000
Generations |  300
Proportion elite |  10%

In this study cMLSGA111The code for cMLSGA is available at
https://www.bitbucket.org/Pag1c18. is selected as it shows the top performance
on a range of evolutionary benchmarking problems (Grudniewski and Sobey, 2021)
and practical problems (Grudniewski and Sobey, 2019). Genetic algorithms are
increasing used to tune neural network hyperparameters including
regularisation parameters for use on new problems (Jin et al., 2004). Many
approaches have multiple genetic algorithm objectives, although these all
minimise an error measure and a measure of network complexity (Wang et al.,
2019) and (Smith and Jin, 2014). The use of multiple different performance
measures as objectives is yet to be explored in the literature.

Table 3: Genetic Algorithm Approaches Approach | Objective(s) | Network Regularisation
---|---|---
GAi | Mean Absolute Error | l1, l2 and dropout
GAii | Mean Absolute Error | l1, l2 and dropout
| Maximum Absolute Error |
GAiii | Mean Fit to Median Error | None
| Mean Absolute Error |
GAiv | Mean Absolute Error | None
| Maximum Absolute Error |

Four approaches are investigated in this study, summarised in Table 3, for
approach (GAi) and (GAii) the genetic algorithm cMLSGA optimises all variables
in Table 2, including the l1 and l2 regularisation rate and the dropout rate
of the networks. Although it is advised that l2 regularisation and dropout are
not used in the same network the genetic algorithms are provided with zero
options for all regularisation parameters, to identify if one is preferable in
this scenario.

Approach (GAi) is a single objective genetic algorithm optimising the Mean
Absolute Error which is compared to a multi-objective formulation where the
(GAii) approach optimises both Mean Absolute Error and Maximum Absolute Error.
For approaches (GAiii) and (GAiv) no network regularisation parameters are
optimised: l1, l2 and dropout rates are all set permanently to zero. They
avoid producing networks that have overfitted by the use of two performance
metrics as multi-objectives, (GAiii) uses the Mean Fit to Median and Mean
Absolute Errors to be minimised and (GAiv) uses the Maximum Absolute and Mean
Absolute. All approaches use 40 CPUs with 2.0 GHz Intel Skylake processors and
192 GB of DDR4 memory, and take less than 3 days, this setup may not be
feasible for widespread industrial application, although it is suggested it is
within reach of some industries.

##  4 Data

Figure 1: The distribution of the observed shaft powers for half knot bins of
speed through the water for dataset where the system is off. In the box and
whisker plots the boxes contain 50% of the distribution and the whiskers
extend to the datum which is at 1.5 times the interquartile range.

The data used in this study are from a large vessel equipped with the
Silverstream® Air Lubrication System. The air lubrication system works through
use of fluid sheering to create an air microbubble carpet directly captured
within the boundary layer on the ship hull bottom. The bubble carpet reduces
the frictional resistance thereby increasing the speed and reducing the shaft
power. Compressors provide a constant supply of air to the hull bottom to
maintain a uniform bubble carpet operated at the optimal compressor power that
maximises the energy balance. The study is performed on both system on and
system off datasets, however for brevity only results for system off are
presented as they show similar performance.This prediction is required for a
baseline determination of how the system is working, but the relationships
between the power, weather, ocean and operating conditions are complex and
difficult to model.

The variables considered in this study are the shaft power, speed through
water, relative wind speed and direction, draught and trim, with shaft power
the target variable. These are selected based on a detailed study into
variable selection for shaft power prediction (Parkes et al., 2019). The speed
through water is selected over the speed over ground, for use as an input
variable, as it is more hydrodymanically relevant and its accuracy is
validated by comparison to the speed over ground. The dataset is cleaned by
removing rows with missing or non-physical values and all datapoints below
0.05 normalised shaft power are removed. The dataset is split into two using
the air lubrication system status: system on and system off, where system on
is defined as air lubrication system power greater than zero. The system on
dataset contains 352,690 datapoints and system off contains 237,962. The data
is split into training, testing and validation sets of 70%, 15% and 15%
respectively. Each network in the genetic algorithm trains on a randomly
sampled 35,000 datapoints from the training set and uses randomly sampled sets
of size 7,500 from validation and testing sets for validation during training,
and testing to produce the fitness of the network for the genetic algorithm.
The errors stated in the paper are from networks on the Pareto fronts of each
approach, which are validated on the full testing set.

The datasets contain large regions of sparse data in all input variable
domains, this is exemplified by the ship speed domain where each half-knot
interval below 16 knots contains less than 0.8% of the data, which accounts
for more than half the speed domain, Figure 1. In addition, the boxplot ranges
and outliers show high heteroscedicity with idiosyncratic noise caused by
situations where the angle of the propeller blades is varied to achieve the
required speed. This highlights the complexity in developing models of the
powering of this vessel, as the dataset also contains the effects from other
latent variables, such as piloting behaviour and route taken.

##  5 Optimisation including regularisation parameters: (GAi) and (GAii)

Previous studies predicting ship powering using neural networks report that
l1, l2 and elastic net increase both test set and off-test set errors and that
optimal values for both l1 and l2 are zero. Therefore the genetic algorithm
setup is biased towards low and zero values of regularisation rates by using a
set of exponentially decreasing values and an explicit zero option.

(a) (GAi)

(b) (GAi)

(c) (GAii)

(d) (GAii)

Figure 2: Distribution of regularisation rates for networks in the last 15
generations of (GAi) cMLSGA with multi-objectives of minimising Maximum and
Mean Absolute Error for (a) l1 and (b) l2 and (GAii) cMLSGA with the single
objective of minimising Mean Absolute Error for (c) l1 and (d) l2

(a)

(b)

Figure 3: (a) Dropout rate for networks in the last 15 generations of cMLSGA
with (GAi) the single objective of minimising Mean Absolute Error and (GAii)
multi-objectives of minimising Maximum and Mean Absolute Error and (b) the
number of neurons in each layer for networks in the last 15 generations of
cMLSGA with (GAi), (GAii),(GAiii) and (GAiv).

The single objective (GAi) fails to identify that zero regularisation rates
produce the lowest errors, favouring networks with the highest possible rate
of l2 (0.01), Figure 2(b). (GAi) produces networks with the highest Mean
Absolute Relative Errors of all the approaches, $(5.19\pm 0.00)$% from Figure
4(a). In contrast, (GAii) favours lower l1 and l2 rates of 0 or 0.00001,
Figures 2(c) and 2(d), which results in networks with the lowest Mean Absolute
Relative Errors of all four approaches, on average, with a value of $(2.87\pm
0.45)\%$, shown in Figure 4(a). This is around 0.5% higher than the lowest
documented error for ship power prediction.

It is posited that the high error for the single objective problem is directly
related to the use of a large l2 regularisation rate, as noted in previous
studies for ship power prediction. It is possible that the use of a multi-
objective search algorithm for a single-objective problem means that the
optimal hyperparameters can’t be found, resulting in large errors. The
implementation also requires restrictions in the number of epochs used for
training in the initial generations, it is possible this biases (GAi) towards
certain size networks, where higher l2 rates are preferable. This hypothesis
is supported by the fact that 74.4% of networks in the first 15 generations of
(GAi) have 1 hidden layer, and that over 99.8% of the networks in the final 15
generations have 1 hidden layer, with $880\pm 17$ neurons in this layer,
Figure 3(b). This is significantly more neurons than those in the hidden layer
of networks in the final 15 generations of (GAii) which range from 3-952 with
a median value of 709, Figure 3(b). The added objective of minimising Maximum
Absolute Error in (GAii) may cause these slightly smaller networks to be more
attractive as they are in a sense regularised by their size, as they have
reduced modelling flexibility therefore are less likely to overfit and produce
high Maximum Absolute Errors.

Another explanation for the difference in l2 rates chosen by (GAi) and (GAii)
is the equivalence of l2 and dropout. Since l2 and dropout are equivalent up
to a Fisher transformation, their use in conjunction is not recommended. The
evidence for this is that (GAi) favours the highest l2 rate and has a median
dropout rate in the final 15 generations of 0.116, whereas (GAii) favours the
zero l2 rate and has a median dropout rate of 0.624, Figure 3(a). This
illustrates that the genetic algorithms will chose either l2 or dropout to
minimise the Mean Absolute Relative Error. The l1 rates also support this
hypothesis, as chosen rates for l1 regularisation in the final 15 generations
are more comparable for (GAi) and (GAii).

##  6 Optimisation using multiple performance measures: (GAiii) and (GAiv)

For approaches (GAiii) and (GAiv) all neural network regularisation parameters
are set to zero. The regularisation is performed by minimising different
network performance measures, the Mean Absolute and Mean Fit to Median for
(GAiii), and the Mean Absolute and Maximum Absolute for (GAiv). The trade-off
between the two objectives produces regularised neural networks, without
explicitly changing the architecture or loss function. The Mean Fit to Median
is chosen as it indicates how close the relationships modelled by a network
are to the conditional averages of the dataset, in many regression examples
this is akin to the ground truth input-output relationships (Parkes et al.,
2021). The Maximum Absolute is chosen as for many industrial applications of
machine learning the maximum prediction error is more pertinent than the mean
error. The Mean Absolute Error is used instead of the Mean Squared Error in
both approaches, as the conditional medians are closer to the ground truth
input-output relationships in these datasets than the conditional means.

Differently shaped networks are favoured by (GAiii) and (GAiv), compared to
(GAi) and (GAii), focusing on networks with 3 hidden layers and on average
less than 400 neurons in each layer, Figure 3(b). These networks have 51 times
the number of connections than the networks chosen in (GAi) and (GAii). Apart
from (GAi), (GAiii) has the most consistently sized networks in the final 15
generations, with an interquartile range of 46 neurons, compared to (GAiv)
which have an interquartile range of 131 neurons. It is suggested that as the
Mean Fit to Median Error biases networks towards specific input-output
relationships, there is a smaller range of potential network architectures
which habitually model these relationships. Whereas networks which minimise
the Maximum Absolute Error are less restricted and can model a wider range of
input and output relationships.

(a)

(b)

Figure 4: Mean Relative Absolute Error (a) and Maximum Absolute Error (b) from
cMLSGA with (GAi) the single objective of minimising Mean Absolute Error and
(GAii) multi-objectives of minimising Maximum and Mean Absolute Error, both
optimising the parameters for l1, l2 regularisation and dropout in the
networks, and (GAiii) and (GAiv) which do not use network regularisation but
minimise Mean Fit to Median and Maximum Absolute Error respectively, alongside
Mean Absolute Error

The Mean Absolute Relative Errors from networks in the Pareto fronts are
$(2.97\pm 0.25)\%$ for (GAiii) and $(3.10\pm 0.28)\%$ for (GAiv). It is
expected that (GAiv) would produce higher Mean Absolute Relative Errors as
discussed above, minimising the Maximum Absolute Error should bias predictions
towards the midpoint of the conditional output distributions, whereas
minimising the Mean Absolute Error should bias predictions towards the median
of these distributions. As it is established that noise in the output
distribution is non-Gaussian, Figure 1, these values will not align so some
sacrifice in Mean Absolute Error is expected from (GAiv). Both (GAiii) and
(GAiv) produce comparable Maximum Absolute Errors, of $(49.5\pm 1.1)$% and
$(51.9\pm 4.5)$%. It is suggested that this is because, although the
conditional median output value and conditional midpoint output value do not
align for the majority of the input domain, they are sufficiently close to
produce comparable Maximum Absolute Errors.

Across all four approaches, the genetic algorithm producing networks with the
highest Mean Absolute Error is the approach which does not provide extra
weighting to sparse areas of data. The approaches minimising Maximum Absolute
Error are implicitly biased away from networks which predict the majority of
the testing datapoints correctly, but predict one datapoint poorly, favouring
networks which predict all testing datapoints to a moderate degree of error.
Approach (GAiii) more explicitly weights prediction in sparse areas of data by
favouring networks which model the conditional median of the dataset across
all input domains, irrespective of the quantity of data across each input
domain. The regression problem of ship power prediction is chosen in part
because of it’s irregular data distribution; more than 9% of the dataset lies
in less than a 0.5 knot interval of ship speed, Figure 1. This provides an
explanation for the high testing errors from (GAi), where only the Mean
Absolute Error is minimised, there is little incentive for the genetic
algorithm to produce networks which generalise across the full range of the
input domain well.

##  7 Comparison of the interpretability

To asses the interpretability of the networks selected by the four different
approaches the learnt relationship between an input, the ship speed, and the
output, shaft power, for the networks in the Pareto front of each approach are
visualised, Figure 5. These are extracted with the following procedure: set
all but one input variable to be constant at the mode; cycle the remaining
variable from its minimum to its maximum recorded values with 150 points
evenly spaced along the domain and run the new dataset through the trained
network.

(a) (GAi)

(b) (GAii)

(c) (GAiii)

(d) (GAiv)

Figure 5: The learnt speed-power curves from 5 networks on the Pareto fronts
of (GAii), (GAiii), (GAiv) and the 5 networks producing lowest Mean Absolute
Relative Error from (GAi).

The approach which produces the most consistent speed-power relationships is
(GAi), with an average variation of 1.8%222Average variation in just the
speed-power curves are discussed in this section, but it is verified that all
input-power curves follow the same trends with variation around 0.5% across
input variables., Figure 5(a). However, the relationship modelled by the 5
networks with the lowest Mean Absolute Relative Error in (GAi) all approximate
a piece-wise linear relationship which clearly underfits the dataset in Figure
1. The expected trend between ship speed through the water and shaft power is
a cubic polynomial, therefore as well as producing the highest Mean Absolute
Relative Errors, networks chosen by (GAi) model the ground truth input-output
relationships the worst out of the four approaches. Both (GAii) and (GAiv)
produce 5 fairly consistent speed-power curves, with average variations of
5.9% and 10% respectively, Figures 5(b) and 5(d). Both approaches approximate
smooth polynomial curves, although the degrees of the polynomials might
differ, as multiple curves intersect at various points along the speed axis.
The spread of learnt relationships is greater at the highest and lowest speeds
for (GAii), with a decrease in spread for speeds of around 15 knots, where
many of the curves intersect. The curves from (GAiv) show equal spread across
the speed domain.

The approach with both accurate and consistent learnt speed-power curves is
(GAiii), with limited intersections of curves and an average spread of 3.0%.
It is suggested that the reason using the Mean Fit to Median Error as an
objective in a multi-objective genetic algorithm produces more interpretable
results, or more consistent learnt relationships, is because instead of
encouraging the networks to model more simple relationships. It encourages the
networks to model the conditional median functions of the dataset, supported
by the increase in network connections. Whereas the other approaches leave
room for networks to fail to model the conditional averages, especially in
irregularly distributed and non-normally distributed datasets. The Mean
Absolute Error values from networks selected by (GAiii) are on average 0.1%
higher than those from (GAii), and the Maximum Absolute Error values are 1.6%
higher.

A limitation of the approach is that the Fit to Median Error measure will
likely perform best at improving interpretability on datasets which violate
the assumptions in Bishop (1995); the ship powering example is chosen to
illustrate this as it provides a clearly heteroscedastic dataset. For
applications where noise profiles are Gaussian, and there is no effect from
latent or interrelated input variables, the Fit to Median Error will not
improve interpretability, but will perform the same as conventional
Minkowski-r metrics, either Mean Squared or Mean Absolute Error depending on
the convexity of input-output relationships.

Interestingly, the approach producing the lowest Mean Absolute and Maximum
Absolute Errors does not model the ground truth the most accurately. This
creates a potential for negative societal impacts, as the standard performance
metrics for regression neural networks do not provide a full picture of
performance or expected behaviour. Interpretability of trained methods is
essential for safe application of machine learning in the real world,
especially when automated methods are used to replace experienced
professionals. (GAii) demonstrates the same accuracy of approach as those with
standard network regularisation, but with a better fit to the ground truth.
This approach bypasses the need to use, and therefore to optimise the
parameters of the regularisation methods. If evolutionary computation is
already being used to optimise network parameters, then compute is saved by
removing the network regularisation parameters l1, l2 and dropout. (GAiii)
completed 300 generations in 46hours whereas (GAii) required 12 hours more
computation to complete 300 generations.

##  8 Conclusion

Interpretable and accurate methods are required for widespread application of
machine learning to real-world regression problems. To automate the training
of neural networks so that the result is interpretable, three different
genetic algorithm approaches are compared: one to minimise the Maximum
Absolute Error of the networks, which includes standard regularization using
l1, l2 and dropout; and two which do not use any network regularisation, one
minimising the Mean Fit to Median and one to minimise the Maximum Absolute
Error. The results show that all three approaches give similar Mean Absolute
Errors from networks on their Pareto fronts, from 2.9% for the approach with
regularisation to 3.1% for the approach minimising Maximum Absolute Error.
However, the Mean Fit to the Median approach shows a considerably better
interpretability, or fit to the ground truth, with a spread in predicted
input-output curves of 3% compared to a spread of 6% for the approach using
regularisation and 10% when minimising the Maximum Absolute Error.

## Acknowledgments and Disclosure of Funding

The work was kindly supported by the Lloyds Register Foundation. The project
is also supported by the KTN with a Knowledge Transfer Partnership between the
University of Southampton and Silverstream Technologies, who the authors would
like to thank for providing the dataset used in this study and their support,
especially David Connolly. The authors acknowledge the use of the IRIDIS High
Performance Computing Facility, and associated support services at the
University of Southampton, in the completion of this work.

## References

  * (1)
  * Battaglia et al. (2018) Battaglia, P. W., Hamrick, J. B., Bapst, V., Sanchez-Gonzalez, A., Zambaldi, V., Malinowski, M., Tacchetti, A., Raposo, D., Santoro, A., Faulkner, R., Gulcehre, C., Song, F., Ballard, A., Gilmer, J., Dahl, G., Vaswani, A., Allen, K., Nash, C., Langston, V., Dyer, C., Heess, N., Wierstra, D., Kohli, P., Botvinick, M., Vinyals, O., Li, Y. and Pascanu, R. (2018), ‘Relational inductive biases, deep learning, and graph networks’.
  * Bishop (1995) Bishop, C. (1995), Neural Networks for Pattern Recognition, Oxford University Press, chapter 6, pp. 194–225.
  * Chollet et al. (2015) Chollet, F. et al. (2015), ‘Keras’, https://keras.io.
  * Dozat (2016) Dozat, T. (2016), ‘Incorporating nesterov momentum into adam’, International Conference on Learning Representations 2016 .
  * Duchi et al. (2011) Duchi, J., Hazan, E. and Singer, Y. (2011), ‘Adaptive subgradient methods for online learning and stochastic optimization’, Journal of Machine Learning Research 12(61), 2121–2159.
  * Grudniewski and Sobey (2019) Grudniewski, P. A. and Sobey, A. J. (2019), Do general genetic algorithms provide benefits when solving real problems?, in ‘2019 IEEE Congress on Evolutionary Computation (CEC)’, IEEE, pp. 1822–1829.
  * Grudniewski and Sobey (2021) Grudniewski, P. A. and Sobey, A. J. (2021), ‘cMLSGA: a co-evolutionary multi-level selection genetic algorithm for multi-objective optimization’.
  * Hanson and Burr (1987) Hanson, S. and Burr, D. (1987), Minkowski-r back-propagation: learning in connectionist models with non-euclidian error signals., in ‘Neural Information Processing Systems (NIPS 1987)’.
  * Hinton et al. (2012) Hinton, G., Srivastava, N. and Swersky, K. (2012), ‘Lecture 6a:overview of mini-batch gradient descent’, http://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf.
  * Jeon et al. (2018) Jeon, M., Noh, Y., Shin, Y., Lim, O., Lee, I. and Cho, D. (2018), ‘Prediction of ship fuel consumption by using an artificial neural network’, Journal of Mechanical Science and Technology 32(12), 5785–5796.
  * Jin et al. (2004) Jin, Y., Okabe, T. and Sendhoff, B. (2004), Neural network regularization and ensembling using multi-objective evolutionary algorithms, in ‘Proceedings of the 2004 Congress on Evolutionary Computation (IEEE Cat. No.04TH8753)’, Vol. 1.
  * Kingma and Ba (2014) Kingma, D. P. and Ba, J. (2014), ‘Adam: A method for stochastic optimization’, arXiv preprint arXiv:1412.6980 .
  * Kumar et al. (2021) Kumar, P., Batra, S. and Raman, B. (2021), ‘Deep neural network hyper-parameter tuning through twofold genetic approach’, Soft Computing pp. 1–25.
  * Le et al. (2020) Le, L., Lee, G., Park, K. and Kim, H. (2020), ‘Neural network-based fuel consumption estimation for container ships in korea’, Maritime Policy & Management pp. 1–18.
  * Leifsson et al. (2008) Leifsson, L., Sævarsdóttir, H., Sigurðsson, S. and Vésteinsson, A. (2008), ‘Grey-box modeling of an ocean vessel for operational optimization’, Simulation Modelling Practice and Theory 16(8), 923–932.
  * Liang et al. (2019) Liang, Y., Niu, D. and Hong, W. (2019), ‘Short term load forecasting based on feature extraction and improved general regression neural network model’, Energy 166, 653–663.
  * Luketina et al. (2016) Luketina, J., Berglund, M., Greff, K. and Raiko, T. (2016), ‘Scalable gradient-based tuning of continuous regularization hyperparameters’.
  * Nowlan and Hinton (1992) Nowlan, S. J. and Hinton, G. E. (1992), ‘Simplifying neural networks by soft weight sharing’, Neural Computation .
  * Parkes et al. (2019) Parkes, A. I., Savasta, T. D., Sobey, A. J. and Hudson, D. A. (2019), Efficient vessel power prediction in operational conditions using machine learning., in ‘Practical Design of Ships and Other Floating Structures(PRADS), September 2019, Yokohama, Japan’.
  * Parkes et al. (2018) Parkes, A. I., Sobey, A. J. and Hudson, D. A. (2018), ‘Physics-based shaft power prediction for large merchant ships using neural networks’, Ocean Engineering 166, 92–104.
  * Parkes et al. (2021) Parkes, A. I., Sobey, A. J. and Hudson, D. A. (2021), ‘Towards error measures which influence a learners inductive bias to the ground truth’.
  * Pedersen and Larsen (2009) Pedersen, B. P. and Larsen, J. (2009), Prediction of full-scale propulsion power using artificial neural networks, in ‘Proceedings of the 8th international conference on computer and IT applications in the maritime industries (COMPIT’09), Budapest, Hungary May’, pp. 10–12.
  * Petersen et al. (2012) Petersen, J. P., Jacobsen, D. J. and Winther, O. (2012), ‘Statistical modelling for ship propulsion efficiency’, Journal of marine science and technology 17(1), 30–39.
  * Smith and Jin (2014) Smith, C. and Jin, Y. (2014), ‘Evolutionary multi-objective generation of recurrent neural network ensembles for time series prediction’, Neurocomputing 143, 302–311.
  * Srivastava et al. (2014) Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I. and Salakhutdinov, R. (2014), ‘Dropout: a simple way to prevent neural networks from overfitting’, The Journal of Machine Learning Research 15(1), 1929–1958.
  * Tani et al. (2021) Tani, L., Rand, D., Veelken, C. and Kadastik, M. (2021), ‘Evolutionary algorithms for hyperparameter optimization in machine learning for application in high energy physics’, The European Physical Journal C 81(2), 1–9.
  * Wager et al. (2013) Wager, S., Wang, S. and Liang, P. (2013), ‘Dropout training as adaptive regularization’, arXiv preprint arXiv:1307.1493 .
  * Wang et al. (2019) Wang, B., Sun, Y., Xue, B. and Zhang, M. (2019), Evolving deep neural networks by multi-objective particle swarm optimization for image classification, in ‘Proceedings of the Genetic and Evolutionary Computation Conference’, GECCO ’19, p. 490–498.
  * Willard et al. (2020) Willard, J., Jia, X., Xu, S., Steinbach, M. and Kumar, V. (2020), ‘Integrating physics-based modeling with machine learning: a survey’.
  * Yang et al. (2021) Yang, S., Tian, Y., He, C., Zhang, X., Tan, K. C. and Jin, Y. (2021), ‘A gradient-guided evolutionary approach to training deep neural networks’, IEEE Transactions on Neural Networks and Learning Systems pp. 1–15.
  * Zeiler (2012) Zeiler, M. D. (2012), ‘Adadelta: an adaptive learning rate method’.


# Kink-antikink collision in the supersymmetric $\phi^{4}$ model

João G. F. Campos
Departamento de Física, Universidade Federal da Paraíba,
João Pessoa - PB - 58051-970, Brazil
<EMAIL_ADDRESS>Azadeh Mohammadi
Departamento de Física, Universidade Federal de Pernambuco,
Av. Prof. Moraes Rego, 1235, Recife - PE - 50670-901, Brazil
<EMAIL_ADDRESS>
###### Abstract

This paper investigates a model containing $\phi^{4}$ kinks interacting with
fermions. The fermion back-reaction is included in the equations of motion,
which affects the kink-antikink collisions. We show that the fermion field
generates a force that can be either attractive or repulsive. Moreover, we
investigate three different scenarios, which exhibit a wide variety of
behaviors, including the usual scenarios observed in the $\phi^{4}$ model as
well as the formation of two oscillons, reflection without contact, one-bounce
resonance windows, and the creation of kink-antikink pairs. We also find
evidence that the fermion field can store part of the energy responsible for
the energy exchange mechanism.

##  1 Introduction

Linear and nonlinear wave field equations are ubiquitous in physics, with
applications in almost all research fields. One fascinating feature of the
nonlinear ones is that they may exhibit solitary waves and solitons. The
solitary waves are localized structures that propagate without losing their
shape. Solitons have the additional property that, after a collision, the
initial constituents reappear unchanged, except for a phase shift. Examples of
localized structures are kinks, vortices, skyrmions, and monopoles [1, 2].

Kinks and antikinks, their $x\to-x$ counterparts, are localized solutions that
appear in systems that effectively have one spatial dimension. In higher
dimensions, the kinks can be extended, forming domain walls. They have
applications in many physical systems such as polymers [3], Josephson
junctions [4], DNA [5], deformations in graphene [6] and surface displacement
in shallow liquids [7]. The two prototypical models which describe kinks are
the sine-Gordon and $\phi^{4}$. The former is integrable and exhibits
solitons, while the latter is non-integrable and exhibits solitary waves.
Thus, the behavior of $\phi^{4}$ kinks is much richer because their
interactions are highly inelastic. Since the triplet of seminal works on kink-
antikink collisions [8, 9, 10], it is known that, for non-integrable models, a
kink-antikink collision can exhibit resonance windows. These are intervals of
the initial kink’s relative velocity where they separate after multiple
bounces. The explanation given for this phenomenon is called the resonant
energy exchange mechanism. It states that, after the first collision, the
translational energy is converted into vibrational energy. The vibrational
energy does not dissipate and can be recovered in subsequent collisions,
allowing the system to finally separate. The set of initial velocities where
resonance windows occur forms a fractal structure. It means that the system is
chaotic and consequently has a high sensitivity to initial conditions. Such
complex behavior can only happen for non-integrable models.

A few examples of recent investigations involving kink-antikink interactions
consider wobbling kinks [11, 12], kinks with quasinormal modes [13, 14], long-
range [15, 16, 17, 18] and short-range [19] kinks. Other interesting works
involve hyperbolic models [20, 21, 22], models with Bogomol’nyi-Prasad-
Sommerfield (BPS) impurities [23, 24, 25, 26], interaction with boundaries
[27, 28, 29] and an interpolation between the $\phi^{4}$ and sine-Gordon
models [30]. Moreover, advances in the application of collective coordinates
to kink-antikink collisions were able to give a more quantitative description
of the resonant energy exchange mechanism [31, 32, 33].

Solitary waves interactions for multiple scalar fields have been investigated
in the literature [34, 35, 36, 37]. However, studies involving collisions of
solitary waves composed of both scalar and fermionic fields are still lacking.
The interaction of topological structures with fermions leads to many novel
properties such as fractional quantum numbers [38, 39], and domain wall and
string superconductivity [40, 41]. The system composed of a fermion
interacting with a single kink via a Yukawa coupling can be solved
analytically if the back-reaction is ignored [42, 43]. As discussed below,
these results are important because the back-reaction vanishes at the limit of
a very massive kink. However, in general, the effect of the back-reaction can
be significant. It can be self-consistently included, as shown in Refs. [44,
45] for kinks, and in Ref. [46] for baby-skyrmions. Interestingly, two kinks
or two baby-skyrmions can bind due to the presence of a fermion field [47,
48].

Here, we will consider collisions between a kink and an antikink interacting
with a fermion field, including the fermion back-reaction on the kinks. In
[49, 50] fermions interacting with a kink-antikink configuration have been
investigated, although with a frozen scalar field configuration, i.e., without
evolving in time and without fermion back-reaction. In [51, 52], the authors
studied the kink-antikink collision, but without the back-reaction as well.
More recently, we have studied fermions interacting with another dynamic
scalar field with the same assumption [53, 54]. Therefore, it is essential to
study the effect created by the fermion on the scalar field evolution, which
is lacking in the literature.

In section 2, we introduce briefly a supersymmetric extension of the
$\phi^{4}$ model, which is the simplest way to add a fermion field with
minimum free parameters. We will define symmetric configurations of interest
and compute the inter-fermion force. In section 3, we will briefly discuss the
stability equation of the coupled system. In section 4, we will discuss the
numerical results for kink-antikink collision and the corresponding fermion
field configurations. Finally, we will summarize our findings in the last
section.

##  2 Model

We are interested in the field theory in (1+1) dimensions described by the
following action

$S=\int
d^{2}y\left(\frac{1}{2}\partial^{\mu}\phi\partial_{\mu}\phi+\frac{1}{2}i\bar{\Psi}\gamma^{\mu}\partial_{\mu}\Psi+\frac{1}{2}F^{2}+W_{\phi}F-\frac{1}{2}W_{\phi\phi}\bar{\Psi}\Psi\right),$
(1)

which can be obtained from the minimal supersymmetric model studied in Refs.
[55, 56]. The model is invariant under a supersymmetry transformation which
guarantees the existence of a supercurrent and a supercharge. For more detail
on supersymmetry algebra and supersymmetric solitons, we refer the reader to
the two references cited above.

Substituting the equation of motion for the auxiliary field $F=-W_{\phi}$, and
choosing

$W_{\phi}=\sqrt{\frac{\lambda}{2}}\left(\phi^{2}-\frac{m^{2}}{\lambda}\right),$
(2)

it becomes

$S=\int
d^{2}y~{}\left(\frac{1}{2}\partial^{\mu}\phi\partial_{\mu}\phi+\frac{1}{2}i\bar{\Psi}\gamma^{\mu}\partial_{\mu}\Psi-\sqrt{\frac{\lambda}{2}}\phi\bar{\Psi}\Psi-V(\phi)\right),$
(3)

with $V(\phi)=\frac{1}{2}W_{\phi}^{2}$. The above action consists of the
$\phi^{4}$ model coupled with a Dirac field via a Yukawa coupling, ensuring
the supersymmetric form. The equations of motion obtained from the action
should be accompanied by the normalization condition

$\int_{-\infty}^{\infty}\Psi^{\dagger}(y,\tau)\Psi(y,\tau)dy=1.$ (4)

Now, we perform the usual rescaling to simplify the $\phi^{4}$ potential,
besides working with dimensionless parameters. We define the kink’s value at
infinity $\phi_{0}=m\lambda^{-\frac{1}{2}}$, and set $\phi\to\phi_{0}\chi$ and
$y^{\mu}\to\phi_{0}^{-1}\lambda^{-\frac{1}{2}}x^{\mu}$. To maintain the
normalization condition we also set
$\Psi\to\phi_{0}^{\frac{1}{2}}\lambda^{\frac{1}{4}}\psi$. This way we work
with the following dimensionless lagrangian density

${\cal
L}=\frac{\phi_{0}^{2}}{2}\partial^{\mu}\chi\partial_{\mu}\chi+\frac{1}{2}i\bar{\psi}\gamma^{\mu}\partial_{\mu}\psi-\frac{1}{\sqrt{2}}\chi\bar{\psi}\psi-\frac{\phi_{0}^{2}}{4}\left(\chi^{2}-1\right)^{2}.$
(5)

supplemented by the normalization condition

$\int_{-\infty}^{\infty}\psi^{\dagger}(x,t)\psi(x,t)dx=1.$ (6)

The equations of motion for this model are

$\displaystyle i\gamma^{\mu}\partial_{\mu}\psi-\sqrt{2}\chi\psi=0,$ (7)
$\displaystyle\partial_{\mu}\partial^{\mu}\chi+\chi\left(\chi^{2}-1\right)+\frac{1}{\sqrt{2}\phi_{0}^{2}}\bar{\psi}\psi=0.$
(8)

The first equation is the Dirac equation for a fermion field coupled to the
scalar field $\chi$. The second is the equation of motion for the rescaled
$\phi^{4}$ model plus the term with the fermion field giving rise to the
fermion back-reaction. Looking at the above equations of motion, one could be
tempted to rescale the fermion field to eliminate the parameter $\phi_{0}$.
The problem is that this rescaling would modify the normalization condition.
In fact, it is not difficult to see that $\phi_{0}$ quantifies the strength of
the back-reaction. Notice that the Dirac equation has fixed parameters because
it does not depend on $\phi_{0}$. In the limit when the parameter $\phi_{0}$
is very large, the back-reaction term in eq. (8) vanishes. Therefore, the
approximation where the fermion back-reaction is neglected is justified for a
very large $\phi_{0}$, weak-coupling regime. In other words, a very massive
kink does not feel the fermion. However, for a small $\phi_{0}$, large-
coupling regime, the impact of the fermion back-reaction on the scalar field
can be considerably large. This is the regime we are most interested in.

###  2.1 Bound states

To move forward, we need to choose a representation for Dirac matrices. We
choose $\gamma^{0}=\sigma_{1}$, and $\gamma^{1}=i\sigma_{3}$, where
$\sigma_{i}$ are the Pauli matrices. Writing $\psi=e^{-iEt}(u\;v)^{T}$, with
real $u$ and $v$, besides considering the scalar field to be static, leads to
the coupled equations of motion as follows

$\displaystyle Eu+v^{\prime}-\sqrt{2}\chi v=0,$ (9) $\displaystyle
Ev-u^{\prime}-\sqrt{2}\chi u=0,$ (10)
$\displaystyle\chi^{\prime\prime}-\chi(\chi^{2}-1)-\frac{\sqrt{2}}{\phi_{0}^{2}}uv=0.$
(11)

These equations are supplemented by the normalization condition

$\int dx\left(u^{2}+v^{2}\right)=1.$ (12)

It is instructive to derive the fermion zero mode solution, which is required
by the index theorem. Setting $E=0$, eqs. 9-10 can be solved in the form

$\displaystyle u$
$\displaystyle=a_{1}e^{-\sqrt{2}\int^{x}\chi(x^{\prime})dx^{\prime}}$ (13)
$\displaystyle v$
$\displaystyle=a_{2}e^{\sqrt{2}\int^{x}\chi(x^{\prime})dx^{\prime}}$ (14)

For a kink solution $\chi(x\to\pm\infty)=\pm 1$. Thus, $v$ diverges, unless
$a_{2}=0$. With the same argument, for an antikink, $a_{1}=0$. In both cases,
the product $uv=0$, and the back-reaction vanishes, meaning that the fermion
does not deform the $\phi^{4}$ kink solution

$\chi_{k}(x)=\tanh\left(\frac{x}{\sqrt{2}}\right).$ (15)

This allows one to find the normalized zero mode solution

$\psi_{0}(x)=\sqrt{\frac{3}{4\sqrt{2}}}\begin{pmatrix}\text{sech}^{2}\left(\frac{x}{\sqrt{2}}\right)\\\
0\end{pmatrix}.$ (16)

For the antikink $\chi_{\bar{k}}(x)=-\chi_{k}(x)$, we have instead

$\psi_{\bar{0}}(x)=\sqrt{\frac{3}{4\sqrt{2}}}\begin{pmatrix}0\\\
\text{sech}^{2}\left(\frac{x}{\sqrt{2}}\right)\end{pmatrix}.$ (17)

As both kink and antikink solutions, $\chi_{k}$ and $\chi_{\bar{k}}$, obey the
BPS equation, they preserve half of the supersymmetry [55].

To find other solutions, eqs. (9)-(11) can be solved self-consistently using
relaxation methods, as done in Refs. [44, 45]. Here, we solve the system in
the interval $0.0<x<50.0$ uniformly separated by 2500 points. The center of
the kink is at $x=0$, and the solution for negative $x$ is obtained using the
solution’s parity properties. Then, the derivatives are discretized using a
five-point stencil approximation. The resulting nonlinear system is solved by
utilizing the optimize.root function in Python’s SciPy library.

In the supersymmetric case, besides the zero mode, there are also two other
bound states with positive and negative energy. They are related by the
energy-reflection symmetry when there is no back-reaction. However, the back-
reaction breaks this symmetry. Unlike the fermion zero mode, these solutions
deform the kink and depend on $\phi_{0}$. We denote them by $\chi_{k,\pm
1}(x;\phi_{0})$ and $\psi_{\pm 1}(x,t;\phi_{0})$. There is a phase freedom in
$\psi_{\pm 1}(x,t;\phi_{0})$, which we fix by choosing the real solution whose
second component is positive. We also take $\chi_{\bar{k},\pm
1}(x;\phi_{0})=-\chi_{k,\pm 1}(x;\phi_{0})$ and
$\psi_{\pm\bar{1}}(x,t;\phi_{0})=i\sigma_{2}\psi_{\pm 1}(x,t;\phi_{0})$, where
the bar above the subscript corresponds to an antikink. It is easy to show
that they are also solutions to the equations of motion. The discrete energy
eigenvalues as a function of $\phi_{0}$ are shown in Fig. 1. Both positive and
negative energy eigenvalues are smaller than $\pm\sqrt{3/2}$, respectively.
These constants are the eigenvalues obtained without back-reaction, which are
approached as $\phi_{0}$ increases. We do not plot the positive energy
eigenvalues for $\phi_{0}\lesssim 0.51$ because the numerical method that we
employed to compute them does not reach the desired accuracy111In Ref. [44],
the authors had similar accuracy issues near the same range of $\phi_{0}$
values.. On the other hand, the negative energy state enters the continuum at
$\phi_{0}\simeq 0.59$

Figure 1: Energy eigenvalues (solid curves) as a function of the parameter
$\phi_{0}$. Dashed lines correspond to the threshold of the continuum spectrum
$E=\sqrt{2}$. Dotted lines correspond to the analytical value obtained in the
absence of back-reaction $E=\sqrt{3/2}$. This figure is a reproduction of
results obtained in [44].

To study the system’s dynamics, we also need boosted solutions, which can be
obtained in the usual way. For the fermion field
$\psi=(\psi_{1}\;\psi_{2})^{T}$, a boost with velocity $v_{i}$ reads

$\displaystyle\psi_{1}^{\prime}(x,t)$
$\displaystyle=\cosh(\nu/2)\psi_{1}(x^{\prime},t^{\prime})-i\sinh(\nu/2)\psi_{2}(x^{\prime},t^{\prime})$
(18) $\displaystyle\psi_{2}^{\prime}(x,t)$
$\displaystyle=i\sinh(\nu/2)\psi_{1}(x^{\prime},t^{\prime})+\cosh(\nu/2)\psi_{2}(x^{\prime},t^{\prime}),$
(19)

where we defined the rapidity $\nu=\tanh^{-1}(v_{i})$, and the transformed
coordinates $x^{\prime}=x-v_{i}t$ and $t^{\prime}=t-v_{i}x$. A boosted
solution will be identified by an explicit dependence on $t$ and $v_{i}$. For
example, $\chi_{k}(x,t;v_{i})$ denotes a boosted kink.

###  2.2 Symmetric collisions

First, let us focus on symmetric collisions where before the kinks superpose,
the fields can be approximately described by

$\displaystyle\chi(x,t)$
$\displaystyle=\chi_{k}(x+x_{0};v_{i})-\chi_{k}(x-x_{0};-v_{i})-1,$ (20)
$\displaystyle\psi(x,t)$
$\displaystyle=\frac{1}{\sqrt{2}}\psi_{0}(x+x_{0};v_{i})+\frac{A}{\sqrt{2}}\psi_{\bar{0}}(x-x_{0};-v_{i}),$
(21)

where $A$ is a complex phase in general and $\psi_{0}(\psi_{\bar{0}})$ is the
fermion zero mode bound to a kink(antikink). We plan to find the value of $A$
for which the fermion will have the largest impact on the scalar field for a
fixed value of $\phi_{0}$. To simplify our analysis we fix $v_{i}=0$. Then

$\psi(x,t)=\sqrt{\frac{3}{8\sqrt{2}}}\begin{pmatrix}\text{sech}^{2}\left(\frac{x+x_{0}}{\sqrt{2}}\right)\\\
A\,\text{sech}^{2}\left(\frac{x-x_{0}}{\sqrt{2}}\right)\end{pmatrix}.$ (22)

The back-reaction is proportional to

$\bar{\psi}(x,t)\psi(x,t)=\frac{3}{8\sqrt{2}}\text{sech}^{2}\left(\frac{x+x_{0}}{\sqrt{2}}\right)\text{sech}^{2}\left(\frac{x-x_{0}}{\sqrt{2}}\right)(A+A^{*}).$
(23)

The absolute value of this expression is largest when $A=\pm 1$. We will study
collisions for these two values of $A$.

The second example is the symmetric configuration where the fermion is in the
excited state with positive energy. It reads

$\displaystyle\chi(x,t)$
$\displaystyle=\chi_{k,+1}(x+x_{0};\sqrt{2}\phi_{0},v_{i})-\chi_{k,+1}(x-x_{0};\sqrt{2}\phi_{0},-v_{i})-1,$
(24) $\displaystyle\psi(x,t)$
$\displaystyle=\frac{1}{\sqrt{2}}\psi_{+1}(x+x_{0},t;\sqrt{2}\phi_{0},v_{i})+\frac{A}{\sqrt{2}}\psi_{+\bar{1}}(x-x_{0},t;\sqrt{2}\phi_{0},-v_{i}).$
(25)

It is necessary to multiply the parameter $\phi_{0}$ by $\sqrt{2}$ to
compensate for the reduction in the strength of the back-reaction. Fixing
$v_{i}=0$ and writing
$\psi_{+1}(x,t;\sqrt{2}\phi_{0})=e^{-iEt}(u(x)\;v(x))^{T}$, the back-reaction
is proportional to

$\bar{\psi}(x,t)\psi(x,t)=\mathcal{B}_{1}+\mathcal{B}_{2},$ (26)

with

$\displaystyle\mathcal{B}_{1}=(u(x+x_{0})v(x+x_{0})-u(x-x_{0})v(x-x_{0})),$
(27)
$\displaystyle\mathcal{B}_{2}=\frac{(A+A^{*})}{2}(v(x-x_{0})v(x+x_{0})-u(x-x_{0})u(x+x_{0})).$
(28)

The first term quantifies the fermion back-reaction on the kink at the same
position. On the other hand, the second one describes the back-reaction due to
the overlap, which depends on the $A$ choice. Again, the fermion back-reaction
is maximized for $A=\pm 1$. In what follows, for brevity, we will analyze only
the positive sign for the excited states. None of the kink-antikink
configurations discussed above is BPS and, therefore, breaks supersymmetry.

###  2.3 Energy-momentum tensor and force

Due to the translational invariance of the model, there are four conserved
currents described by the energy-momentum tensor. One can show that the tensor
is given by

$T^{\mu\nu}=T^{\mu\nu}_{\text{scalar}}+T^{\mu\nu}_{\text{Dirac}},$ (29)

where

$\frac{T^{\mu\nu}_{\text{scalar}}}{\phi_{0}^{2}}=\partial^{\mu}\chi\partial^{\nu}\chi-\eta^{\mu\nu}\left(\frac{1}{2}\partial^{\sigma}\chi\partial_{\sigma}\chi-\frac{1}{4}(\chi^{2}-1)^{2}\right),$
(30)

and

$T^{\mu\nu}_{\text{Dirac}}=\frac{i}{4}\left[\bar{\psi}\gamma^{\mu}\left(\partial^{\nu}\psi\right)-\left(\partial^{\nu}\bar{\psi}\right)\gamma^{\mu}\psi\right].$
(31)

Suppose that we have a kink-antikink configuration where the two are
symmetrically placed around the origin and bound to fermions. Examples are
given in eqs. (20) and (21) as well as in eqs. (24) and (25). We can easily
compute the force between the two composite systems for such cases. The
momentum on the left side is

$P=\int_{-\infty}^{0}dxT^{01}(x),$ (32)

and the force is

$F=\frac{dP}{dt}=\int_{-\infty}^{0}dx\partial_{t}T^{01}(x)=-\int_{-\infty}^{0}dx\partial_{x}T^{11}(x)=-T^{11}(0).$
(33)

If we consider eqs. (20) and (21) with $v_{i}=0$ and $A=\pm 1$, the scalar
field contribution is the usual one $F_{\text{scalar}}\simeq
16\phi_{0}^{2}e^{-2\sqrt{2}x_{0}}$, while the fermion field causes a force
$F_{\text{Dirac}}\simeq\pm 6e^{-2\sqrt{2}x_{0}}$, respectively. Therefore,
$A=+1$ corresponds to the attractive case and $A=-1$ to the repulsive one. The
total force is $F=F_{\text{scalar}}+F_{\text{Dirac}}$. One can see that, the
fermion force surpasses the scalar field one as the parameter $\phi_{0}$
decreases.

We can also compute the force for the configuration described by eqs. (24) and
(25) with $v_{i}=0$ and $A=+1$. We write the excited state as
$\psi_{+1}(x,t;\phi_{0})=e^{-iEt}(u(x)\;v(x))^{T}$. It is easy to show that
the asymptotic form of the state for $x\gg 1$ is

$\psi_{+1}(x,t;\phi_{0})\simeq
e^{-iEt}\begin{pmatrix}\frac{1}{E}(\sqrt{2}+\kappa)v(x)\\\ v(x)\end{pmatrix},$
(34)

with $v(x)=Ce^{-\kappa x}$ knowing that $\kappa\equiv\sqrt{2-E^{2}}$ and the
constant $C$ is real. Besides that, the parity symmetry implies $v(x)=v(-x)$
and $u(x)=-u(-x)$. Therefore, the force can be approximated by

$F=\frac{1}{2}[u(x_{0})+v(x_{0})][-u^{\prime}(x_{0})+v^{\prime}(x_{0})]\simeq
C^{2}\,\frac{\kappa^{2}(\sqrt{2}+\kappa)}{E^{2}}\,e^{-2\kappa x_{0}},$ (35)

Therefore, we find that the force is again attractive.

##  3 Stability equation

It is also important to look at the stability of the fermion-kink system. To
do so, we write $\chi(x,t)=\chi_{k}(x)+\eta(x,t)$,
$\dot{\chi}(x,t)=\dot{\eta}(x,t)\equiv\zeta(x,t)$ and

$\psi(x,t)=e^{-iEt}\left[\begin{pmatrix}u(x)\\\
v(x)\end{pmatrix}+\begin{pmatrix}\alpha_{1}(x,t)+i\alpha_{2}(x,t)\\\
\beta_{1}(x,t)+i\beta_{2}(x,t)\end{pmatrix}\right],$ (36)

where $\eta$, $\zeta$, $\alpha_{1}$, $\alpha_{2}$, $\beta_{1}$ and $\beta_{2}$
are small and real. The second column matrix in eq. (36) will be denoted by
$\delta\psi$. Substituting the above expressions in the equations of motion
leads to the following linear system

$\begin{pmatrix}\dot{\eta}\\\ \dot{\zeta}\\\ \dot{\alpha}_{1}\\\
\dot{\alpha}_{2}\\\ \dot{\beta}_{1}\\\
\dot{\beta}_{2}\end{pmatrix}=\begin{pmatrix}0&1&0&0&0&0\\\
\frac{d^{2}}{dx^{2}}-V^{\prime\prime}(\chi_{k})&0&-\frac{\sqrt{2}}{\phi_{0}^{2}}v&0&-\frac{\sqrt{2}}{\phi_{0}^{2}}u&0\\\
0&0&0&-E&0&-\frac{d}{dx}+\sqrt{2}\chi_{k}\\\
-\sqrt{2}v&0&E&0&\frac{d}{dx}-\sqrt{2}\chi_{k}&0\\\
0&0&0&\frac{d}{dx}+\sqrt{2}\chi_{k}&0&-E\\\
-\sqrt{2}u&0&-\frac{d}{dx}-\sqrt{2}\chi_{k}&0&E&0\end{pmatrix}\begin{pmatrix}\eta\\\
\zeta\\\ \alpha_{1}\\\ \alpha_{2}\\\ \beta_{1}\\\ \beta_{2}\end{pmatrix}$ (37)

The solution is stable if there are no eigenvalues $\lambda$ with a positive
real part. Moreover, an oscillatory solution is purely imaginary and comes
with its complex conjugate. Therefore, the number of modes is half the number
of eigenvalues. To find the spectrum, we discretize the derivatives using a
five-point stencil approximation, and the resulting matrix is diagonalized
using the NumPy library in Python. This method has a few difficulties that
make the analysis very challenging. First, the operator is not Hermitian. This
means that the resulting matrix is not symmetric, and its eigenvectors are not
orthogonal. Thus, they cannot be interpreted as independent internal modes.
The other issues will be explained below.

For the fermions at the zero mode, $E=0$ and

$\begin{pmatrix}u\\\ v\end{pmatrix}=\psi_{0},$ (38)

one can find two modes with $\lambda=0$ analytically. The first one is given
by the derivative of the unperturbed solution with respect to $x$

$\displaystyle\eta_{0,s}=\mathcal{N}\text{sech}^{2}(x/\sqrt{2}),\quad\zeta_{0,s}=0,$
(39)
$\displaystyle\delta\psi_{0,s}=\mathcal{N}\sqrt{\frac{3}{4\sqrt{2}}}\begin{pmatrix}-2\,\text{sech}^{2}(x/\sqrt{2})\tanh(x/\sqrt{2})\\\
0\end{pmatrix}.$ (40)

This solution is guaranteed to solve the linearized system by the
translational symmetry of the model. From this construction, it should be
clear that this mode receives no back-reaction, as well. The second mode with
$\lambda=0$ is easily found in the form

$\eta_{0,f}=0,\quad\zeta_{0,f}=0,\quad\delta\psi_{0,f}=\psi_{0}.$ (41)

Then, the other modes can be computed numerically. We found three vibrational
modes with purely imaginary eigenvalues and no $\phi_{0}$ dependence. Thus,
all three modes have vanishing back-reaction with absolute value
$|\lambda|=\sqrt{3/2}$ within our numerical precision. We did not find any
vibrational mode with non-vanishing back-reaction using this numerical method.
The spectrum of the linear matrix is closely related to the combined spectrum
of scalar perturbations around an isolated kink and fermionic bound states on
a kink background without back-reaction. In this case, both excitation
frequencies are equal as the model is supersymmetric.

In general, our numerical analysis shows that the composite system is stable
because we did not find any eigenvalues with a positive real part, even for a
fermion in the excited state. However, we avoid a more detailed stability
analysis of the latter case. The main issue is that the excited state is close
to the threshold between discrete and continuum eigenvalues, making it
difficult to distinguish between the discrete vibrational modes and continuum
states.

##  4 Results

###  4.1 Dynamics

To evolve the full equations of motion in time, we need to define the real and
imaginary components of the fermionic field as follows
$\psi=(u_{1}+iu_{2}\;v_{1}+iv_{2})^{T}$. Thus, the equations of motion, eqs.
(7) and (8), result in

$\displaystyle\ddot{\chi}$
$\displaystyle=\partial_{x}^{2}\chi-\chi(\chi^{2}-1)-\frac{\sqrt{2}}{\phi_{0}^{2}}(u_{1}v_{1}+u_{2}v_{2}),$
(42) $\displaystyle\dot{u}_{1}$ $\displaystyle=-\partial_{x}v_{2}+\sqrt{2}\chi
v_{2},$ (43) $\displaystyle\dot{u}_{2}$
$\displaystyle=\partial_{x}v_{1}-\sqrt{2}\chi v_{1},$ (44)
$\displaystyle\dot{v}_{1}$ $\displaystyle=\partial_{x}u_{2}+\sqrt{2}\chi
u_{2},$ (45) $\displaystyle\dot{v}_{2}$
$\displaystyle=-\partial_{x}u_{1}-\sqrt{2}\chi u_{1}.$ (46)

Clearly, if $u_{2}(x,t=0)=v_{1}(x,t=0)=0$, then both fields will vanish at all
times and consequently the back-reaction will be zero at all times. Similar
considerations apply when $u_{1}(x,t=0)=v_{2}(x,t=0)=0$.

Suppose that we initialize the scalar field at a boosted kink-antikink
configuration

$\chi(x,t)=\chi_{k}(x+x_{0},t;v)-\chi_{k}(x-x_{0},t;-v)-1.$ (47)

Moreover, we consider a single fermion at the zero mode bound to the kink,
$\psi(x,t)=\psi_{0}(x+x_{0};v)$. This configuration obeys
$u_{2}(x,t=0)=v_{1}(x,t=0)=0$ and the back-reaction vanishes at all times. In
order to turn on the back-reaction, we can consider the symmetric
configurations discussed in section 2.2. These are the attractive and
repulsive cases described by eq. (21) with $A=\pm 1$, and the case with the
excited fermion is described by eq. (25) with $A=+1$. The latter is also
attractive.

We integrate the equations of motion in a box in the interval $-100.0<x<100.0$
by discretizing space in intervals with width $\Delta x=0.05$. The partial
derivatives with respect to $x$ were approximated using a five-point stencil.
The resulting set of equations was integrated using a fifth-order Runge-Kutta
method with adaptive step size. Moreover, we considered periodic boundary
conditions. When $\phi_{0}$ is very small, the error in the energy
conservation may increase due to the increased strength of the nonlinear back-
reaction term. The error tolerance in energy conservation in our numerical
simulation is below $10^{-4}$. However, for large $\phi_{0}$, the error is
much less significant.

###  4.2 Fermions at the zero mode

In this section, we consider the configuration described by eqs. (20) and
(21). The results for the attractive and repulsive cases are summarized in
Figs. 2 and 3, respectively. The colors represent the value of the scalar
field at the origin and $t=40/v_{i}$. This value is plotted as a function of
the velocity $v_{i}$ and the parameter $\phi_{0}$. Particular cases of the
field’s evolution in spacetime are shown in Figs. 4 and 5.

Figure 2: Final value of the scalar field at the center of the collision as a
function of $\phi_{0}$ and $v$. We consider the attractive case where the
fermions start at the zero mode. Figure 3: Final value of the scalar field at
the center of the collision as a function of $\phi_{0}$ and $v$. We consider
the repulsive case where the fermions start at the zero mode.

The kink and the antikink interact and separate in the blue region in Figs. 2
and 3, leaving the vacuum $\phi\simeq+1.0$ at the center. This happens in
three scenarios: resonance windows, reflection, and repulsion without contact.
Resonance windows and reflection happen in both repulsive and attractive
cases. Fig. 4(b) gives an example for two-bounce resonance windows for the
attractive case. A three-bounce example for the repulsive case is shown in
Fig. 5(c). A reflection between the kink and the antikink is shown in Fig.
4(c) for the attractive case.

(a) $\phi_{0}=0.8$, $v_{i}=0.2$

(b) $\phi_{0}=0.8$, $v_{i}=0.235$

(c) $\phi_{0}=0.8$, $v_{i}=0.32$

(d) $\phi_{0}=0.3$, $v_{i}=0.3$

(e) $\phi_{0}=0.475$, $v_{i}=0.3$

(f) $\phi_{0}=0.545$, $v_{i}=0.3$

Figure 4: Scalar field evolution in spacetime. We consider the attractive case
where the fermions start at the zero mode.

(a) $\phi_{0}=0.4$, $v_{i}=0.3$

(b) $\phi_{0}=0.54$, $v_{i}=0.47$

(c) $\phi_{0}=1.083$, $v_{i}=0.2005$

Figure 5: Scalar field evolution in spacetime. We consider the repulsive case
where the fermions start at the zero mode.

For $\phi_{0}=1.2$, the structure of two-bounce resonance windows is similar
to that in the absence of fermions. This is true in both scenarios. As
$\phi_{0}$ decreases and the back-reaction becomes more significant, the
resonance structure is deformed. In both cases, the critical velocity
increases, and the windows become thinner. Moreover, a new window appears
approximately in the interval $0.7\lesssim\phi_{0}\lesssim 0.8$ for the
attractive case. For the repulsive case, it appears for a slightly larger
$\phi_{0}$. This window has one less oscillation between bounces than the
lowest two-bounce resonance windows in the original $\phi^{4}$ model. An
example of the evolution of the field in such a window is shown in Fig. 4(f).
There are thinner windows near the large one, corresponding to separation
after three bounces.

Reflection without contact between the kink and the antikink happens only in
the repulsive case for $\phi_{0}\lesssim 0.6$. The repulsive force is strong
enough to cause the kinks to recede before they collide in this region. An
example is shown in Fig. 5(a). The region where this behavior occurs merges
with the reflection region in the upper right corner of Fig. 3. In between, we
have a very soft bounce which is enough to reflect the kinks due to the
repulsion of the fermionic field. An example is shown in Fig. 5(b). The yellow
and green regions correspond to bion formation. It happens in both scenarios.
An example is shown in Fig. 4(a). Finally, the kinks tend to form two
oscillons after the collision when $\phi_{0}$ decreases in the attractive
case. It leaves a vacuum $\phi\simeq-1.0$ at the center and corresponds to the
red regions in Fig. 2. These regions are called pseudo-windows. There is a
large pseudo-window on the left side of Fig. 2. The behavior in this region is
illustrated in Fig. 4(d). Moreover, there are also many thinner pseudo-
windows. They correspond to the formation of oscillons after more than a
single bounce. We illustrate this behavior in Fig. 4(e).

In both Figs. 2 and 3, one can see that the critical velocity increases as
$\phi_{0}$ is lowered. This means that the kinks need higher translational
energy to overcome the mutual attraction after the collision. There are two
possible explanations for this result. The first is that radiative losses
increase as the back-reaction becomes more significant. The second one is that
fermions are excited and store part of the energy. Unfortunately, it is hard
to isolate the amount of energy lost as radiation and tell which effect is
more important.

We also analyzed the evolution of the fermionic field during the collision.
This is illustrated in Fig. 6 for the same cases discussed in Fig. 4. In
general, after the collision, the fermion stays bound to the kink, the bion,
or the oscillon. Moreover, the fermion density exhibits oscillations
afterward. They are expected for two reasons. First, we know that the scalar
field internal modes are excited after the collision. Thus, the fermion
interacts with an oscillating field [54]. Second, fermion internal modes are
also excited after the collision. However, the oscillations do not have a
small amplitude and cannot be described as a linear combination of vibrational
modes due to the highly nonlinear character of the equations of motion.

(a) $\phi_{0}=0.8$, $v_{i}=0.2$

(b) $\phi_{0}=0.8$, $v_{i}=0.235$

(c) $\phi_{0}=0.8$, $v_{i}=0.32$

(d) $\phi_{0}=0.3$, $v_{i}=0.3$

(e) $\phi_{0}=0.475$, $v_{i}=0.3$

(f) $\phi_{0}=0.545$, $v_{i}=0.3$

Figure 6: Evolution of the fermion density in spacetime. We consider the
attractive case where the fermions start at the zero mode.

The fractal character of bounce windows is an important property of resonant
systems. To analyze further how this property is affected, we fix
$\phi_{0}=0.8$. For this value, the back-reaction is already significant.
Then, we measure the field at the collision center as a function of time and
initial velocity. This is shown in Fig. 7 (a). The instants where bounces
occur can be easily spotted because they occur when the field crosses the
value $-1$, marked by a red color. At bounce windows, the kinks separate, and
the field at the center approaches the value $1$, marked in blue. Therefore,
the resonance windows occur for values of $v_{i}$ where the field eventually
reaches a constant blue color, creating a vertical stripe. In Fig. 7 (b), we
measure the final velocity $v_{f}$ as a function of the initial one. We
searched for bounce windows in the range $0.1<v_{i}<0.5$ with velocity
increments $\Delta v=10^{-5}$. We did the simulations up to a maximum time of
$t=1400.0$. Whenever $v_{f}$ is non-vanishing, there is a bounce window. If
$v_{f}$ is very small, the instant the separation occurs is large, and there
is small imprecision in our measurement due to interaction with returning
emitted radiation. The color indicates the number of bounces before
separation. Interestingly, it is possible that $v_{f}>v_{i}$ in this model
because the fermion can acquire negative energy and, thus, increase the kink’s
energy thanks to energy conservation. In both panels of Fig. 7, we see that
the structure of higher bounce windows is still very rich. Therefore, the
inclusion of back-reaction in our model does not seem to destroy the fractal
property. However, as shown in Figs. 2 and 3, if $\phi_{0}$ is decreased
further, all resonance windows eventually disappear for both attractive and
repulsive cases.

(a)

(b)

Figure 7: (a) Field at the collision center as a function of time and initial
velocity. (b) Kink’s final velocity as a function of the initial one. Dashed
line corresponds to the identity curve. We fix $\phi_{0}=0.8$.

An interesting observation is that when we calculate the scalar field and the
fermion density at the center of the collision, both have the same number of
oscillations between bounces, as shown in Fig. 8. This means that the fermion
and kink oscillate at approximately the same frequency. This is expected when
the fermions start at the zero-mode due to the degenerate excitation spectrum
of the system, calculated in section 3.

(a)

(b)

Figure 8: Evolution of (a) scalar field and (b) the fermion density at the
center of the collision. We consider the attractive case where the fermions
start at the zero mode. Parameters are $v_{i}=0.235$ and $\phi_{0}=0.8$.

###  4.3 Fermions at the excited state

Now, let us turn our attention to the case where the fermions are excited.
There are many possibilities to explore in this case, but we will study only a
single case for brevity. It is given by eq. (25) with $A=1$. The possible
behaviors of the system are shown in Fig. 9. Once more, the colors represent
the value of the scalar field at the origin and $t=40/v_{i}$. The blue region
corresponds to reflection or resonance windows with one or multiple bounces.
The yellow and green regions correspond to bion formation. Finally, the red
region corresponds to pseudo-windows, where two oscillons or an extra kink-
antikink pair are formed.

Figure 9: Final value of the scalar field at the center of the collision as a
function of $\phi_{0}$ and $v_{i}$. We consider the case where the fermions
start at the excited state.

For large $\phi_{0}$, we see a sequence of two-bounce windows at the boundary
between the reflection and annihilation regions. As $\phi_{0}$ is lowered,
some resonance windows are suppressed, and others survive. For instance, for
$\phi_{0}=0.76$ and $v_{i}=0.14$ the resonance window is shown in Fig. 10(a).
This corresponds to the third two-bounce window relative to the original
$\phi^{4}$ model. We associate the suppression with the fact that the fermion
field does not have a well-defined oscillation between bounces, unlike the one
shown in Fig. 8 for example.

(a) $\phi_{0}=0.76$, $v_{i}=0.14$

(b) $\phi_{0}=0.55$, $v_{i}=0.23$

(c) $\phi_{0}=0.575$, $v_{i}=0.23$

(d) $\phi_{0}=0.47$, $v_{i}=0.24$

(e) $\phi_{0}=0.47$, $v_{i}=0.31$

(f) $\phi_{0}=0.425$, $v_{i}=0.43$

Figure 10: Scalar field evolution in spacetime. We consider the case where the
fermions start at the excited state.

On the far left of Fig. 9, pseudo-windows and resonance windows alternate. One
finds two almost vertical stripes between this region and the yellow and green
ones. One stripe is red, and the other is blue, corresponding to two-bounce
pseudo- and resonance windows, respectively, with no oscillations between
bounces. We illustrate this behavior in Figs. 10(b) and (c). On the left of
the stripes, we have an alternating structure of one-bounce pseudo- and
resonance windows. They are illustrated in Figs. 10(d) and (e). The boundary
between the two behaviors approximately forms a spine structure [12], which
can be related to the oscillating character of the fermion at the excited
state. This also suggests that the fermion is storing energy that can be
recovered by the kinks.

The appearance of one-bounce resonance windows is very similar to what is
observed in collisions between wobbling kinks [11]. However, the two results
are different. In the present work, the one-bounce resonance windows are
separated by pseudo-windows, while for wobbling kinks, they are separated by
bions and higher-bounce resonance windows. Even though the model we are
dealing with here is supersymmetric and, in principle the fermionic
excitations could be mapped into the bosonic ones by supersymmetry, the kink-
antikink configurations break this symmetry leaving the results quite distinct
from wobbling kinks’ collisions.

Finally, at the tip of the red region, we find kink-antikink reflection with
the formation of an extra kink-antikink pair. The field evolution, in this
case, is shown in Fig. 10(f) as an example. The energy source to create the
extra pair comes from the fermion, which becomes more energetic relative to
the kink for small $\phi_{0}$. Let us compute the initial energy contained in
the kink-antikink pair to prove this assertion. The energy of a single kink
without back-reaction is $E_{cl}=\frac{2\sqrt{2}}{3}\phi_{0}^{2}$ [40].
Including a boost and a constant $C$, which is different from unity due to the
fermion back-reaction, we obtain

$E_{2}=C\frac{4\sqrt{2}}{3\sqrt{1-v_{i}^{2}}}\phi_{0}^{2}.$ (48)

for the kink-antikink pair with $C\approx 1.166$ which we find numerically.
Plugging in $v_{i}=0.43$ and $\phi_{0}=0.425$, we obtain $E_{2}\approx 0.440$,
which is much less than the minimum energy to create two kink-antikink pairs
at rest $E_{4,min}=\frac{8\sqrt{2}}{3}\phi_{0}^{2}\approx 0.681$.

The behavior of fermions during collisions is very similar to what we found in
the previous cases. After the collision, a significant fraction of the fermion
density stays bound to the kinks, the oscillons, or the bion. Moreover,
oscillations with a nonlinear character are also observed for small
$\phi_{0}$.

##  5 Conclusion

In this work, we studied a supersymmetric model containing fermion-kink
interactions. We showed that the equations of motion contain fermion back-
reaction to the kink. The back-reaction becomes more relevant for light kinks,
which possess a small value at infinity, $\phi_{0}$.

We computed both zero-energy and excited fermion bound states and used boosted
configurations to construct symmetric collisions. We have shown that two
parameters could control the strength of the fermion back-reaction on the
kinks in a symmetric collision, the parameter $\phi_{0}$ and the phase between
the fermions on the kink and antikink $A$. We focused on cases with maximum
back-reaction, $A=\pm 1$ for a fixed value of $\phi_{0}$. Moreover, we have
estimated the force between the fermions computing it from the energy-momentum
tensor. The force can be either attractive or repulsive. As $\phi_{0}$ is
lowered, the interaction between the fermions can become larger than the
scalar field attraction.

Then, we studied the stability equation of the kink-fermion system. The
stability equation is given by a complicated linear operator, which needs to
be analyzed carefully. We employed a simple discretization method and then
diagonalized the matrix numerically. For the composite system with the fermion
at the zero mode, we found a spectrum containing two modes with eigenvalue
$\lambda=0$, which can be found analytically, and degenerate excited states.
Moreover, the analysis suggests that the composite system is stable in all
cases because we did not find positive and real eigenvalues in our numerical
analysis.

Next, we analyzed the equations of motion and discussed the necessary
condition to have a non-vanishing back-reaction. We integrated the equations
of motion in time and found a wide variety of behaviors expected by the back-
reaction’s nonlinear character. For the attractive case with the fermion at
the zero mode, we observed bion formation, resonance windows, reflection, and
pseudo-windows, where two oscillons are formed. We also observed reflection
without contact or with a soft collision for the repulsive case. In both
cases, we found a rich structure of resonance windows up to very low values of
$\phi_{0}$.

After colliding, a significant fraction of the fermion density is bound at the
localized structures in the scalar field. Moreover, the density shows
complicated oscillation patterns that a linear combination of vibrational
modes cannot describe. Interestingly, we observed that the fermions exhibit
the same number of oscillations between bounces as the scalar field. However,
this only occurs for the zero-mode initial condition.

Many fascinating behaviors appeared in the collision between composite systems
with excited fermions. First, some resonance windows are suppressed as
$\phi_{0}$ is lowered. This seems to be related to fact that the resonant
behavior of the fermion field shown in Fig. 8 does not occur anymore. For even
lower values of $\phi_{0}$, we found an alternating structure of pseudo-
windows and one-bounce resonance windows related to the oscillating character
of the excited fermion states. This is similar to what happens in collisions
between wobbling kinks [11, 12], but in this case, the oscillation is
performed by the fermion. Therefore, it appears that the fermion is storing
part of the resonant energy, but a more detailed analysis is needed, as will
be discussed below.

The resonant energy exchange mechanism is clearly in play because resonance
windows are observed. For large $\phi_{0}$, the exchange is conducted by the
kink’s shape mode. The reason is that the back-reaction is negligible, and we
recover the original $\phi^{4}$ in this limit. When $\phi_{0}$ is lowered, it
is natural to ask if the fermion excitation also conducts the energy exchange.
Even though we found evidence that this is the case, we cannot give a definite
answer in the present work because both fermion and kinks excitations have the
same frequency, and it is not possible to distinguish between them. It would
also be interesting to investigate whether a fermion excited state alone could
create a resonant structure. This and other pending issues discussed in this
section are the subjects of an ongoing investigation.

## Acknowledgments

We acknowledge financial support from the Brazilian agencies CAPES and the
National Council for Scientific and Technological Development - CNPq. AM also
thanks financial support from Universidade Federal de Pernambuco Edital Qualis
A. We thank Mauro Copelli for the access to his lab’s computer cluster, which
was essential for conducting the research reported in this paper. Part of the
simulations was also conducted in the SDumont cluster from the Brazilian
laboratory LNCC (Laboratório Nacional de Computação Científica).

## References

  * [1] Ramamurti Rajaraman.  Solitons and instantons.  North Holland, 1982.
  * [2] Nicholas Manton and Paul Sutcliffe.  Topological solitons.  Cambridge University Press, 2004.
  * [3] WP Su, JR Schrieffer, and AJ Heeger.  Solitons in polyacetylene.  Physical review letters, 42(25):1698, 1979.
  * [4] AV Ustinov.  Solitons in josephson junctions.  Physica D: Nonlinear Phenomena, 123(1-4):315–329, 1998.
  * [5] Ludmila V Yakushevich.  Nonlinear physics of DNA.  John Wiley & Sons, 2006.
  * [6] RD Yamaletdinov, VA Slipko, and YV Pershin.  Kinks and antikinks of buckled graphene: A testing ground for the $\varphi^{4}$ field model.  Physical Review B, 96(9):094306, 2017.
  * [7] Bruce Denardo, William Wright, Seth Putterman, and Andrés Larraza.  Observation of a kink soliton on the surface of a liquid.  Physical review letters, 64(13):1518, 1990.
  * [8] David K Campbell, Jonathan F Schonfeld, and Charles A Wingate.  Resonance structure in kink-antikink interactions in $\varphi^{4}$ theory.  Physica D: Nonlinear Phenomena, 9(1-2):1–32, 1983.
  * [9] Michel Peyrard and David K Campbell.  Kink-antikink interactions in a modified sine-gordon model.  Physica D: Nonlinear Phenomena, 9(1-2):33–51, 1983.
  * [10] David K Campbell, Michel Peyrard, and Pasquale Sodano.  Kink-antikink interactions in the double sine-gordon equation.  Physica D: Nonlinear Phenomena, 19(2):165–205, 1986.
  * [11] A Alonso Izquierdo, J Queiroga-Nunes, and LM Nieto.  Scattering between wobbling kinks.  Physical Review D, 103(4):045003, 2021.
  * [12] João GF Campos and Azadeh Mohammadi.  Wobbling double sine-gordon kinks.  Journal of High Energy Physics, 2021(9):67, 2021.
  * [13] Patrick Dorey and Tomasz Romańczukiewicz.  Resonant kink–antikink scattering through quasinormal modes.  Physics Letters B, 779:117–123, 2018.
  * [14] João GF Campos and Azadeh Mohammadi.  Quasinormal modes in kink excitations and kink–antikink interactions: a toy model.  The European Physical Journal C, 80:1–15, 2020.
  * [15] Ivan C Christov, Robert J Decker, A Demirkaya, Vakhid A Gani, PG Kevrekidis, and RV Radomskiy.  Long-range interactions of kinks.  Physical Review D, 99(1):016010, 2019.
  * [16] Ivan C Christov, Robert J Decker, A Demirkaya, Vakhid A Gani, PG Kevrekidis, Avinash Khare, and Avadh Saxena.  Kink-kink and kink-antikink interactions with long-range tails.  Physical review letters, 122(17):171601, 2019.
  * [17] Ivan C Christov, Robert J Decker, A Demirkaya, Vakhid A Gani, PG Kevrekidis, and Avadh Saxena.  Kink-antikink collisions and multi-bounce resonance windows in higher-order field theories.  Communications in Nonlinear Science and Numerical Simulation, 97:105748, 2021.
  * [18] João GF Campos and Azadeh Mohammadi.  Interaction between kinks and antikinks with double long-range tails.  Physics Letters B, 818:136361, 2021.
  * [19] D Bazeia, Adalto R Gomes, and Fabiano C Simas.  Semi-compactness and multiple oscillating pulses in kink scattering.  The European Physical Journal C, 81(6):1–12, 2021.
  * [20] Dionisio Bazeia, Ekaterina Belendryasova, and Vakhid A Gani.  Scattering of kinks of the sinh-deformed $\varphi^{4}$ model.  The European Physical Journal C, 78:1–14, 2018.
  * [21] D Bazeia, Adalto R Gomes, KZ Nobrega, and Fabiano C Simas.  Kink scattering in a hybrid model.  Physics Letters B, 793:26–32, 2019.
  * [22] D Bazeia, Adalto R Gomes, KZ Nobrega, and Fabiano C Simas.  Oscillons in hyperbolic models.  Physics Letters B, 803:135291, 2020.
  * [23] C Adam, K Oles, T Romanczukiewicz, and A Wereszczynski.  Spectral walls in soliton collisions.  Physical review letters, 122(24):241601, 2019.
  * [24] Christoph Adam, K Oles, JM Queiruga, T Romanczukiewicz, and Andrzej Wereszczynski.  Solvable self-dual impurity models.  Journal of High Energy Physics, 2019(7):1–25, 2019.
  * [25] C Adam, T Romanczukiewicz, and A Wereszczynski.  The $\phi^{4}$ model with the bps preserving defect.  Journal of High Energy Physics, 2019(3):131, 2019.
  * [26] Christof Adam, K Oles, T Romanczukiewicz, and A Wereszczynski.  Kink-antikink collisions in a weakly interacting $\phi^{4}$ model.  Physical Review E, 102(6):062214, 2020.
  * [27] Robert Arthur, Patrick Dorey, and Robert Parini.  Breaking integrability at the boundary: the sine-gordon model with robin boundary conditions.  Journal of physics A: mathematical and theoretical, 49(16):165205, 2016.
  * [28] Patrick Dorey, Aliaksei Halavanau, James Mercer, Tomasz Romanczukiewicz, and Yasha Shnir.  Boundary scattering in the $\phi^{4}$ model.  Journal of High Energy Physics, 2017(5):107, 2017.
  * [29] Fred C Lima, Fabiano C Simas, KZ Nobrega, and Adalto R Gomes.  Boundary scattering in the $\phi^{6}$ model.  Journal of High Energy Physics, 2019(10):1–30, 2019.
  * [30] Patrick Dorey, Anastasia Gorina, Ilya Perapechka, Tomasz Romańczukiewicz, and Yakov Shnir.  Resonance structures in kink-antikink collisions in a deformed sine-gordon model.  Journal of High Energy Physics, 2021(9):1–38, 2021.
  * [31] NS Manton, Katarzyna Oleś, Tomasz Romańczukiewicz, and Andrzej Wereszczyński.  Collective coordinate model of kink-antikink collisions in $\phi^{4}$ theory.  Physical Review Letters, 127(7):071601, 2021.
  * [32] NS Manton, K Oleś, T Romańczukiewicz, and A Wereszczyński.  Kink moduli spaces: Collective coordinates reconsidered.  Physical Review D, 103(2):025024, 2021.
  * [33] C Adam, NS Manton, K Oles, T Romanczukiewicz, and A Wereszczynski.  Relativistic moduli space for kink collisions.  Physical Review D, 105(6):065012, 2022.
  * [34] A Halavanau, T Romanczukiewicz, and Ya Shnir.  Resonance structures in coupled two-component $\phi^{4}$ model.  Physical Review D, 86(8):085027, 2012.
  * [35] A Alonso-Izquierdo.  Reflection, transmutation, annihilation, and resonance in two-component kink collisions.  Physical Review D, 97(4):045016, 2018.
  * [36] Alberto Alonso-Izquierdo.  Non-topological kink scattering in a two-component scalar field theory model.  Communications in Nonlinear Science and Numerical Simulation, 85:105251, 2020.
  * [37] A Alonso-Izquierdo, MA González León, J Martin Vaquero, and M de la Torre Mayado.  Kink scattering in a generalized wess-zumino model.  Communications in Nonlinear Science and Numerical Simulation, 103:106011, 2021.
  * [38] R Jackiw and Paolo Rossi.  Zero modes of the vortex-fermion system.  Nuclear Physics B, 190(4):681–691, 1981.
  * [39] Jeffrey Goldstone and Frank Wilczek.  Fractional quantum numbers on solitons.  Physical Review Letters, 47(14):986, 1981.
  * [40] Tanmay Vachaspati.  Kinks and domain walls: An introduction to classical and quantum solitons.  Cambridge University Press, 2006.
  * [41] Alexander Vilenkin and E Paul S Shellard.  Cosmic strings and other topological defects.  Cambridge University Press, 2000.
  * [42] Farid Charmchi and Siamak S Gousheh.  Complete spectral analysis of the jackiw-rebbi model, including its zero mode.  Physical Review D, 89(2):025002, 2014.
  * [43] Farid Charmchi and Siamak S Gousheh.  Massive jackiw–rebbi model.  Nuclear Physics B, 883:256–266, 2014.
  * [44] André Amado and Azadeh Mohammadi.  Coupled fermion–kink system in jackiw–rebbi model.  The European Physical Journal C, 77(7):465, 2017.
  * [45] Vladislav Klimashonok, Ilya Perapechka, and Yakov Shnir.  Fermions on kinks revisited.  Physical Review D, 100(10):105003, 2019.
  * [46] Ilya Perapechka, Nobuyuki Sawado, and Ya Shnir.  Soliton solutions of the fermion-skyrmion system in (2+ 1) dimensions.  Journal of High Energy Physics, 2018(10):81, 2018.
  * [47] Ilya Perapechka and Yakov Shnir.  Kinks bounded by fermions.  Physical Review D, 101(2):021701, 2020.
  * [48] I Perapechka and Ya Shnir.  Fermion exchange interaction between magnetic skyrmions.  Physical Review D, 99(12):125001, 2019.
  * [49] Yi-Zen Chu and Tanmay Vachaspati.  Fermions on one or fewer kinks.  Physical Review D, 77(2):025006, 2008.
  * [50] Y Brihaye and T Delsate.  Remarks on bell-shaped lumps: stability and fermionic modes.  Physical Review D, 78(2):025014, 2008.
  * [51] Gary Gibbons, Kei-ichi Maeda, and Yu-ichi Takamizu.  Fermions on colliding branes.  Physics Letters B, 647(1):1–7, 2007.
  * [52] Paul M Saffin and Anders Tranberg.  Particle transfer in braneworld collisions.  Journal of High Energy Physics, 2007(08):072, 2007.
  * [53] João GF Campos and Azadeh Mohammadi.  Fermion transfer in the $\phi^{4}$ model with a half-bps preserving impurity.  Physical Review D, 102(4):045003, 2020.
  * [54] João GF Campos and Azadeh Mohammadi.  Fermions on wobbling kinks: normal versus quasinormal modes.  Journal of High Energy Physics, 2021(9):103, 2021.
  * [55] M Shifman, A Vainshtein, and M Voloshin.  Anomaly and quantum corrections to solitons in two-dimensional theories with minimal supersymmetry.  Physical Review D, 59(4):045016, 1999.
  * [56] K Shizuya.  Superfield formulation of central charge anomalies in two-dimensional supersymmetric theories with solitons.  Physical Review D, 69(6):065021, 2004.