11institutetext:  European Space Agency, Noordwijk, 2201 AZ, The Netherlands
11email: {marcus.maertens<EMAIL_ADDRESS>
# Symbolic Regression for Space Applications: Differentiable Cartesian Genetic
Programming Powered by Multi-objective Memetic Algorithms

Marcus Märtens  0000-0003-1950-7111 Dario Izzo  0000-0002-9846-8423

###### Abstract

Interpretable regression models are important for many application domains, as
they allow experts to understand relations between variables from sparse data.
Symbolic regression addresses this issue by searching the space of all
possible free form equations that can be constructed from elementary algebraic
functions. While explicit mathematical functions can be rediscovered this way,
the determination of unknown numerical constants during search has been an
often neglected issue. We propose a new multi-objective memetic algorithm that
exploits a differentiable Cartesian Genetic Programming encoding to learn
constants during evolutionary loops. We show that this approach is competitive
or outperforms machine learned black box regression models or hand-engineered
fits for two applications from space: the Mars express thermal power
estimation and the determination of the age of stars by gyrochronology.

###### Keywords:

Cartesian Genetic Programming Symbolic Regression Memetic Algorithms.

##  1 Introduction

Regression models, frequently based on machine learning concepts [7] like
support vector machines [4], decision trees [14], random forests [3] or neural
networks [13] are common methods to predict or estimate properties of
interest. However, a difficult balance needs to be found when regression
arises as part of a real-world engineering problem: while complex models may
provide a high accuracy, they can be opaque. For example, while certain neural
network models are powerful enough to approximate any function [8], their
nonlinear and heavily parameterized layered structure is prohibitive to
interpretation. For many space-related problems, a lack of interpretability is
a show stopper, as predictability and robustness of the systems are of
critical importance.

Consequently, simpler models like linear regressors [24] are favored in such
cases, as they allow to study the behavior of the model analytically. The
price of assuming a linear model is often a suboptimal accuracy if non-linear
relationships are dominant in the underlying data. A human expert in the loop
might instead impose certain types of functions (e.g., logarithmic,
polynomial, etc.) to improve the fit, but this pragmatic approach introduces
human bias and scalalibity problems.

Symbolic regression approaches aim to substitute this expert by an impartial
system that will discover the best function by evolutionary recombination of
basic mathematical operators, like addition, multiplication, trigonometrical
functions, exponentials and similar. This allows various approaches [22, 20,
15] the rediscovery of explicit mathematical functions from sparse samples. In
contrast, work that applies the same techniques to practical engineering
problems is, putting a few noteworthy industrial exceptions [23] aside,
seemingly rare, raising the question why this is the case?

A potential explanation could be that numerical constants are essential to
describe any reasonably noisy and complex application. Constants are, however,
rarely found in the expressions deployed to test the effectiveness of symbolic
regression techniques. While those techniques excel at recombining the basic
operators, they are poorly equipped to invent constants, either relying on
providing the correct constants already as a fixed (sometimes mutatable) input
or producing them by idiosyncratic combinations of operators.111For example,
the constant $3$ might be constructed from an arbitrary input variable $x$ by
evolving the expression $\frac{x}{x}+\frac{x}{x}+\frac{x}{x}$.

In this report, we introduce a new multi-objective memetic algorithm based on
differentiable Cartesian Genetic Programming (dCGP) [11, 10] that is capable
of finding interpretable mathematical expressions for two different real world
problems coming from the space domain without the need for specific expert
knowledge. We show that by exploiting the differentiability of a dCGP
expression, the best fitting constants can be learned directly. The evolved
expressions are shown to be either competitive or outperforming state of the
art reference models for the tasks under consideration, including machine
learned regression models. Due to its multi-objective nature, a decision maker
is able to balance interpretability (i.e. expression complexity) with model
accuracy.

Our report is structured as follows: Section 2 provides some context of
previous symbolic regression works and relates them to ours. Section 3 briefly
explains dCGP and introduces our memetic and multi-objective algorithm.
Section 4 and Section 5 describe two case studies, one from the Mars express
mission and one from star dating, to which our algorithm is applied and
evaluated. We conclude our work in Section 6.

##  2 Related Work

Early works on symbolic regression were inspired by the Genetic Programming
approach of Koza [12] and deployed encodings based on arithmetic trees or
grammars to evolve mathematical expressions. The encoding of our chromosomes
is based on an extension of Cartesian Genetic Programming (CGP), which was
originally developed by Miller [16].A comprehensive overview about recent
genetic symbolic regression techniques and their applications is given by Wang
et. al [23]. Most influential is the seminal work of Michael Schmidt and Hod
Lipson [19] who extracted conservation laws from observing physical systems
like the double pendulum. Their evolutionary technique is available in the
commercial “Eureqa” software. Eureqa uses an evolutionary technique to
determine the numerical value of constants.

There also exist approaches (the deterministic FFX algorithm [15]) that do not
rely on genetic programming. Most recently, Udrescu and Tegmark [22]
introduced the “AI Feynman algorithm” which relies on a series of regressions
techniques, brute-forcing and neural networks instead of evolutionary
approaches. However, AI Feynman requires constants to be provided as input a
priori or must otherwise construct them by operator recombination.

##  3 Symbolic Regression with dCGP

###  3.1 dCPG

Differentiable Cartesian Genetic Programming (dCGP) is a recent development in
the field of Genetic Programming (GP) adding to genetic programs their any-
order differential representation to be used in learning tasks [10]. The use
of low-order differentials, i.e. gradients, to learn parameters in Genetic
Programs is rare in previous works and can be found, for example, in [21] and
[5]. In past research, attempts to use genetic operators such as crossover and
mutation [9] as well as meta-heuristics such as simulated annealing [6] were
made to adapt such parameters. It is clear how access to the model’s
differential information, when available, may bring substantial gains if used
during learning as proved, for example, by the enormous success of stochastic
gradient descent in machine learning and in particular in learning deep neural
network parameters. Key to the successful use of the differential information
is its availability at a low computational cost, a feature that for first
order derivatives, i.e. the gradient, is guaranteed by the backward automated
differentiation technique, while for higher order derivatives is more
difficult to obtain. Thanks to the built-in efficient implementation of the
algebra of Taylor truncated polynomials, evaluating a dCGP program allows
seamless access to gradients and Hessians (as well as higher order
differential information). In the following, we describe how our new Multi
Objective Memetic Evolutionary Strategy (MOMES) leverages the loss gradient
and its Hessian to evolve the dCGP program and learn the ephemeral constants
simultaneously.

###  3.2 Multi Objective Memetic Evolutionary Strategy (MOMES)

Let us denote with $\boldsymbol{\eta}=(\boldsymbol{\eta}_{\\_}u,\mathbf{c})$ a
generic chromosome encoding the symbolic expression
$\hat{y}=f_{\\_}{\boldsymbol{\eta}_{\\_}u}(\mathbf{x},\mathbf{c})$. We
recognize in $\boldsymbol{\eta}$ two distinct parts:
$\boldsymbol{\eta}_{\\_}u$ is the integer part of the chromosome making use of
the classical Cartesian Genetic Programming [16] encoding. Thus,
$\boldsymbol{\eta}_{\\_}u$ represents the mathematical form of $\mathbf{f}$,
while $\mathbf{c}$ represents the actual values of the ephemeral constants in
a continuous direct encoding. We may then introduce the gradient and Hessian
as:

$\begin{array}[]{l}\nabla\hat{y}=[\frac{\partial\hat{y}}{\partial
c_{\\_}1},...,\frac{\partial\hat{y}}{\partial c_{\\_}n}]\hskip
42.67912pt\nabla^{2}\hat{y}=\begin{bmatrix}\frac{\partial\hat{y}^{2}}{\partial
c_{\\_}1\partial c_{\\_}1}&\dots&\frac{\partial\hat{y}^{2}}{\partial
c_{\\_}1\partial c_{\\_}n}\\\ \vdots&\ddots&\\\
\frac{\partial\hat{y}^{2}}{\partial c_{\\_}n\partial
c_{\\_}1}&&\frac{\partial\hat{y}^{2}}{\partial c_{\\_}n\partial
c_{\\_}n}\end{bmatrix}\end{array}$

Assuming to work on labelled data denoted by
$\mathcal{D}=(\mathbf{x}_{\\_}i,y_{\\_}i)$, $i=1..N$ we may thus compute, for
any chromosome $\boldsymbol{\eta}$, the mean square error (MSE) loss:

$\ell=\sum_{\\_}i(y_{\\_}i-\hat{y}_{\\_}i)^{2}$

its full Hessian ${\mathbf{H}}$ and its full gradient ${\mathbf{G}}$. We also
define the complexity $\mu$ of the chromosome $\boldsymbol{\eta}$ as the
number of active nodes in its CGP. Let us now introduce the concept of active
constants $\tilde{\mathbf{c}}$ as the selection from $\mathbf{c}$ of all
components with nonzero entries in the gradient ${\mathbf{G}}$. Zero entries
correspond, most of the time, to ephemeral constants that are not expressed by
the corresponding CGP and are thus inactive. Similarly, we call
$\tilde{\mathbf{H}}$ and $\tilde{\mathbf{G}}$, respectively, the loss active
Hessian and the loss active gradient when considering only the active
constants. Following these notations and definitions, we may now describe the
MOMES algorithm.

MOMES is an evolutionary strategy where all members in a population of $NP$
individuals, each represented by a chromosome $\boldsymbol{\eta}_{\\_}i$,
$i=1..NP$, undergo a mutation step at each generation, acting on all the genes
(active and inactive) in $\boldsymbol{\eta}_{\\_}u$. The mutated chromosome is
then subject to a lifelong learning process consisting of one single step of
Newton’s method and acting on the (active) continuous part
$\tilde{\mathbf{c}}$ of the chromosome:

$\tilde{\mathbf{c}}_{\\_}{new}=\tilde{\mathbf{c}}_{\\_}{old}-\tilde{\mathbf{H}}^{-1}\tilde{\mathbf{G}}$

The new individual fitness is then evaluated and, to preserve diversity, the
chromosome is added to the $NP$ parents only if its fitness is not already
present in the candidate pool. Non-dominated sorting over the objectives
$(\ell,\mu)$ is then applied to the candidate pool selecting $NP$ individuals
to be carried over to the new generation. The exact implementation details of
MOMES are documented in the open source project dCGP [10].

The design of MOMES required the use of a few ideas that were the results of
experiments and iterations that are worth discussing briefly at this point.
Firstly, the mutation of the CGP expression encoded in
$\boldsymbol{\eta}_{\\_}u$ acts randomly on all genes, also the ones that are
not expressed by the resulting expression. This type of mutation, already
noted as beneficial when used in the standard evolutionary strategy technique
used for CGP [16], has an added benefit here as it allows for a mutated
chromosome to express the same formula as the parent. In these cases, the
subsequent lifelong learning ensures the possibility to apply further
iterations of the Newton’s method to the same expression, and the mechanism
will thus ensure that, eventually, a full local descent and not only one
single iteration is made over the most promising individuals. Additionally,
this simple trick acts as a necessary diversity preservation mechanism as
otherwise the very same expression (possibly having a different complexity as
counted by the number of active nodes in the CGP) would flood the non-
dominated front and lead to a degraded convergence. Secondly, the choice of
using a single step of a second order method for the memetic part of MOMES was
found to be most effective when coupled with MSE loss. It has been noted (see
[11]) that, whenever constants appear linearly in the expressed CGP, one
single step of a second order method is enough to get to the optimal choice if
the loss is expressed as MSE. This mechanism adds a small, but beneficial bias
towards simple expressions during evolution, since these are evaluated with
optimal constants values at each generation, while more complex appearances of
constants would require more generations to reach their optimal values.
Lastly, the introduction of active constants, Hessians and gradients allows
the Hessian matrix to admit an inverse and thus enable an effective lifelong
learning.

##  4 Case 1: Mars Express Thermal Data

###  4.1 Background and Dataset

Mars Express (MEX) is a spacecraft of the European Space Agency (ESA) that has
been monitoring the “red planet” continuously since 2004. To remain operable
in the challenging environment of space, the probe is equipped with a thermal
regulation system to which a certain power has to be allocated. An accurate
power budget is essential for the mission, as underestimates might lead to
malfunctions and overestimates will waste energy that should have been
allocated to scientific tasks instead. The MEX data for this work [18] is
split in a training set MEX1 and a test set MEX2. MEX1 consists of 5795
completed orbits in the time frame from 01.01.2015 to 18.08.2019 and MEX2 to
1507 orbits from 18.08.2019 to 31.10.2020. Measurements of potential relevance
for the prediction of the thermal power budget have been aggregated for each
orbit into several thermal contributors. In particular, the task is to find
the thermal power demand $P_{\\_}{th}$ as a function

$P_{\\_}{th}(LVAH,SH,D_{\\_}{ecl},TX,FO,NS).$ (1)

A reference model based on linear regression

$P_{\\_}{th}=c\cdot LVAH+d\cdot SH+e\cdot D_{\\_}{ecl}+f\cdot TX+g\cdot
FO+h\cdot NS+i$ (2)

was fitted to the MEX1 data. Table 1 gives a short description of each thermal
contributor and the value of the fitted coefficients. Figure 1 shows the MEX
data in relation to orbital conditions and actual power consumption. The
evaluation metric for all MEX data is the Root Mean Square Error (RMSE).

thermal contri- butor |  description |  related coefficient |  dCGP variable
---|---|---|---
$LVAH$ | Launch Vehicle Adapter heating | $c=-7.666\cdot 10^{1}$ | $x_{\\_}0$
$SH$ | solar heating | $d=-1.764\cdot 10^{-1}$ | $x_{\\_}1$
$D_{\\_}{ecl}$ | heating due to eclipses | $e=-3.387\cdot 10^{-3}$ | $x_{\\_}2$
$TX$ | transmitter activities | $f=-6.898\cdot 10^{0}$ | $x_{\\_}3$
$FO$ | flag indicating power off | $g=+1.107\cdot 10^{1}$ | $x_{\\_}4$
$NS$ | guidance flag | $h=+6.820\cdot 10^{0}$ | $x_{\\_}5$
- | (baseline) | $i=+2.267\cdot 10^{2}$ |
Table 1: The reference model THP for MEX data is a linear regression model
with the highlighted coefficients (compare Equation 2). Last column shows
corresponding variables for dCGP expressions. Figure 1: MEX data, including
orbital and operative parameters. Goal is to predict the average power
consumption from thermal contributors, compare Equation 1.

###  4.2 Results

For this experiment, MEX1 data is used for training and MEX2 data for
additional testing. Our features are the thermal contributors of the reference
THP model (see Table 1), with two differences: $SH$ and $D_{\\_}{ecl}$ are
standardized by subtracting the mean and dividing by the standard deviation of
the training set. The determination of optimal hyperparameters for dCGP and
MOMES is not within the scope of this work, as it would require us to
thoroughly explain, characterize and statistically analyze each of them. In
practice, evolution converges fast enough to perform grid searches over some
reasonable parameter guesses that can be further calibrated over a few
iterations of the optimization pipeline. Table 2 reports the hyperparameters
that we found to be effective enough to make our point with this experiment.

Figure 2: Top: Non-dominated front found by MOMES, 18 out of 40 individuals, MEX1 data. Blue line is the reference model THP. The expressions M1, M2 and M3 are highlighted for further analysis. Bottom: Predictions of THP and dCGP expression M1 in direct comparison with target power consumption. Shaded region is test data MEX2. parameter |  value
---|---
basic operators | $+,-,\times,\div,\log$
rows | 2
columns | 20
levels back | 20
maximum number of constants | 5
maximum mutations | 4
generations | 50000
population size | 40
Table 2: Hyperparameters deployed for dCGP and MOMES for the MEX1 dataset.

As common with evolutionary techniques, we multi-start MOMES 200 times to
account for its stochasticity. Out of those 200 results, we select the one
that has the strongest extreme point (i.e. lowest RMSE) and show the complete
Pareto front in Figure 2. For further comparison with the THP reference model,
we select three different formulas: M1 is at the extreme point, M2 is the
expression with minimum complexity that is still above the reference and M3 is
a low complexity expression. Table 3 shows the corresponding RMSE and the
average over- and underestimations.

MEX1

Model |  RMSE |  complx. |  avg. over- estimate |  avg. under- estimate
---|---|---|---|---
THP | 8.672 | - | 6.665 | 6.774
M1 | 8.478 | 61 | 6.657 | 6.442
M2 | 8.647 | 34 | 6.826 | 6.619
M3 | 9.223 | 22 | 7.184 | 7.232

MEX2

Model |  RMSE |  complx. |  avg. over- estimate |  avg. under- estimate
---|---|---|---|---
THP | 12.454 | - | 10.684 | 4.549
M1 | 11.390 | 61 | 9.753 | 4.461
M2 | 12.828 | 34 | 10.924 | 4.998
M3 | 14.991 | 22 | 12.928 | 6.745

Table 3: RMSE for 3 selected expression from Figure 2, complexities and values
for average over- and underestimation of thermal power consumption. THP is the
linear regression reference model.

We observe that the RMSE for MEX2 is considerably higher than for MEX1 for
each of the models, with Expression M1 generalizing slightly better to the
unseen conditions of MEX2 than THP. We show the actual predictions of both,
THP and Expression M1 in Figure 2, which highlights the differences between
the models. From a mathematical point of view, M1 includes several products of
variables, resulting in a degree 2 polynomial. Despite having access to the
division and logarithm operator, evolution did not favor them in the
population. However, other non-dominated fronts from the 200 runs (not shown)
included sometimes division, constructing rational functions of similar RMSE.
M3 is a linear model without the thermal contributors $x_{\\_}2$ and
$x_{\\_}3$ that is worse than THP but could provide an interesting alternative
if certain measurements would not be readily available. While some of the
numerical constants are similar across different expression, it is clear that
the memetic algorithm is optimizing them for each expression independently to
further reduce the error. Furthermore, the maximum of 5 constants is not
(always) expressed, highlighting that MOMES is selective and parsimonious. In
comparison, the simple linear regression model THP relies already on 7
explicit constants.

##  5 Case 2: Star Dating Gyrochronology

###  5.1 Background and Dataset

This dataset was released by Moya et al. [17] as a benchmark for AI research
on star dating. Many physical models for astronomical problems (e.g. the
search for habitable exoplanets) depend on the age of a star, which is however
not directly measurable. Consequently, the age of a star needs to be inferred
from observable features. The public benchmark describes 1464 stars whose
accurate ages have been determined by asteroseismology and association to star
clusters, providing for each star the stellar features described in Table 4.

In this context an approach termed “gyrochronology” is of particular interest,
which hypothesizes a functional relationship between star age, stellar mass
and rotational period. Some empirical relations (assuming a fixed functional
form) have been proposed [1] but the accuracy of gyrochronology remains
controversial, with some works suggesting that linear relationships might not
be sufficient at all or just applicable for stars of a certain age group [2].
Consequently, the authors of the dataset propose to study star dating directly
as regression problem using supervised machine learning models.

Feature |  Description |  dCGP variable
---|---|---
$M$ | Stellar Mass | $x_{\\_}0$
$R$ | Stellar Radius | $x_{\\_}1$
$T_{\\_}{eff}$ | Stellar effective temperature | $x_{\\_}2$
$L$ | Luminosity | $x_{\\_}3$
$[Fe/H]$ | Metallicity | $x_{\\_}4$
$\log g$ | Logarithm of surface gravity | $x_{\\_}5$
$P_{\\_}{rot}$ | Stellar rotational period | $x_{\\_}6$
Table 4: Features available for every star in the dataset.

Several reference models including linear regression, decision tree regressor,
random forest regressor, support vector regressor, Gaussian process, kNN,
neural networks and ensembles are provided within the benchmark. Additionally,
four different splits of the data for training and test are available:

  * •

A a random 80/20 split

  * •

B1 stars of age $[0.4,4.2]$ Gyr are used for training, ages $[4.2,13.8]$ for
testing

  * •

B2 stars dated by cluster belonging used for training, remainder for testing

  * •

C independent control with 32 new stars for testing, including our own sun

The evaluation metric for all gyrochronology models is the Mean Absolute Error
(MAE). Additionally, the dataset provides error bounds on the age of the stars
allowing to define a precision metric by the fraction of star age predictions
that fall within the corresponding confidence interval.

###  5.2 Results

For this experiment, we deployed MOMES on all training sets of the given
benchmarks and re-evaluated the found expressions on the test sets. We
selected the 6 features as described in Table 4 and scaled them by dividing
each of them by the standard deviation of the corresponding training set.
Similar to Section 4, a grid-search was used to find suitable hyperparameters
for dCGP and MOMES, reported in Table 5.

parameter |  value
---|---
basic operators | $+,-,\times,\div,\log,\sin$
rows | 2
columns | 16
levels back | 16
number of constants | 5
maximum mutations | 6
generations | 500000
population size | 40
Table 5: Hyperparameters deployed for dCGP and MOMES for star dating.

A multi-start setup of 400 independent runs (100 for each benchmark) has been
deployed to our servers and completed in about 12 hours. We selected the run
that resulted in the lowest MAE (extreme point) on the training set for
further analysis. Table 6 shows the MAE and the precision of the lowest MAE
expression found in the non-dominated front of MOMES in comparison with the
reference models. To exemplify the diversity and complexity trade-offs in the
found expressions, we show the approximated Pareto fronts in Figure 3 for
Benchmark A and B1 (for reasons of space).

Interpreting the results, the dCGP based models are competitive with some of
the machine learning models for benchmark A, but have a higher error than
neural networks, Gaussian processes, kNN and random forests. Benchmark B1
represents a domain gap study, which is (without the deployment of additional
techniques or assumptions) hard for most machine learning models to bridge.
Since dCGP looks to capture the physical relations between the different
features in algebraic expressions, it leads itself to better generalization.

Investigating the Pareto front of dCGP on the test-set for B2, the simple
expression $x_{\\_}6-\sin(x_{\\_}1-x_{\\_}4)+6.987$ could be shown to
outperforms all other machine learning models by some margin. On a closer look
however, it turned out that the correlation between the training error and
test error in this benchmark was rather low when looking at the entirety of
our MOMES runs. Although this particular result for B2 has to be taken with a
grain of salt, it nevertheless demonstrates the existence of surprisingly
simple expressions that provide an accurate explanation of the data. Thanks to
the interpretability of dCGP expressions, these type of insights may uncover
potential issues related to overfitting of machine learning models or issues
related to the data split that would go unnoticed otherwise. Lastly, on
benchmark C the best dCGP expressions achieves a low MAE, only surpassed by
the Gaussian process model.

Figure 3: Non-dominated front found by MOMES for Benchmark A (top) and B1
(bottom). Shown is the non-dominated front after the complete population was
re-evaluated on the test-set of its corresponding benchmark. Blue lines
represent reference machine learned models as described by Moya et al. [17]

Precision MAE

|  dcgp |  nnet |  lr |  dtr |  rf |  svm |  bayes |  knn |  gp
---|---|---|---|---|---|---|---|---|---
A | 0.62 | 0.41 | 0.95 | 0.70 | 0.57 | 0.68 | 0.95 | 0.53 | 0.45
B1 | 1.41 | 1.48 | 1.92 | 1.62 | 1.71 | 1.72 | 1.94 | 1.61 | 1.67
B2 | 1.39 | 1.85 | 1.56 | 1.90 | 1.85 | 1.51 | 1.48 | 1.84 | 1.62
C | 0.97 | 1.10 | 1.30 | 2.11 | 1.01 | 1.28 | 1.30 | 1.73 | 0.84
|  dcgp |  nnet |  lr |  dtr |  rf |  svm |  bayes |  knn |  gp
A | 50.39 | 73.23 | 35.43 | 63.78 | 60.63 | 46.46 | 32.28 | 77.95 | 62.20
B1 | 51.97 | 59.84 | 29.92 | 51.18 | 48.03 | 37.79 | 28.35 | 55.91 | 51.97
B2 | 42.26 | 28.45 | 30.96 | 21.34 | 25.94 | 33.05 | 29.71 | 25.52 | 30.96
C | 15.63 | 34.37 | 21.87 | 9.37 | 21.87 | 12.50 | 21.87 | 15.62 | 40.62

Table 6: MAE and precision on test set for all star dating benchmarks. Each
model was trained for each benchmark separately. The dCGP expressions used are
always the ones that are at the extreme end of the MAE objective of the non-
dominated front created by MOMES.

##  6 Conclusions

Following our experiments on the Mars express and star dating datasets, we
have demonstrated that MOMES is capable of finding algebraic expressions that
explain unseen data and have the capability to generalize well to it. The fact
that the numerical constants for each expression are learned at each step of
evolution guides the memetic algorithm towards expressions that are both
accurate and of low complexity. The expressions with the lowest error notably
contain products of sub-expressions or non-linear operators like the
$\sin$-function and are thus more expressive than linear regressors while
still amenable to mathematical analysis. Consequently, no expert is required
to make potentially biased guesses in finding non-linear expressions. The task
that still remains for the expert is the engineering and scaling of features
(e.g. the thermal contributors in case of MEX) to obtain the best possible
results, as well as the determination of suitable basic operators and
efficient hyperparameters.

As shown in the case of benchmark B2 for the star dating case study, deploying
symbolic regression in general is beneficial as sometimes simple expressions
appear to be more general than complex black box models. Given the practical
benefits of interpretability that comes with algebraic expressions together
with the low requirements for their inference make algorithms like MOMES
valuable methods of knowledge discovery. In particular, applications with
sparse data in extreme environments like space are likely to benefit from it.

#### Acknowledgements

The authors are grateful to Thomas Dreßler, who helped making the MEX data
public and explained it to us. Similarly, the authors want to thank Roberto J.
López-Sastre for providing support and insight into the gyrochronology data.

## References

  * [1] Angus, R., Aigrain, S., Foreman-Mackey, D., McQuillan, A.: Calibrating gyrochronology using Kepler asteroseismic targets. Monthly Notices of the Royal Astronomical Society 450(2), 1787–1798 (04 2015). https://doi.org/10.1093/mnras/stv423
  * [2] Barnes, S.A.: A simple nonlinear model for the rotation of main-sequence cool stars. i. introduction, implications for gyrochronology, and color–period diagrams. The Astrophysical Journal 722(1), 222 (2010). https://doi.org/10.1088/0004-637X/722/1/222
  * [3] Breiman, L.: Random forests. Machine learning 45(1), 5–32 (2001). https://doi.org/10.1023/A:1010933404324
  * [4] Chang, C.C., Lin, C.J.: Libsvm: a library for support vector machines. ACM transactions on intelligent systems and technology (TIST) 2(3), 1–27 (2011). https://doi.org/10.1145/1961189.1961199
  * [5] Emigdio, Z., Trujillo, L., Schütze, O., Legrand, P., et al.: A local search approach to genetic programming for binary classification. In: Proceedings of the 2015 on Genetic and Evolutionary Computation Conference-GECCO’15 (2015). https://doi.org/10.1145/2739480.2754797
  * [6] Esparcia-Alcazar, A.I., Sharman, K.C.: Genetic programming techniques that evolve recurrent neural network architectures for signal processing. In: Neural Networks for Signal Processing [1996] VI. Proceedings of the 1996 IEEE Signal Processing Society Workshop. pp. 139–148. IEEE (1996). https://doi.org/10.1109/NNSP.1996.548344
  * [7] Fernández-Delgado, M., Sirsat, M.S., Cernadas, E., Alawadi, S., Barro, S., Febrero-Bande, M.: An extensive experimental survey of regression methods. Neural Networks 111, 11–34 (2019). https://doi.org/10.1016/j.neunet.2018.12.010
  * [8] Hornik, K., Stinchcombe, M., White, H.: Multilayer feedforward networks are universal approximators. Neural networks 2(5), 359–366 (1989). https://doi.org/10.1016/0893-6080(89)90020-8
  * [9] Howard, L.M., D’Angelo, D.J.: The GA-P: A genetic algorithm and genetic programming hybrid. IEEE Expert 10(3), 11–15 (1995)
  * [10] Izzo, D., Biscani, F.: dcgp: Differentiable cartesian genetic programming made easy. Journal of Open Source Software 5(51), 2290 (2020). https://doi.org/10.21105/joss.02290
  * [11] Izzo, D., Biscani, F., Mereta, A.: Differentiable genetic programming. In: European conference on genetic programming. pp. 35–51. Springer (2017). https://doi.org/10.1007/978-3-319-55696-33
  * [12] Koza, J.R.: Genetic programming as a means for programming computers by natural selection. Statistics and computing 4(2), 87–112 (1994). https://doi.org/10.1007/BF00175355
  * [13] Lathuilière, S., Mesejo, P., Alameda-Pineda, X., Horaud, R.: A comprehensive analysis of deep regression. IEEE transactions on pattern analysis and machine intelligence 42(9), 2065–2081 (2019). https://doi.org/10.1109/TPAMI.2019.2910523
  * [14] Loh, W.Y.: Classification and regression trees. Wiley interdisciplinary reviews: data mining and knowledge discovery 1(1), 14–23 (2011). https://doi.org/10.1002/widm.8
  * [15] McConaghy, T.: FFX: fast, scalable, deterministic symbolic regression technology. In: Genetic Programming Theory and Practice IX, pp. 235–260. Springer (2011)
  * [16] Miller, J.F.: Cartesian genetic programming. In: Cartesian Genetic Programming, pp. 17–34. Springer (2011). https://doi.org/10.1007/978-3-642-17310-32
  * [17] Moya, A., Recio-Martinez, J., López-Sastre, R.J.: Ai for dating stars: a benchmarking study for gyrochronology. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 1971–1981 (2021). https://doi.org/10.1109/CVPRW53098.2021.00225
  * [18] Märtens, M., Dreßler, T.: Mars express thermal power dataset (Apr 2022). https://doi.org/10.5281/zenodo.6417900
  * [19] Schmidt, M., Lipson, H.: Distilling free-form natural laws from experimental data. science 324(5923), 81–85 (2009). https://doi.org/10.1126/science.1165893
  * [20] Schmidt, M.D., Lipson, H.: Solving iterated functions using genetic programming. In: Proceedings of the 11th Annual Conference Companion on Genetic and Evolutionary Computation Conference: Late Breaking Papers. pp. 2149–2154 (2009). https://doi.org/10.1145/1570256.1570292
  * [21] Topchy, A., Punch, W.F.: Faster genetic programming based on local gradient search of numeric leaf values. In: Proceedings of the 3rd Annual Conference on Genetic and Evolutionary Computation. pp. 155–162. Morgan Kaufmann Publishers Inc. (2001)
  * [22] Udrescu, S.M., Tegmark, M.: AI Feynman: A physics-inspired method for symbolic regression. Science Advances 6(16) (2020). https://doi.org/10.1126/sciadv.aay2631
  * [23] Wang, Y., Wagner, N., Rondinelli, J.M.: Symbolic regression in materials science. MRS Communications 9(3), 793–805 (2019). https://doi.org/10.1557/mrc.2019.85
  * [24] Weisberg, S.: Applied linear regression, vol. 528. John Wiley & Sons (2005)


-title

11institutetext: Institute of Electrophysics & Radiation Technologies,
NAS of Ukraine, Kharkov, Ukraine 22institutetext: NSC Kharkov Institute of
Physics & Technology
NAS of Ukraine, Kharkov, Ukraine

# The method of unitary clothing transformations in the theory of
nucleon–nucleon scattering

I. Dubovyk 11<EMAIL_ADDRESS>A. Shebeko 22<EMAIL_ADDRESS>
###### Abstract

The clothing procedure, put forward in quantum field theory (QFT) by Greenberg
and Schweber, is applied for the description of nucleon–nucleon ($N$–$N$)
scattering. We consider pseudoscalar ($\pi$ and $\eta$), vector ($\rho$ and
$\omega$) and scalar ($\delta$ and $\sigma$) meson fields interacting with
$1/2$ spin ($N$ and $\bar{N}$) fermion ones via the Yukawa–type couplings to
introduce trial interactions between ”bare” particles. The subsequent unitary
clothing transformations (UCTs) are found to express the total Hamiltonian
through new interaction operators that refer to particles with physical
(observable) properties, the so–called clothed particles. In this work, we are
focused upon the Hermitian and energy–independent operators for the clothed
nucleons, being built up in the second order in the coupling constants. The
corresponding analytic expressions in momentum space are compared with the
separate meson contributions to the one–boson–exchange potentials in the meson
theory of nuclear forces. In order to evaluate the $T$ matrix of the $N$–$N$
scattering we have used an equivalence theorem that enables us to operate in
the clothed particle representation (CPR) instead of the bare particle
representation (BPR) with its huge amount of virtual processes. We have
derived the Lippmann–Schwinger(LS)–type equation for the CPR elements of the
$T$–matrix for a given collision energy in the two–nucleon sector of the
Hilbert space $\mathcal{H}$ of hadronic states and elaborated a code for its
numerical solution in momentum space.

##  1 Introductory remarks

We know that there are a number of high precision, boson–exchange models of
the two nucleon force $V_{NN}$, such as Paris Lacom80 , Bonn MachHolElst87 ,
Nijmegen Stocks94 , Argonne WirStocksSchia95 , CD Bonn Mach01  potentials and
a fresh family of covariant one–boson–exchange (OBE) ones GrossStad08 . Note
also successful treatments based on chiral effective field theory
OrdoRayKolck94 ; EpelGloeMeiss00 , for a review see Epel05 .

In this talk, we would like to draw attention to the first application of
unitary clothing transformations SheShi01 ; KorCanShe07  in describing the
nucleon-nucleon ($N$–$N$) scattering. Recall that such transformations $W$,
being aimed at the inclusion of the so–called cloud or persistent effects,
make it possible the transition from the bare–particle representation (BPR) to
the clothed–particle representation (CPR) in the Hilbert space $\mathcal{H}$
of meson–nucleon states. In this way, a large amount of virtual processes
induced with the meson absorption/emission, the $N\overline{N}$–pair
annihilation/production and other cloud effects can be accumulated in the
creation (destruction) operators $\alpha_{c}$ for the ”clothed” (physical)
mesons and nucleons. Such a bootstrap reflects the most significant
distinction between the concepts of clothed and bare particles.

In the course of the clothing procedure all the generators of the Poincaré
group get one and the same sparse structure on $\mathcal{H}$ SheShi01 . Here
we will focus upon one of them, viz., the total Hamiltonian

$H=H_{F}(\alpha)+H_{I}(\alpha)\equiv H(\alpha)$ (1)

with

$H_{I}(\alpha)=V(\alpha)+\mbox{mass and vertex counterterms},$ (2)

where free part $H_{F}(\alpha)$ and interaction $V(\alpha)$ depend on creation
(destruction) operators $\alpha^{{\dagger}}(\alpha)$ in the BPR , i.e.,
referred to bare particles with physical masses KorCanShe07 , where they have
been introduced via the mass–changing Bogoliubov–type UTs. To be more
definite, let us consider fermions (nucleons and antinucleons) and bosons
($\pi$–, $\eta$–, $\rho$–, $\omega$–mesons, etc.) interacting via the Yukawa-
type couplings for scalar (s), pseudoscalar (ps) and vector (v) mesons (see,
e.g., MachHolElst87 ). Then, using a trick prompted by the derivation of Eq.
(7.5.22) in WeinbergBook1995  to eliminate in a proper way the vector–field
component $\varphi_{\rm{v}}^{0}$ , we have $V(\alpha)=V_{s}+V_{ps}+V_{\rm{v}}$
with

$V_{s}=g_{s}\int
d\vec{x}\,\bar{\psi}(\vec{x})\psi(\vec{x})\varphi_{s}(\vec{x})$ (3)
$V_{ps}=ig_{ps}\int
d\vec{x}\,\bar{\psi}(\vec{x})\gamma_{5}\psi(\vec{x})\varphi_{ps}(\vec{x})$ (4)
$V_{\rm{v}}=\int
d\vec{x}\,\left\\{g_{\rm{v}}\bar{\psi}(\vec{x})\gamma_{\mu}\psi(\vec{x})\varphi_{\rm{v}}^{\mu}(\vec{x})\phantom{\frac{f}{4}}\right.\\\
\left.+\frac{f_{\rm{v}}}{4m}\bar{\psi}(\vec{x})\sigma_{\mu\nu}\psi(\vec{x})\varphi_{\rm{v}}^{\mu\nu}(\vec{x})\right\\}\\\
+\int
d\vec{x}\,\left\\{\frac{g_{\rm{v}}^{2}}{2m_{\rm{v}}^{2}}\bar{\psi}(\vec{x})\gamma_{0}\psi(\vec{x})\bar{\psi}(\vec{x})\gamma_{0}\psi(\vec{x})\right.\\\
\left.+\frac{f_{\rm{v}}^{2}}{4m^{2}}\bar{\psi}(\vec{x})\sigma_{0i}\psi(\vec{x})\bar{\psi}(\vec{x})\sigma_{0i}\psi(\vec{x})\right\\},$
(5)

with the boson fields $\varphi_{b}$ and the fermion field $\psi$, where
$\varphi_{\rm{v}}^{\mu\nu}(\vec{x})=\partial^{\mu}\varphi_{\rm{v}}^{\nu}(\vec{x})-\partial^{\nu}\varphi_{\rm{v}}^{\mu}(\vec{x})$
the tensor of the vector field included. The mass (vertex) counterterms are
given by Eqs. (32)–(33) of Ref. KorCanShe07  (the difference $V_{0}(\alpha)$
\- $V(\alpha)$ where a primary interaction $V_{0}(\alpha)$ is derived from
$V(\alpha)$ replacing the ”physical” coupling constants by ”bare” ones).

The corresponding set $\alpha$ involves operators $a^{{\dagger}}(a)$ for the
bosons, $b^{{\dagger}}(b)$ for the nucleons and $d^{{\dagger}}(d)$ for the
antinucleons. In their terms, e.g., we have the free pion and fermion parts

$H_{F}(\alpha)=\int
d\vec{k}\,\omega_{\vec{k}}a^{\dagger}(\vec{k})a(\vec{k})\\\ +\int
d\vec{p}\,E_{\vec{p}}\sum\limits_{\mu}\left[b^{\dagger}(\vec{p},\mu)b(\vec{p},\mu)+d^{\dagger}(\vec{p},\mu)d(\vec{p},\mu)\right]$
(6)

and the primary trilinear interaction

$V(\alpha)\sim
ab^{{\dagger}}b+ab^{{\dagger}}d^{{\dagger}}+adb+add^{{\dagger}}+H.c.$ (7)

with the three-legs vertices. Here
$\omega_{\vec{k}}=\sqrt{m^{2}_{b}+\vec{k}^{2}}\,\,\,$($E_{\vec{p}}=\sqrt{m^{2}+\vec{p}^{2}}$)
the pion (nucleon) energy with physical mass $m_{b}(m)$, $\mu$ the fermion
polarization index.

We have tried to draw parallels with that field–theoretic background which has
been employed in boson–exchange models. First of all, we imply the approach by
the Bonn group MachHolElst87 ; Mach01 , where, following the idea by Schütte
Schutte , the authors started from the total Hamiltonian (in our notations),

$H=H_{F}(\alpha)+V(\alpha)$ (8)

with the boson-nucleon interaction

$V(\alpha)\sim ab^{{\dagger}}b+H.c.$ (9)

##  2 Analytic expressions for the quasipotentials in momentum space

As shown in SheShi01 , after eliminating the so-called bad terms111By
definition, they prevent the bare vacuum $\Omega_{0}$
($a|\Omega_{0}\rangle=b|\Omega_{0}\rangle=\ldots=0$) and the bare one–particle
states $|1bare\rangle\equiv a^{{\dagger}}|\Omega_{0}\rangle$
($b^{{\dagger}}|\Omega_{0}\rangle,\ldots$) to be $H$ eigenstates. from
$V(\alpha)$ the primary Hamiltonian $H(\alpha)$ can be represented in the
form,

$H(\alpha)=K_{F}(\alpha_{c})+K_{I}(\alpha_{c})\equiv K(\alpha_{c})$ (10)

The free part of the new decomposition is determined by

$K_{F}(\alpha_{c})=\int
d\vec{k}\,\omega_{\vec{k}}a_{c}^{\dagger}(\vec{k})a_{c}(\vec{k})\\\ +\int
d\vec{p}\,E_{\vec{p}}\sum\limits_{\mu}\left[b_{c}^{\dagger}(\vec{p},\mu)b_{c}(\vec{p},\mu)+d_{c}^{\dagger}(\vec{p},\mu)d_{c}(\vec{p},\mu)\right]$
(11)

while $K_{I}$ contains only interactions responsible for physical processes,
these quasipotentials between the clothed particles, e.g.,

$K_{I}^{(2)}(\alpha_{c})=K(NN\to NN)+K(\bar{N}\bar{N}\to\bar{N}\bar{N})\\\
+K(N\bar{N}\to N\bar{N})+K(bN\to bN)+K(b\bar{N}\to b\bar{N})\\\ +K(bb\to
N\bar{N})+K(N\bar{N}\to bb)$ (12)

In accordance with the clothing procedure developed in SheShi01  they obey the
following requirements:

i) The physical vacuum (the $H$ lowest eigenstate) must coincide with a new
no–particle state $\Omega$, i.e., the state that obeys the equations

$a_{c}(\vec{k})\left|\Omega\right\rangle=b_{c}(\vec{p},\mu)\left|\Omega\right\rangle=d_{c}(\vec{p},\mu)\left|\Omega\right\rangle=0,$
$\forall\,\,{\vec{k},\,\vec{p},\,\mu}\,\,\,\,\,\left(\left\langle\Omega|\Omega\right\rangle=1\right).$
(13)

ii) New one-clothed-particle states $|\vec{k}\rangle_{c}\equiv
a_{c}^{{\dagger}}(\vec{k})\Omega$ etc. are the eigenvectors both of $K_{F}$
and $K$,

$K(\alpha_{c})|\vec{k}\rangle_{c}=K_{F}(\alpha_{c})|\vec{k}\rangle_{c}=\omega_{k}|\vec{k}\rangle_{c}$
(14) $K_{I}(\alpha_{c})|\vec{k}\rangle_{c}=0$ (15)

iii) The spectrum of indices that enumerate the new operators must be the same
as that for the bare ones .

iv) The new operators $\alpha_{c}$ satisfy the same commutation rules as do
their bare counterparts $\alpha$, since the both sets are connected to each
other via the similarity transformation

$\alpha_{c}=W^{{\dagger}}\alpha W,$ (16)

with a unitary operator $W$ to be obtained as in SheShi01 .

It is important to realize that operator $K(\alpha_{c})$ is the same
Hamiltonian $H(\alpha)$. Accordingly [10,11] the $N$–$N$ interaction operator
in the CPR has the following structure:

$K(NN\rightarrow NN)=\sum\limits_{b}K_{b}(NN\rightarrow NN),$
$K_{b}(NN\rightarrow
NN)=\int\sum\limits_{\mu}d\vec{p}^{\prime}_{1}\,d\vec{p}^{\prime}_{2}\,d\vec{p}_{1}\,d\vec{p}_{2}\\\
\times
V_{b}(1^{\prime},2^{\prime};1,2)b_{c}^{\dagger}(1^{\prime})b_{c}^{\dagger}(2^{\prime})b_{c}(1)b_{c}(2),$
(17)

where the symbol $\sum\limits_{\mu}$ denotes the summation over nucleon spin
projections, $1=\\{\vec{p}_{1},\mu_{1}\\}$, etc.

For our evaluations of the c–number matrices $V_{b}$ we have employed some
experience from Refs. SheShi01 ; KorCanShe07  to get in the second order in
the coupling constants

$V_{b}(1^{\prime},2^{\prime};1,2)=\frac{1}{(2\pi)^{3}}\frac{m^{2}}{\sqrt{E_{\vec{p}^{\prime}_{1}}E_{\vec{p}^{\prime}_{2}}E_{\vec{p}_{1}}E_{\vec{p}_{2}}}}\\\
\times\delta\left(\vec{p}^{\prime}_{1}+\vec{p}^{\prime}_{2}-\vec{p}_{1}-\vec{p}_{2}\right)v_{b}(1^{\prime},2^{\prime};1,2),$
(18) $v_{s}(1^{\prime},2^{\prime};1,2)\\\
=-\frac{g_{s}^{2}}{2}\bar{u}(\vec{p}^{\prime}_{1})u(\vec{p}_{1})\frac{1}{(p_{1}-p^{\prime}_{1})^{2}-m_{s}^{2}}\bar{u}(\vec{p}^{\prime}_{2})u(\vec{p}_{2}),$
(19) $v_{ps}(1^{\prime},2^{\prime};1,2)\\\
=\frac{g_{ps}^{2}}{2}\bar{u}(\vec{p}^{\prime}_{1})\gamma_{5}u(\vec{p}_{1})\frac{1}{(p_{1}-p^{\prime}_{1})^{2}-m_{ps}^{2}}\bar{u}(\vec{p}^{\prime}_{2})\gamma_{5}u(\vec{p}_{2}),$
(20)
$v_{\rm{v}}(1^{\prime},2^{\prime};1,2)=\frac{1}{2}\frac{1}{(p^{\prime}_{1}-p_{1})^{2}-m_{\rm{v}}^{2}}\\\
\times\left[\bar{u}(\vec{p}^{\prime}_{1})\left\\{(g_{\rm{v}}+f_{\rm{v}})\gamma_{\nu}-\frac{f_{\rm{v}}}{2m}(p^{\prime}_{1}+p_{1})_{\nu}\right\\}u(\vec{p}_{1})\right.\\\
\times\bar{u}(\vec{p}^{\prime}_{2})\left\\{(g_{\rm{v}}+f_{\rm{v}})\gamma^{\nu}-\frac{f_{\rm{v}}}{2m}(p^{\prime}_{2}+p_{2})^{\nu}\right\\}u(\vec{p}_{2})\\\
-\bar{u}(\vec{p}^{\prime}_{1})\left\\{(g_{\rm{v}}+f_{\rm{v}})\gamma_{\nu}-\frac{f_{\rm{v}}}{2m}(p^{\prime}_{1}+p_{1})_{\nu}\right\\}u(\vec{p}_{1})\\\
\times\bar{u}(\vec{p}^{\prime}_{2})\frac{f_{\rm{v}}}{2m}\left\\{(\hat{p}_{1}^{\prime}+\hat{p}_{2}^{\prime}-\hat{p}_{1}-\hat{p}_{2})\gamma^{\nu}\right.\\\
\left.\phantom{\frac{f_{\rm{v}}}{2m}}\left.-(p^{\prime}_{1}+p^{\prime}_{2}-p_{1}-p_{2})^{\nu}\right\\}u(\vec{p}_{2})\right],$
(21)

where $m_{b}$ the mass of the clothed boson (its physical value) and
$\hat{q}=q_{\mu}\gamma^{\mu}$. In the framework of the isospin formalism one
needs to add the factor $\vec{\tau}(1)\vec{\tau}(2)$ in the corresponding
expressions.

At this point, our derivation of the vector-boson contribution (21) is to be
specifically commented. Actually, it is the case, where for a
Lorentz–invariant Lagrangian it is not necessarily to have ”… the interaction
Hamiltonian as the integral over space of a scalar interaction density; we
also need to add non–scalar terms to the interaction density …” (quoted from
p.292 of Ref. WeinbergBook1995 ). Let us recall that the density in question
has the property,

$U_{F}(\Lambda,a)\mathscr{H}(x)U_{F}^{-1}(\Lambda,a)=\mathscr{H}(\Lambda
x+a),$ (22)

where the operators $U_{F}(\Lambda,a)$ realize a unitary irreducible
representation of the Poincaré group in the Hilbert space of states for free
(non–interacting) fields.

By definition, the first clothing transformation $W^{(1)}=\exp[R^{(1)}]$
(${R^{(1)}}^{{\dagger}}=-R^{(1)}$) eliminates all interactions linear in the
coupling constants, viz.,

$V^{(1)}=V_{s}+V_{ps}+V_{\rm{v}}^{(1)},$

with

$V_{\rm{v}}^{(1)}=\int
d\vec{x}\left\\{g_{\rm{v}}\bar{\psi}(\vec{x})\gamma_{\mu}\psi(\vec{x})\varphi_{\rm{v}}^{\mu}(\vec{x})\phantom{\frac{f_{\rm{v}}}{4m}}\right.\\\
\left.+\frac{f_{\rm{v}}}{4m}\bar{\psi}(\vec{x})\sigma_{\mu\nu}\psi(\vec{x})\varphi_{\rm{v}}^{\mu\nu}(\vec{x})\right\\}.$
(23)

Following Ref.SheShi01  we have

$R^{(1)}=-i\lim\limits_{\varepsilon\rightarrow
0+}\int\limits_{0}^{\infty}V_{D}^{(1)}(t)e^{-\varepsilon t}dt$ (24)

if $m_{b}<2m$. Here

$V_{D}^{(1)}\equiv\exp[iH_{F}t]V^{(1)}\exp[-iH_{F}t]=\int\mathscr{H}^{(1)}(x)d\vec{x},$

where $\mathscr{H}^{(1)}(x)$ is the Lorentz scalar.

The corresponding interaction operator in the CPR (12) can be written as

$K_{I}^{(2)}(\alpha_{c})=\frac{1}{2}\left[R^{(1)}(\alpha_{c}),V^{(1)}(\alpha_{c})\right]+V^{(2)}(\alpha_{c}),$
(25)

where we have kept only the contributions of the second order in the coupling
constants, so

$V^{(2)}=\int
d\vec{x}\left\\{\frac{g_{\rm{v}}^{2}}{2m_{\rm{v}}^{2}}\bar{\psi}(x\vec{)}\gamma_{0}\psi(\vec{x})\bar{\psi}(\vec{x})\gamma_{0}\psi(\vec{x})\right.\\\
\left.+\frac{f_{\rm{v}}^{2}}{4m^{2}}\bar{\psi}(\vec{x})\sigma_{0i}\psi(\vec{x})\bar{\psi}(\vec{x})\sigma_{0i}\psi(\vec{x})\right\\}.$
(26)

We point out that all quantities in the r.h.s. of Eq.(25) depend on the new
creation(destruction) operators $\alpha_{c}$. In particular, it means that in
the standard Fourier expansions of the fields involved in the definitions of
$V^{(1)}$ and $V^{(2)}$ one should replace the set $\\{\alpha\\}$ by the set
$\\{\alpha_{c}\\}$. Thus, there is an essential distinction between
$V^{(1)}$($V^{(2)}$), on the one hand, and the first(second) integral in the
r.h.s. of Eq.(5), on the other hand.

For this exposition we do not intend to derive all interactions between the
clothed mesons and nucleons, allowed by formula (25). Our aim is more humble,
viz., to find in the r.h.s. of Eq.(25) terms of the type (17), responsible for
the $N$–$N$ interaction. Meanwhile, in case of the vector mesons we encounter
an interplay between the commutator $[R^{(1)}_{\rm{v}},V_{\rm{v}}^{(1)}]/2$
and the integral (26). Indeed, after a simple algebra we find

$\frac{1}{2}\left[R^{(1)},V^{(1)}\right]_{\rm{v}}(NN\rightarrow NN)\\\
=K_{\rm{v}}(NN\rightarrow NN)+K_{cont}(NN\rightarrow NN),$

where the first term has the structure of Eq.(17) with the coefficients by
(21). At the same time the second term $K_{cont}$ completely cancels the
non–scalar operator $V^{(2)}$. The latter may be associated with a contact
interaction since it does not contain any propagators (cf. the approach by the
Osaka group TamuraSato88 ), being expressed through the
$b^{{\dagger}}_{c}(b_{c})$. In other words, the first UCT enables us to remove
the non–invariant terms directly in the Hamiltonian. In our opinion, such a
cancellation, first discussed here, is a pleasant feature of the CPR.

Moreover, as it was shown in Ref.SheShi01 , for each boson included the
corresponding relativistic and properly symmetrized $N$–$N$ interaction, the
kernel of integral equations for the $N$–$N$ bound and scattering states, is
determined by

$\left\langle
b_{c}^{\dagger}(\vec{p}^{\prime}_{1})b_{c}^{\dagger}(\vec{p}^{\prime}_{2})\Omega\right|K_{b}(NN\rightarrow
NN)\left|b_{c}^{\dagger}(\vec{p}_{1})b_{c}^{\dagger}(\vec{p}_{2})\Omega\right\rangle\\\
=V_{b}^{dir}(1^{\prime},2^{\prime};1,2)-V_{b}^{exc}(1^{\prime},2^{\prime};1,2),$
(27)

where we have separated the so–called direct

$V_{b}^{dir}(1^{\prime},2^{\prime};1,2)=-V_{b}(1^{\prime},2^{\prime};1,2)-V_{b}(2^{\prime},1^{\prime};2,1)$
(28)

and exchange

$V_{b}^{exc}(1^{\prime},2^{\prime};1,2)=V_{b}^{dir}(2^{\prime},1^{\prime};1,2)$
(29)

terms. For example, the one–pion–exchange contribution can be divided into the
two parts:

$V_{\pi}^{dir}(1^{\prime},2^{\prime};1,2)=-\frac{g_{\pi}^{2}}{(2\pi)^{3}}\frac{m^{2}}{\sqrt{E_{\vec{p}^{\prime}_{1}}E_{\vec{p}^{\prime}_{2}}E_{\vec{p}_{1}}E_{\vec{p}_{2}}}}\\\
\times\delta\left(\vec{p}^{\prime}_{1}+\vec{p}^{\prime}_{2}-\vec{p}_{1}-\vec{p}_{2}\right)\bar{u}(\vec{p}^{\prime}_{1})\gamma_{5}u(\vec{p}_{1})\bar{u}(\vec{p}^{\prime}_{2})\gamma_{5}u(\vec{p}_{2})\\\
\times\frac{1}{2}\left\\{\frac{1}{(p_{1}-p^{\prime}_{1})^{2}-m_{\pi}^{2}}+\frac{1}{(p_{2}-p^{\prime}_{2})^{2}-m_{\pi}^{2}}\right\\}$
(30)

and

$V_{\pi}^{exc}(1^{\prime},2^{\prime};1,2)=-\frac{g_{\pi}^{2}}{(2\pi)^{3}}\frac{m^{2}}{\sqrt{E_{\vec{p}^{\prime}_{1}}E_{\vec{p}^{\prime}_{2}}E_{\vec{p}_{1}}E_{\vec{p}_{2}}}}\\\
\times\delta\left(\vec{p}^{\prime}_{1}+\vec{p}^{\prime}_{2}-\vec{p}_{1}-\vec{p}_{2}\right)\bar{u}(\vec{p}^{\prime}_{1})\gamma_{5}u(\vec{p}_{2})\bar{u}(\vec{p}^{\prime}_{2})\gamma_{5}u(\vec{p}_{1})\\\
\times\frac{1}{2}\left\\{\frac{1}{(p_{2}-p^{\prime}_{1})^{2}-m_{\pi}^{2}}+\frac{1}{(p_{1}-p^{\prime}_{2})^{2}-m_{\pi}^{2}}\right\\}$
(31)

to be depicted in Fig.1, where the dashed lines correspond to the following
Feynman–like ”propagators”:

$\frac{1}{2}\left\\{\frac{1}{(p_{1}-p^{\prime}_{1})^{2}-m_{\pi}^{2}}+\frac{1}{(p_{2}-p^{\prime}_{2})^{2}-m_{\pi}^{2}}\right\\}$

on the left panel and

$\frac{1}{2}\left\\{\frac{1}{(p_{2}-p^{\prime}_{1})^{2}-m_{\pi}^{2}}+\frac{1}{(p_{1}-p^{\prime}_{2})^{2}-m_{\pi}^{2}}\right\\}$

on the right panel. Other distinctive features of the result (27) have been
discussed in SheShi01 ; KorCanShe07 .

Figure 1: The Feynman–like diagrams for the direct and exchange contributions
in the r.h.s. of Eq.(27).

##  3  The field–theoretic description of the elastic N–N scattering

###  3.1  The $T$–matrix in the CPR

In order to evaluate the $N$–$N$ scattering amplitude for the collision energy
$E$ we will regard a field operator $T$ that meets the equation

$T(E+i0)=H_{I}+H_{I}(E+i0-H_{F})^{-1}T(E+i0)$ (32)

and whose matrix elements $\langle NN|T(E+i0)|NN\rangle$ on the energy shell
$E=E_{1}+E_{2}=E_{1}^{{}^{\prime}}+E_{2}^{{}^{\prime}}$ can be expressed
through the phase shifts and mixing parameters.

Unlike nonrelativistic quantum mechanics (NQM) in relativistic QFT the
interaction $H_{I}$ does not conserve the particle number, being the spring of
particle creation and destruction. The feature makes the problem of finding
the $N$–$N$ scattering matrix much more complicated than in the framework of
nonrelativistic approach.

Such a general field–theoretic consideration can be simplified with the help
of an equivalence theorem ShebekoFB17  according to which the $S$ matrix
elements in the Dirac (D) picture, viz.,

$S_{fi}\equiv\langle\alpha^{{\dagger}}...\Omega_{0}|S(\alpha)|\alpha^{{\dagger}}...\Omega_{0}\rangle$
(33)

are equal to the corresponding elements

$S_{fi}^{c}\equiv\langle\alpha_{c}^{{\dagger}}...\Omega|S(\alpha_{c})|\alpha_{c}^{{\dagger}}...\Omega\rangle$
(34)

of the $S$ matrix in the CPR once the UCTs

$W_{D}(t)=\exp(iK_{F}t)W\exp(-iK_{F}t)$

obey the condition

$W_{D}(\pm\infty)=1$ (35)

The $T$ operator in the CPR satisfies the equation

$T_{cloth}(E+i0)=K_{I}\\\ +K_{I}(E+i0-K_{F})^{-1}T_{cloth}(E+i0)$ (36)

and the matrix

$T_{fi}\equiv\langle f;b|T(E+i0)|i;b\rangle\\\ =\langle
f;c|T_{cloth}(E+i0)|i;c\rangle\equiv T_{fi}^{c},$ (37)

where $|;b\rangle$ ( $|;c\rangle$ ) are the $H_{F}$ ( $K_{F}$ ) eigenvectors,
may be evaluated relying upon properties of the new interaction
$K_{I}(\alpha_{c})$.

If in Eq.(36) we approximate $K_{I}$ by $K_{I}^{(2)}$, then initial task of
evaluating the CPR matrix elements can be reduced to solving the equation

$\langle 1^{\prime},2^{\prime}|T_{NN}(E)|1,2\rangle=\langle
1^{\prime},2^{\prime}|K_{NN}|1,2\rangle\\\ +\langle
1^{\prime},2^{\prime}|K_{NN}(E+i0-K_{F})^{-1}T_{NN}(E)|1,2\rangle.$ (38)

###  3.2  The $R$–matrix equation and its angular–momentum decomposition

For practical applications one prefers to work with the corresponding
$R$–matrix that meets the set of equations

$\left\langle{1^{\prime}2^{\prime}}\right|\bar{R}(E)\left|{12}\right\rangle=\left\langle{1^{\prime}2^{\prime}}\right|\bar{K}_{NN}\left|{12}\right\rangle\\\
+\int\limits_{34}\\!\\!\\!\\!\\!\\!\\!\\!\sum{\left\langle{1^{\prime}2^{\prime}}\right|\bar{K}_{NN}\left|{34}\right\rangle\frac{{\left\langle{34}\right|\bar{R}(E)\left|{12}\right\rangle}}{{E-E_{3}-E_{4}}}}$
(39)

with $\bar{R}(E)=R(E)/2$ and $\bar{K}_{NN}=K_{NN}/2$, where the operation
$\int\limits_{34}\\!\\!\\!\\!\\!\\!\sum$ means the summation over nucleon
polarizations and the $p.v.$ integration over nucleon momenta. The kernel of
Eq.(39) is

$\left\langle{1^{\prime}2^{\prime}}\left|\bar{K}_{NN}\right|{12}\right\rangle=\delta\left(\vec{p}^{\prime}_{1}+\vec{p}^{\prime}_{2}-\vec{p}_{1}-\vec{p}_{2}\right)\left\langle{1^{\prime}2^{\prime}}\left|\bar{V}\right|{12}\right\rangle\\\
\equiv\delta\left(\vec{p}^{\prime}_{1}+\vec{p}^{\prime}_{2}-\vec{p}_{1}-\vec{p}_{2}\right)\\\
\times\left\langle{\vec{p}^{\prime}_{1}\mu^{\prime}_{1}\tau^{\prime}_{1},\vec{p}^{\prime}_{2}\mu^{\prime}_{2}\tau^{\prime}_{2}\left|{\bar{V}}\right|\vec{p}_{1}\mu_{1}\tau_{1},\vec{p}_{2}\mu_{2}\tau_{2}}\right\rangle$

The subsequent calculations are essentially simplified in the center–of–mass
system (c.m.s) in which

$\left\langle\vec{p}^{\prime}\mu^{\prime}_{1}\mu^{\prime}_{2},\tau^{\prime}_{1}\tau^{\prime}_{2}\left|{\bar{R}(E)}\right|\vec{p}\mu_{1}\mu_{2},\tau_{1}\tau_{2}\right\rangle\\\
=\left\langle{\vec{p}^{\prime}\mu^{\prime}_{1}\mu^{\prime}_{2},\tau^{\prime}_{1}\tau^{\prime}_{2}\left|{\bar{V}}\right|\vec{p}\mu_{1}\mu_{2},\tau_{1}\tau_{2}}\right\rangle\\\
+\sum p.v.\int
d\vec{q}\left\langle\vec{p}^{\prime}\mu^{\prime}_{1}\mu^{\prime}_{2},\tau^{\prime}_{1}\tau^{\prime}_{2}\left|{\bar{V}}\right|\vec{q}\mu_{3}\mu_{4},\tau_{3}\tau_{4}\right\rangle\\\
\times\frac{\left\langle\vec{q}\mu_{3}\mu_{4},\tau_{3}\tau_{4}\left|{\bar{R}(E)}\right|\vec{p}\mu_{1}\mu_{2},\tau_{1}\tau_{2}\right\rangle}{E-2E_{\vec{q}}}$
(40)

Here the quantum numbers $\mu(\tau)$ are the individual spin (isospin)
projections.

Accordingly Eq. (27)

$\left\langle{1^{\prime}2^{\prime}}\right|{\bar{V}}\left|{12}\right\rangle=\frac{1}{2(2\pi)^{3}}\frac{m^{2}}{E_{\vec{p}^{\prime}}E_{\vec{p}}}\\\
\times\sum\limits_{b}[v_{b}^{dir}(1^{\prime},2^{\prime};1,2)-v_{b}^{exc}(1^{\prime},2^{\prime};2,1)]$
(41)

with

$v_{b}^{dir}(1^{\prime},2^{\prime};1,2)=-v_{b}(1^{\prime},2^{\prime};1,2)-v_{b}(2^{\prime},1^{\prime};2,1)$
(42)

and

$v_{b}^{exc}(1^{\prime},2^{\prime};2,1)=v_{b}^{dir}(2^{\prime},1^{\prime};1,2),$

where the separate boson contributions are determined by Eqs. (19)–(21) with
$\vec{p}_{1}=\vec{p}=-\vec{p}_{2}$ and
$\vec{p}^{\prime}_{1}=\vec{p}^{\prime}=-\vec{p}^{\prime}_{2}$.

Following a common practice we are interested in the angular–momentum
decomposition of Eq.(40) assuming a nonrelativistic analog of relativistic
partial wave expansions (see Werle66  and refs. therein) for two–particle
states. For example, the clothed two–nucleon state (the so–called two–nucleon
plane wave) can be represented as

$\left|\vec{p}\mu_{1}\mu_{2},\tau_{1}\tau_{2}\right\rangle=\sum\left(\tfrac{1}{2}\mu_{1}\tfrac{1}{2}\mu_{2}\left|SM_{S}\right.\right)\left(\tfrac{1}{2}\tau_{1}\tfrac{1}{2}\tau_{2}\left|TM_{T}\right.\right)\\\
\left(lm_{l}SM_{S}\left|JM_{J}\right.\right)Y_{lm_{l}}^{*}\left(\vec{p}/p\right)\left|pJ(lS)M_{J},TM_{T}\right\rangle,$
(43)

where $J$, $S$ and $T$ are, respectively, total angular momentum, spin and
isospin of the $NN$ pair, being the eigenvalues of the operators
$\vec{J}_{ferm}$, $\vec{S}_{ferm}$ and $\vec{T}_{ferm}$. Here

$\vec{J}_{ferm}=\vec{L}_{ferm}+\vec{S}_{ferm},$ (44)

where $\vec{L}_{ferm}$ ($\vec{S}_{ferm}$) the orbital (spin) momentum of the
fermion field,

$\vec{L}_{ferm}=\frac{i}{2}\sum\limits_{\mu}\int
d\vec{p}\,\,\vec{p}\times\left[\frac{\partial
b_{c}^{\dagger}(\vec{p}\mu)}{\partial\vec{p}}b_{c}(\vec{p}\mu)\right.\\\
-b_{c}^{\dagger}(\vec{p}\mu)\frac{\partial
b_{c}(\vec{p}\mu)}{\partial\vec{p}}+\frac{\partial
d_{c}^{\dagger}(\vec{p}\mu)}{\partial\vec{p}}d_{c}(\vec{p}\mu)\\\
\left.-d_{c}^{\dagger}(\vec{p}\mu)\frac{\partial
d_{c}(\vec{p}\mu)}{\partial\vec{p}}\right]$ (45)

and

$\vec{S}_{ferm}=\frac{1}{2}\sum\limits_{\mu\mu^{\prime}}\int
d\vec{p}\,\,\chi^{{\dagger}}_{\mu^{\prime}}\mathbf{\sigma}\chi_{\mu}\left\\{b_{c}^{\dagger}(\vec{p}\mu^{\prime})b_{c}(\vec{p}\mu)\right.\\\
\left.-d_{c}^{\dagger}(\vec{p}\mu^{\prime})d_{c}(\vec{p}\mu)\right\\},$ (46)

where $\chi_{\mu^{\prime}}(\chi_{\mu})$ are the Pauli spinors. For brevity, we
do not show the isospin operator $\vec{T}_{ferm}$.

The corresponding eigenvalue equations look as

${\vec{J}}^{\,2}_{ferm}\left|{pJ(lS)M_{J}}\right\rangle=J(J+1)\left|{pJ(lS)M_{J}}\right\rangle$
$J_{ferm}^{z}\left|{pJ(lS)M_{J}}\right\rangle=M_{J}\left|{pJ(lS)M_{J}}\right\rangle$
(47)

and

${\vec{S}}^{\,2}_{ferm}\left|{\vec{p}SM_{S}}\right\rangle=S(S+1)\left|{\vec{p}SM_{S}}\right\rangle$
$S_{ferm}^{z}\left|{\vec{p}SM_{S}}\right\rangle=M_{S}\left|{\vec{p}SM_{S}}\right\rangle$
(48)

Doing so, we have introduced the vectors222For a moment, the isospin quantum
numbers are suppressed.

$\left|{\vec{p}SM_{S}}\right\rangle=\sum{\left({\left.\frac{1}{2}\mu_{1}\frac{1}{2}\mu_{2}\right|SM_{S}}\right)}\left|{\vec{p}\mu_{1}\mu_{2}}\right\rangle$
(49)

and

$\left|{pJ(lS)M_{J}}\right\rangle\\\
=\int{d\hat{\vec{p}}}\,Y_{lm_{l}}\left(\vec{p}/p\right)\left|{\vec{p}SM_{S}}\right\rangle\left({lm_{l}SM_{S}\left|{JM_{J}}\right.}\right)$
(50)

A simple way of deriving Eqs.(47)–(48) is to use the transformation

$U_{F}^{c}(R)|\vec{p}SM_{S}\rangle=|R\vec{p}SM_{S}^{\prime}\rangle
D_{M_{S}^{\prime}M_{S}}^{(S)}(R)$ (51) $\forall\,R\,\in\,\mbox{the rotation
group}$

One should note that in our case the separable ansatz

$\left|{\vec{p}_{1}\vec{p}_{2}\mu_{1}\mu_{2}}\right\rangle=\left|{\vec{p}_{1}\mu_{1}}\right\rangle\left|{\vec{p}_{2}\mu_{2}}\right\rangle$

often exploited in relativistic quantum mechanics (RQM) (see, e.g., Werle66
and KeisPoly91 ) does not work. However, one can employ the similarity
transformation 333Sometimes it is convenient to handle the operators
$b_{c}^{{\dagger}}(p\mu)=\sqrt{p_{0}}b_{c}^{{\dagger}}(\vec{p}\mu)$ and their
adjoints $b_{c}(p\mu)$ that meet covariant relations
$\left\\{{b_{c}^{\dagger}(p^{\prime}\mu^{\prime}),b_{c}(p\mu)}\right\\}=p_{0}\delta(\vec{p}^{\prime}-\vec{p})\delta_{\mu^{\prime}\mu}$

$U_{F}^{c}\left({\Lambda,a}\right)b_{c}^{\dagger}\left({p\mu}\right)U_{F}^{c}{\dagger}\left({\Lambda,a}\right)\\\
=e^{i\Lambda p\cdot a}b_{c}^{\dagger}\left({\Lambda
p\mu^{\prime}}\right)D_{\mu^{\prime}\mu}^{(\frac{1}{2})}\left({W\left({\Lambda,p}\right)}\right)$
(52)

with the Wigner rotation $W(\Lambda,p)$ (e.g., for rotations $W(R,p)=R$ ) and
the property of the physical vacuum $\Omega$ to be invariant with respect to
unitary transformations $U_{F}^{c}$ in the CPR (some details can be found in a
separate paper).

The use of expansion (43) gives rise to the well known JST representation, in
which

$\left\langle
p^{\prime}J^{\prime}(l^{\prime}S^{\prime})M^{\prime}_{J},T^{\prime}M^{\prime}_{T}\left|\bar{R}(E)\\{\bar{V}\\}\right|pJ(lS)M_{J},TM_{T}\right\rangle\\\
=\bar{R}(E)\\{\bar{V}\\}_{l^{\prime}l}^{JST}(p^{\prime},p)\delta_{J^{\prime}J}\delta_{M^{\prime}_{J}M_{J}}\delta_{S^{\prime}S}\delta_{T^{\prime}T}\delta_{M^{\prime}_{T}M_{T}},$
(53)

so Eq.(40) reduces to the set of simple integral equations,

$\bar{R}_{l^{\prime}\,l}^{JST}(p^{\prime},p)=\bar{V}_{l^{\prime}\,l}^{JST}(p^{\prime},p)\\\
+\sum\limits_{l^{\prime\prime}}{p.v.}\int\limits_{0}^{\infty}\frac{q^{2}\,dq}{2(E_{p}-E_{q})}\bar{V}_{l^{\prime}\,l^{\prime\prime}}^{JST}(p^{\prime},q)\bar{R}_{l^{\prime\prime}{\rm{}}l}^{JST}(q,p)$
(54)

to be solved for each submatrix $\overline{R}^{JST}$ composed of the elements

$\bar{R}_{l^{\prime}l}^{JST}(p^{\prime},p)\equiv\bar{R}_{l^{\prime}l}^{JST}(p^{\prime},p;2E_{p}),$
(55)

where $E_{p}=\sqrt{\vec{p}^{2}+m^{2}}$ the collision energy in the c.m.s.. One
should note that in view of the charge independence assumed in this work one
has to solve two separate equations for isospin values $T=0$ and $T=1$.

##  4 Results of numerical calculations and their discussion

In the course of our computations we have used the so–called matrix inversion
method (MIM) (see BJ76  and refs. therein). Since we deal with the
relativistic dispersion law for the particle energies, the well known
substraction procedure within the MIM leads to equations

$R_{l^{\prime}{\rm{}}l}^{JST}(p^{\prime},p)=V_{l^{\prime}{\rm{}}l}^{JST}(p^{\prime},p)\\\
+\frac{1}{2}\sum\limits_{l^{\prime\prime}}\int\limits_{0}^{\infty}\frac{dq}{p^{2}-q^{2}}\left\\{q^{2}(E_{p}+E_{q})V_{l^{\prime}l^{\prime\prime}}^{JST}(p^{\prime},q)R_{l^{\prime\prime}l}^{JST}(q,p)\right.\\\
\left.-2p^{2}E_{p}V_{l^{\prime}l^{\prime\prime}}^{JST}(p^{\prime},p)R_{l^{\prime\prime}l}^{JST}(p,p)\right\\}.$
(56)

Table 1: The best–fit parameters for the two models. The third (fourth) column taken from Table A.1 Mach89  (obtained by solving Eqs.(56) with a least squares fitting to OBEP values in Table 2). All masses are in $MeV$, and $n_{b}=1$ except for $n_{\rho}=n_{\omega}=2.$ Meson |  |  Potential B |  UCT
---|---|---|---
$\pi$ |  $g^{2}_{\pi}/4\pi$ |  14.4 |  14.5
|  $\Lambda_{\pi}$ |  1700 |  2200
|  $m_{\pi}$ |  138.03 |  138.03
$\eta$ |  $g^{2}_{\eta}/4\pi$ |  3 |  2.8534
|  $\Lambda_{\eta}$ |  1500 |  1200
|  $m_{\eta}$ |  548.8 |  548.8
$\rho$ |  $g^{2}_{\rho}/4\pi$ |  0.9 |  1.3
|  $\Lambda_{\rho}$ |  1850 |  1450
|  $f_{\rho}/g_{\rho}$ |  6.1 |  5.85
|  $m_{\rho}$ |  769 |  769
$\omega$ |  $g^{2}_{\omega}/4\pi$ |  24.5 |  27
|  $\Lambda_{\omega}$ |  1850 |  2035.59
|  $m_{\omega}$ |  782.6 |  782.6
$\delta$ |  $g^{2}_{\delta}/4\pi$ |  2.488 |  1.6947
|  $\Lambda_{\delta}$ |  2000 |  2200
|  $m_{\delta}$ |  983 |  983
$\sigma,\,T=0$ |  $g^{2}_{\sigma}/4\pi$ |  18.3773 |  19.4434
|  $\Lambda_{\sigma}$ |  2000 |  1538.13
|  $m_{\sigma}$ |  720 |  717.7167
$\sigma,\,T=1$ |  $g^{2}_{\sigma}/4\pi$ |  8.9437 |  10.8292
|  $\Lambda_{\sigma}$ |  1900 |  2200
|  $m_{\sigma}$ |  550 |  568.8612

To facilitate comparison with some derivations and calculations from Refs.
MachHolElst87 , Mach89 , we introduce the notation

$\left\langle{\vec{p}^{\prime}\,\mu^{\prime}_{1}\mu^{\prime}_{2}}\right|v^{UCT}_{b}\left|{\vec{p}\,\mu_{1}\mu_{2}}\right\rangle\\\
\equiv-F^{2}_{b}(p^{\prime},p)\left[v_{b}(1^{\prime},2^{\prime};1,2)+v_{b}(2^{\prime},1^{\prime};2,1)\right]$

for the regularized UCT quasipotentials in the c.m.s. As in Ref.MachHolElst87
, we put that invariants $F_{b}(p^{\prime},p)=F_{b}(\Lambda p^{\prime},\Lambda
p)$ have a phenomenological form,

$F_{b}(p^{\prime},p)=\left[\frac{\Lambda_{b}^{2}-m_{b}^{2}}{\Lambda_{b}^{2}-(p^{\prime}-p)^{2}}\right]^{n_{b}}\equiv
F_{b}[(p^{\prime}-p)^{2}]$

Doing so, we have

$\left\langle\vec{p}^{\,\prime}\,\mu^{\prime}_{1}\mu^{\prime}_{2}\right|v^{UCT}_{s}\left|\vec{p}\,\mu_{1}\mu_{2}\right\rangle\\\
=g_{s}^{2}\bar{u}(\vec{p}^{\,\prime})u(\vec{p})\frac{F^{2}_{s}[(p^{\prime}-p)^{2}]}{(p^{\prime}-p)^{2}-m_{s}^{2}}\bar{u}(-\vec{p}^{\,\prime})u(-\vec{p}),$
(57)
$\left\langle\vec{p}^{\,\prime}\,\mu^{\prime}_{1}\mu^{\prime}_{2}\right|v^{UCT}_{ps}\left|\vec{p}\,\mu_{1}\mu_{2}\right\rangle\\\
=-g_{ps}^{2}\bar{u}(\vec{p}^{\,\prime})\gamma_{5}u(\vec{p})\frac{F^{2}_{ps}[(p^{\prime}-p)^{2}]}{(p^{\prime}-p)^{2}-m_{ps}^{2}}\bar{u}(-\vec{p}^{\,\prime})\gamma_{5}u(-\vec{p})$
(58)

Table 2: Neutron–proton phase shifts (in degrees) for various laboratory energies (in MeV). The OBEP(OBEP∗)–rows taken from Table 5.2 Mach89  (calculated by solving Eqs.(61) with the model parameters from the third column in Table 1). The UCT∗(UCT)–rows calculated by solving Eqs.(56) with the parameters from the third (fourth) column in Table 1. As in MachHolElst87 , we have used the bar convention Stapp  for the phase parameters. State |  Potential |  25 |  50 |  100 |  150 |  200 |  300
---|---|---|---|---|---|---|---
|  OBEP |  50.72 |  39.98 |  25.19 |  14.38 |  5.66 |  -8.18
${{}^{1}}S_{0}$ |  OBEP∗ |  50.71 |  39.98 |  25.19 |  14.37 |  5.66 |  -8.18
|  UCT∗ |  66.79 |  53.01 |  36.50 |  25.27 |  16.54 |  3.12
|  UCT |  50.03 |  39.77 |  25.55 |  15.20 |  6.92 |  -6.07
|  OBEP |  -7.21 |  -11.15 |  -16.31 |  -20.21 |  -23.47 |  -28.70
${{}^{1}}P_{1}$ |  OBEP∗ |  -7.17 |  -11.15 |  -16.32 |  -20.21 |  -23.48 |  -28.71
|  UCT∗ |  -7.40 |  -11.70 |  -17.73 |  -22.63 |  -26.98 |  -34.54
|  UCT |  -7.15 |  -10.95 |  -15.62 |  -18.90 |  -21.49 |  -25.41
|  OBEP |  0.68 |  1.58 |  3.34 |  4.94 |  6.21 |  7.49
${{}^{1}}D_{2}$ |  OBEP∗ |  0.68 |  1.58 |  3.34 |  4.94 |  6.21 |  7.49
|  UCT∗ |  0.68 |  1.59 |  3.40 |  5.10 |  6.52 |  8.20
|  UCT |  0.68 |  1.56 |  3.22 |  4.68 |  5.77 |  6.68
|  OBEP |  9.34 |  12.24 |  9.80 |  4.57 |  -1.02 |  -11.48
${{}^{3}}P_{0}$ |  OBEP∗ |  9.34 |  12.24 |  9.80 |  4.57 |  -1.02 |  -11.48
|  UCT∗ |  9.48 |  12.53 |  10.32 |  5.27 |  -0.15 |  -10.27
|  UCT |  9.30 |  12.16 |  9.81 |  4.73 |  -0.68 |  -10.76
|  OBEP |  -5.33 |  -8.77 |  -13.47 |  -17.18 |  -20.49 |  -26.38
${{}^{3}}P_{1}$ |  OBEP∗ |  -5.33 |  -8.77 |  -13.47 |  -17.18 |  -20.48 |  -26.38
|  UCT∗ |  -5.27 |  -8.62 |  -13.09 |  -16.56 |  -19.63 |  -25.06
|  UCT |  -5.28 |  -8.58 |  -12.85 |  -16.06 |  -18.86 |  -23.79
|  OBEP |  3.88 |  9.29 |  17.67 |  22.57 |  24.94 |  25.36
${{}^{3}}D_{2}$ |  OBEP∗ |  3.89 |  9.29 |  17.67 |  22.57 |  24.94 |  25.36
|  UCT∗ |  3.86 |  9.15 |  17.12 |  21.51 |  23.47 |  23.48
|  UCT |  3.89 |  9.25 |  17.31 |  21.77 |  23.75 |  23.61
|  OBEP |  80.32 |  62.16 |  41.99 |  28.94 |  19.04 |  4.07
${{}^{3}}S_{1}$ |  OBEP∗ |  80.31 |  62.15 |  41.98 |  28.93 |  19.03 |  4.06
|  UCT∗ |  92.30 |  72.71 |  51.44 |  38.10 |  28.20 |  13.70
|  UCT |  79.60 |  61.53 |  41.57 |  28.75 |  19.08 |  4.60
|  OBEP |  -2.99 |  -6.86 |  -12.98 |  -17.28 |  -20.28 |  -23.72
${{}^{3}}D_{1}$ |  OBEP∗ |  -2.99 |  -6.87 |  -12.99 |  -17.28 |  -20.29 |  -23.72
|  UCT∗ |  -2.74 |  -6.43 |  -12.36 |  -16.54 |  -19.47 |  -22.78
|  UCT |  -3.00 |  -6.90 |  -13.12 |  -17.66 |  -21.11 |  -26.03
|  OBEP |  1.76 |  2.00 |  2.24 |  2.58 |  3.03 |  4.03
$\varepsilon_{1}$ |  OBEP∗ |  1.76 |  2.00 |  2.24 |  2.58 |  3.03 |  4.03
|  UCT∗ |  0.02 |  -0.12 |  -0.17 |  0.04 |  0.41 |  1.40
|  UCT |  1.80 |  2.01 |  2.19 |  2.50 |  2.90 |  3.83
|  OBEP |  2.62 |  6.14 |  11.73 |  14.99 |  16.65 |  17.40
${{}^{3}}P_{2}$ |  OBEP∗ |  2.62 |  6.14 |  11.73 |  14.99 |  16.65 |  17.39
|  UCT∗ |  2.80 |  6.61 |  12.71 |  16.28 |  18.10 |  18.91
|  UCT |  2.57 |  6.00 |  11.32 |  14.18 |  15.37 |  15.07
|  OBEP |  0.11 |  0.34 |  0.77 |  1.04 |  1.10 |  0.52
${{}^{3}}F_{2}$ |  OBEP∗ |  0.11 |  0.34 |  0.77 |  1.04 |  1.10 |  0.52
|  UCT∗ |  0.11 |  0.34 |  0.77 |  1.05 |  1.13 |  0.64
|  UCT |  0.11 |  0.34 |  0.75 |  1.00 |  1.03 |  0.41
|  OBEP |  -0.86 |  -1.82 |  -2.84 |  -3.05 |  -2.85 |  -2.02
$\varepsilon_{2}$ |  OBEP∗ |  -0.86 |  -1.82 |  -2.84 |  -3.05 |  -2.85 |  -2.02
|  UCT∗ |  -0.87 |  -1.83 |  -2.82 |  -2.99 |  -2.75 |  -1.88
|  UCT |  -0.86 |  -1.83 |  -2.84 |  -3.05 |  -2.89 |  -2.18
Figure 2: Neutron-proton phase parameters for the uncoupled partial waves,
plotted versus the nucleon kinetic energy in the lab. system. Dashed[solid]
curves calculated with Potential B parameters (Table 1) by solving Eqs.
(56)[(61)]. Dotted represent the solutions of Eqs. (56) with UCT parameters
(Table 1). The rhombs show original OBEP results (see Table 2).

$\left\langle{\vec{p}^{\,\prime}\,\mu^{\prime}_{1}\mu^{\prime}_{2}}\right|v^{UCT}_{\rm{v}}\left|{\vec{p}\,\mu_{1}\mu_{2}}\right\rangle=-\frac{F^{2}_{\rm
v}[(p^{\prime}-p)^{2}]}{\left({p^{\prime}-p}\right)^{2}-m_{\rm{v}}^{2}}\\\
\times\left\\{\bar{u}(\vec{p}^{\,\prime})\left[\left(g_{\rm{v}}+f_{\rm{v}}\right)\gamma_{\nu}-\frac{f_{\rm{v}}}{2m}\left(p^{\prime}+p\right)_{\nu}\right.\right.\\\
\left.-\frac{f_{\rm{v}}}{2m}(E_{\vec{p}^{\prime}}-E_{\vec{p}})[\gamma_{0}\gamma_{\nu}-g_{0\nu}]\right]u\left(\vec{p}\right)\\\
\times\bar{u}\left(-\vec{p}^{\,\prime}\right)\left[\left({g_{\rm{v}}+f_{\rm{v}}}\right)\gamma^{\nu}-\frac{f_{\rm{v}}}{2m}\overline{\left(p^{\prime}+p\right)}^{\nu}\right.\\\
\left.-\frac{f_{\rm{v}}}{2m}(E_{\vec{p}^{\prime}}-E_{\vec{p}})[\gamma^{0}\gamma^{\nu}-g^{0\nu}]\right]u(-\vec{p})$
$-\frac{{f_{\rm
v}}^{2}}{4m^{2}}(E_{p^{\prime}}-E_{p})^{2}\bar{u}(\vec{p}^{\,\prime})[\gamma_{0}\gamma_{\nu}-g_{0\nu}]u(\vec{p})\\\
\left.\phantom{\frac{f_{\rm{v}}}{2m}}\times\bar{u}(-\vec{p}^{\,\prime})[\gamma^{0}\gamma^{\nu}-g^{0\nu}]u(-\vec{p})\right\\},$
(59)

where
$\overline{(p^{\prime}+p)}^{\nu}=(E_{\vec{p}^{\prime}}+E_{\vec{p}},-(\vec{p}^{\prime}+\vec{p}))$.

At first sight, such a regularization can be achieved via a simple
substitution $g_{b}\rightarrow g_{b}F_{b}(p^{\prime},p)$ with some cutoff
functions $F_{b}(p^{\prime},p)$ depending on the 4–momenta $p^{\prime}$ and
$p$. However, the principal moment is to satisfy the requirement (22) for the
Hamiltonian invariant under space inversion, time reversal and charge
conjugation. In this context, let us remind that the baryon–nucleon–nucleon
form factors are expressed through the matrix elements $\langle
p^{\prime}|j_{b}(0)|p\rangle$ of the corresponding baryon current density
$j_{b}(x)$ at $x=0$ between physical(clothed) one–nucleon states Gas . Such
matrix elements might be evaluated in terms of the cutoffs
$F_{b}(p^{\prime},p)$ using some idea from SheShi00  (cf. the clothed particle
representation of a current therein).

Figure 3: The same in Fig. 2 but for the coupled waves.

Replacing in equations (57)–(59)

$\frac{1}{(p^{\prime}-p)^{2}-m_{b}^{2}}\,F^{2}_{b}[(p^{\prime}-p)^{2}]$

by

$\frac{-1}{(\vec{p}^{\prime}-\vec{p})^{2}+m_{b}^{2}}\,F^{2}_{b}[-(\vec{p}^{\prime}-\vec{p})^{2}]$

and neglecting the tensor-tensor term

$\frac{{f_{\rm v}}^{2}}{4m^{2}}(E_{p^{\prime}}-E_{p})^{2}\\\
\times\bar{u}(\vec{p}^{\,\prime})[\gamma_{0}\gamma_{\nu}-g_{0\nu}]u(\vec{p})\bar{u}(-\vec{p}^{\,\prime})[\gamma^{0}\gamma^{\nu}-g^{0\nu}]u(-\vec{p})$
(60)

in (59), we obtain approximate expressions that with the common factor

$(2\pi)^{-3}m^{2}/E_{p^{\prime}}E_{p}$

instead of

$(2\pi)^{-3}m/\sqrt{E_{p^{\prime}}E_{p}}$

are equivalent to Eqs. (E.21)–(E.23) from MachHolElst87 . Such an equivalence
becomes coincidence if in our formulae instead of the canonical two-nucleon
basis $\left|{\vec{p}\,\mu_{1}\mu_{2}}\right\rangle$ one uses the helicity
basis as in MachHolElst87 .

In the context, we have considered the set of equations

${{}^{B}}R_{l^{\prime}{\rm{}}l}^{JST}(p^{\prime},p)={{}^{B}}V_{l^{\prime}{\rm{}}l}^{JST}(p^{\prime},p)\\\
+m\sum\limits_{l^{\prime\prime}}\int\limits_{0}^{\infty}\frac{dq}{p^{2}-q^{2}}\left\\{q^{2}\,{{}^{B}}V_{l^{\prime}l^{\prime\prime}}^{JST}(p^{\prime},q){{}^{B}}R_{l^{\prime\prime}l}^{JST}(q,p)\right.\\\
\left.-p^{2}\,{{}^{B}}V_{l^{\prime}l^{\prime\prime}}^{JST}(p^{\prime},p){{}^{B}}R_{l^{\prime\prime}l}^{JST}(p,p)\right\\},$
(61)

where the superscript $B$ refers to the partial matrix elements of the
potential $B$ defined in Mach89  with the just mentioned interchange of the
bases.

Our calculations of the $R$ matrices that meet the equations (56) and (61) are
twofold. On the one hand, we will check reliability of our numerical procedure
(in particular, its code). On the other hand, we would like to show
similarities and discrepancies between our results and those by the Bonn group
both on the energy shell and beyond it. These results are depicted in Figs.
2–3 and collected in Table 2.

As seen in Figs. 2–3, the most appreciable distinctions between the UCT and
OBEP curves take place for the phase shifts with the lowest $l-$values. As the
orbital angular momentum increases the difference between the solid and dashed
curves decreases. Such features may be explained if one takes into account
that the approximations under consideration affect mainly high–momentum
components of the UCT quasipotentials (their behavior at ”small” distances).
With the $l$–increase the influence of small distances is suppressed by the
centrifugal barrier repulsion.

Of course, it would be more instructive to compare the corresponding
half–off–energy–shell $R$–matrices (see definition (55)). Their
$p^{\prime}$–dependencies not shown here have been prepared for a separate
publication. They are necessary to know when calculating the $\psi^{(\pm)}$
scattering states for a two–nucleon system. In the context, one should
emphasize that hitherto we have explored the OBEP and UCT $R$–matrices in the
c.m.s., where the both approaches yield most close results. It is not the case
in those situations when the c.m.s. cannot be referred to everywhere (e.g., in
the reactions $NN\rightarrow\gamma NN$ and $\gamma d\rightarrow pn$). In this
respect our studies of the differences between UCT and OBE approaches are
under way.

##  5 Summary

The present work has been made to develop a consistent field–theoretical
approach in the theory of nucleon–nucleon scattering. It has been shown that
the method of UCT’s, based upon the notion of clothed particles, is proved to
be appropriate in achieving this purpose.

Using the unitary equivalence of the CPR to the BPR, we have seen how in the
approximation $K_{I}=K_{I}^{(2)}$ the extremely complicated scattering problem
in QFT can be reduced to the three–dimensional LS–type equation for the
$T$–matrix in momentum space.The equation kernel is given by the clothed
two–nucleon interaction of the class [2.2]. Such a conversation becomes
possible owing to the property of $K_{I}^{(2)}$ to leave the two–nucleon
sector and its separate subsectors to be invariant.

Special attention has been paid to the elimination of auxiliary field
components. We encounter such a necessity for interacting vector and fermion
fields when in accordance with the canonical formalism the interaction
Hamiltonian density embodies not only a scalar contribution but nonscalar
terms too. It has proved (at least, for the primary $\rho N$ and $\omega N$
couplings) that the UCT method allows us to remove such noncovariant terms
directly in the Hamiltonian. To what extent this result will take place in
higher orders in coupling constants it will be a subject of further
explorations.

## References

  * (1) M. Lacombe et al., Phys. Rev. C21, (1980) 861
  * (2) R. Machleidt, K. Holinde and C. Elster, Phys. Rep. 149, (1987) 1
  * (3) V.G.J. Stocks et al., Phys. Rev. C49, (1994) 2950
  * (4) R.B. Wiringa, V.G.J. Stocks and R. Schiavilla, Phys. Rev. C51, (1995) 38
  * (5) R. Machleidt, Phys. Rev. C63, (2001) 024001
  * (6) F. Gross and A. Stadler, Few Body Syst. 44, (2008) 295
  * (7) C. Ordonez, L. Ray and U. van Kolck, Phys. Rev. Lett. 72, (1994) 1982
  * (8) E. Epelbaum, W. Glöckle and U.-G. Meissner, Nucl. Phys. A671, (2000) 295
  * (9) E. Epelbaum, Prog. Part. Nucl. Phys. 57, (2006) 654
  * (10) A.V. Shebeko and M.I. Shirokov, Phys. Part. Nucl. 32, (2001) 31
  * (11) V. Korda, L. Canton and A. Shebeko, Ann. Phys. 322, (2007) 736
  * (12) S. Weinberg, The Quantum Theory of Fields, (University Press, Cambridge, 1995) Vol. 1
  * (13) D.A. Schütte, Nucl. Phys. A221, (1974) 450
  * (14) K. Tamura, T. Niva, T. Sato and H. Ohtsubo, Prog. Theor. Phys. 80, (1988) 138.
  * (15) A.V. Shebeko, Nucl. Phys. A737, (2004) 252
  * (16) J.Werle, Relativistic Theory of Reactions (PWN – Polish Scientific Publishers, Warszawa, 1966)
  * (17) B.D.Keister, W.N.Polyzou, Adv. Nucl. Phys. 20, (1991) 266
  * (18) G.E. Brown and A.D Jackson, Nucleon–Nucleon Interaction (Amsterdam: North–Holland Publ. Co. 1976)
  * (19) H.Stapp et al., Phys. Rev. 105, (1957) 302
  * (20) R. Machleidt, Adv. Nucl. Phys. 19, (1989) 189
  * (21) S. Gasiorowicz, Elementary Particle Physics, (John Wiley & Sons, New York, 1966)
  * (22) A.V. Shebeko and M.I. Shirokov, Prog. Part. Nucl. Phys. 44, (2000) 75.