# Exchange effects in nucleus-nucleus reactions

J. Dohet-Eraly<EMAIL_ADDRESS>Physique Quantique, C.P. 165/82, Université
Libre de Bruxelles (ULB), B 1050 Brussels, Belgium  Physique Nucléaire
Théorique et Physique Mathématique, C.P. 229, Université Libre de Bruxelles
(ULB), B 1050 Brussels, Belgium  P. Descouvemont<EMAIL_ADDRESS>Physique
Nucléaire Théorique et Physique Mathématique, C.P. 229, Université Libre de
Bruxelles (ULB), B 1050 Brussels, Belgium

###### Abstract

We present a scattering model for nuclei with similar masses. In this three-
body model, the projectile has a core+valence structure, whereas the target is
identical to the core nucleus. The three-body wave functions must be
symmetrized for the exchange of the cores. This property gives rise to non-
local potentials, which are computed without approximation. The present model
is an extension of the Continuum Discretized Coupled Channel (CDCC) formalism,
with an additional treatment of core exchange. We solve the coupled-channel
system, including non-local terms, by the $R$-matrix method using Lagrange
functions. This model is applied to the ${}^{13}{\rm C}+^{12}$C, ${}^{13}{\rm
N}+^{12}$C and ${}^{16}{\rm O}+^{12}$C systems. Experimental scattering cross
sections are fairly well reproduced without any parameter fitting. The
backward-angle enhancement of the elastic cross sections is due to the non-
local potential. We discuss in more detail the various non-local contributions
and present effective local potentials.

##  I Introduction

Nucleus-nucleus reactions represent an important topic in nuclear physics. In
particular, they constitute the only way to investigate exotic nuclei. With
the development of radioactive beams, more and more data become available.
Accurate theoretical models are needed to interpret these data and to extract
the relevant properties of exotic nuclei.

A popular approach is the optical model Feshbach (1958); Dickhoff and Charity
(2019), where the structure of the colliding nuclei is neglected. Microscopic
effects and absorption channels are simulated by complex potentials. This
approach is very simple, but usually involves several parameters. The
information about the structure of the nuclei is therefore limited.

Three-body models represent a step further in the description of nucleus-
nucleus collisions. One of the participating nuclei is described by a two-body
structure, and the main part of the absorption is simulated by breakup effects
in this two-body nucleus. This approach is referred to as the Continuum
Discretized Coupled Channel (CDCC) method Rawitscher (1974); Kamimura  _et
al._ (1986); Austern  _et al._ (1987); Yahiro  _et al._ (2012), and has been
extended to systems involving four-body systems Matsumoto  _et al._ (2004);
Descouvemont (2018). It is well adapted to nuclei with a low separation
energy, where breakup effects are expected to be important. The CDCC method
was originally developed to describe deuteron scattering Rawitscher (1974),
but many applications have been performed recently for reactions involving
exotic nuclei (see, e.g., Refs. de Diego  _et al._ (2014); Pesudo  _et al._
(2017) for recent works).

In its present form, the CDCC method neglects possible exchange effects
between the projectile and the target. A typical example is the
$\alpha+^{8}$Be reaction Ogata  _et al._ (2009), where the symmetrization
between the colliding $\alpha$ particle and the $\alpha$’s involved in 8Be is
not taken into account. A more recent example is the $d+^{11}$Be system
Descouvemont (2017), where $d$ and 11Be are described by $p+n$ and
${}^{10}{\rm Be}+n$ structures, without antisymmetrization between the
neutrons of $d$ and of 11Be.

An obvious situation where exchange effects are important is when the
colliding nuclei have similar masses. Representative examples are the
${}^{13}{\rm C}+^{12}$C and ${}^{17}{\rm O}+^{16}$O reactions. In such a case
the system can be described by a three-body structure involving two cores and
an exchanged particle (typically a nucleon or an $\alpha$ particle). The
symmetrization of the wave function for the core exchange is then crucial. In
the literature, several works have been done in this direction, with various
approximations of the exchange effects von Oertzen (1970); Imanishi and von
Oertzen (1987); Buttle and Goldfarb (1966); Sparenberg  _et al._ (2000).

In the present work, we use a three-body model, and treat exchange effects
exactly. This procedure gives rise to non-local potentials in a coupled-
channel formalism, but does not require any parameter fit. As in the
traditional CDCC approach, the only inputs are the two-body interactions
between the constituents. A first important step is to determine the non-local
potentials, stemming from exchange effects. In a second step, one has to solve
a coupled-channel integro-differential system. This is in general a
complicated task, but can be simplified with the help of the $R$-matrix
formalism Descouvemont and Baye (2010) associated with the Lagrange-mesh
technique Baye (2015).

The paper is organized as follows. In Sec. II, we present the model, with
emphasis on the calculation of the non-local terms. Section III is devoted to
some applications. We present results on ${}^{13}{\rm C}+^{12}$C, ${}^{13}{\rm
N}+^{12}$C and ${}^{16}{\rm O}+^{12}$C scattering. In Sec. IV, we discuss non-
local effects in more detail. We focus on the long-range part of the non-local
kernels. We also present equivalent local potentials. Concluding remarks and
outlook are presented in Sec. V.

##  II The three-body model

###  II.1 Total wave functions

We consider the three-body system presented in Fig. 1. The projectile is
formed by a core ($C$) + valence ($v$) system, and the target is identical to
the core. A typical example is the ${}^{13}{\rm C}+^{12}$C system, where the
core is 12C and the valence particle a neutron. For the sake of simplicity, we
assume that the spin of the core is zero. The Hamiltonian of this system is
defined as

$\displaystyle H=$ $\displaystyle
T_{\boldsymbol{r}}+T_{\boldsymbol{R}}+V_{Cv}(\boldsymbol{r})$
$\displaystyle+V_{Cv}(|\alpha\boldsymbol{r}-\boldsymbol{R}|)+V_{CC}(|\beta\boldsymbol{r}+\boldsymbol{R}|),$
(1)

where $T_{\boldsymbol{r}}$ and $T_{\boldsymbol{R}}$ are the kinetic energies
associated with $\boldsymbol{r}$ and $\boldsymbol{R}$, and $V_{Cv}$ and
$V_{CC}$ are core-valence and core-core potentials. The coordinates
$\boldsymbol{R}$ and $\boldsymbol{R}^{\prime}$ are the relative coordinates
between the projectile and the target before and after symmetrization,
respectively. In definition (1), $\alpha$ and $\beta$ are positive
coefficients given by

$\displaystyle\alpha=\frac{A_{C}}{A_{C}+A_{v}},\
\beta=\frac{A_{v}}{A_{C}+A_{v}}=1-\alpha,$ (2)

where $A_{C}$ and $A_{v}$ are the masses of the core and of the valence
particle, respectively. We also define

$\displaystyle\gamma=\frac{1}{1-\alpha^{2}},$ (3)

which will be used later.

The present approach is based on the CDCC formalism, but we include the
symmetrization of the wave function with respect to core exchange. Notice that
two possible choices exist for the potential $V_{Cv}$: either it reproduces
the spectroscopic properties of the $C+v$ system, or it is fitted on elastic
scattering.

Figure 1: Coordinates $(\boldsymbol{r},\boldsymbol{R})$ and
$(\boldsymbol{r}^{\prime},\boldsymbol{R}^{\prime})$. $C$ and $v$ represent the
core and valence particles.

Let us define the core-valence Hamiltonian by

$\displaystyle H_{0}=T_{\boldsymbol{r}}+V_{Cv}(\boldsymbol{r}),$ (4)

which is diagonalized from

$\displaystyle H_{0}\phi^{\ell jm}_{n}(\boldsymbol{r})=E^{\ell
j}_{n}\phi^{\ell jm}_{n}(\boldsymbol{r}).$ (5)

The two-body wave functions are factorized as

$\displaystyle\phi^{\ell jm}_{n}(\boldsymbol{r})=\frac{u^{\ell
j}_{n}(r)}{r}[Y_{\ell}(\Omega_{r})\otimes\chi^{v}]^{jm},$ (6)

where $\chi^{v}$ is a spinor associated with the valence particle. The radial
eigenfunctions $u^{\ell j}_{n}(r)$ are expanded over a basis, chosen here from
Lagrange functions Baye (2015). Lagrange-Laguerre functions regularized by $r$
and by $\sqrt{r}$ Baye (2015); Dohet-Eraly (2017) have been considered. Both
provide the same level of accuracy. As usual in CDCC calculations, energies
$E^{\ell j}_{n}<0$ correspond to physical states, whereas states with $E^{\ell
j}_{n}>0$, referred to as pseudostates, simulate the projectile continuum.

The total wave functions, associated with (1) are written before
symmetrization as

$\displaystyle\Phi^{JM\pi}(\boldsymbol{R},\boldsymbol{r})=\frac{1}{R}\sum_{cL}g^{J\pi}_{cL}(R)\,\varphi^{JM\pi}_{cL}(\Omega_{R},\boldsymbol{r}),$
(7)

where index $c$ stands for $c=(\ell jn)$ and where $g^{J\pi}_{cL}$ are radial
wave functions to be determined. The relative angular momentum is denoted as
$L$. The channel functions $\varphi^{JM\pi}_{cL}$ are defined as

$\displaystyle\varphi^{JM\pi}_{cL}(\Omega_{R},\boldsymbol{r})=\bigl{[}Y_{L}(\Omega_{R})\otimes\phi^{\ell
j}_{n}(\boldsymbol{r})\bigr{]}^{JM}.$ (8)

The total wave function (7) must be symmetrized with respect to the exchange
of the cores. For spin-zero cores, this is achieved with the exchange operator
$P$

$\displaystyle\Psi^{JM\pi}(\boldsymbol{R},\boldsymbol{r})=(1+P)\Phi^{JM\pi}(\boldsymbol{R},\boldsymbol{r}),$
(9)

where $P$ permutes the coordinates of the cores. More precisely, we have

$\displaystyle
P\Phi^{JM\pi}(\boldsymbol{R},\boldsymbol{r})=\Phi^{JM\pi}(\boldsymbol{R}^{\prime},\boldsymbol{r}^{\prime}),$
(10)

with

$\displaystyle\boldsymbol{r}^{\prime}=\alpha\boldsymbol{r}-\boldsymbol{R},$
$\displaystyle\boldsymbol{R}^{\prime}=(\alpha^{2}-1)\boldsymbol{r}-\alpha\boldsymbol{R}.$
(11)

In the $(\boldsymbol{R}^{\prime},\boldsymbol{r}^{\prime})$ system, Hamiltonian
(1) can be written, in the so-called “post” form, as

$\displaystyle
H=T_{\boldsymbol{r}^{\prime}}+T_{\boldsymbol{R}^{\prime}}+V_{Cv}(\boldsymbol{r}^{\prime})$
$\displaystyle\hskip
28.45274pt+V_{Cv}(|\alpha\boldsymbol{r}^{\prime}-\boldsymbol{R}^{\prime}|)+V_{CC}(|\beta\boldsymbol{r}^{\prime}+\boldsymbol{R}^{\prime}|).$
(12)

In the next step, we consider the three-body Schrödinger equation

$\displaystyle H\Psi^{JM\pi}=E\Psi^{JM\pi},$ (13)

and use (9) with expansion (7). After projection on the channel functions,
this procedure provides the integro-differential system

$\displaystyle(T_{R}+E_{c}-E)g^{J\pi}_{cL}(R)+\sum_{c^{\prime}L^{\prime}}V^{J\pi}_{cL,c^{\prime}L^{\prime}}(R)g^{J\pi}_{c^{\prime}L^{\prime}}(R)$
$\displaystyle+\sum_{c^{\prime}L^{\prime}}\int
W^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})g^{J\pi}_{c^{\prime}L^{\prime}}(R^{\prime})dR^{\prime}=0,$
(14)

where

$\displaystyle
T_{R}=-\frac{\hbar^{2}}{2\mu}\biggl{[}\frac{d^{2}}{dR^{2}}-\frac{L(L+1)}{R^{2}}\biggr{]},$
(15)

$\mu$ being the reduced mass. The first two terms of Eq. (14) correspond to
the standard CDCC system Austern  _et al._ (1987). The coupling potentials are
defined by

$\displaystyle
V^{J\pi}_{cL,c^{\prime}L^{\prime}}(R)=\langle\varphi^{JM\pi}_{cL}|V_{Cv}+V_{CC}|\varphi^{JM\pi}_{c^{\prime}L^{\prime}}\rangle,$
(16)

where the integration is performed over $\Omega_{R}$ and $\boldsymbol{r}$. The
last term of (14) is non-local, and arises from the symmetrization operator
$P$. The non-local potential
$W^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})$ can be decomposed as

$\displaystyle W^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})=(E_{c}-E){\cal
N}^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})$ $\displaystyle\hskip
28.45274pt+{\cal T}^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})+{\cal
V}^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime}),$ (17)

which explicitly shows overlap ${\cal N}^{J\pi}_{cL,c^{\prime}L^{\prime}}$,
kinetic energy ${\cal T}^{J\pi}_{cL,c^{\prime}L^{\prime}}$ and potential
${\cal V}^{J\pi}_{cL,c^{\prime}L^{\prime}}$ terms. The overlap kernel is
defined from

$\displaystyle\int{\cal
N}^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})g^{J\pi}_{c^{\prime}L^{\prime}}(R^{\prime})dR^{\prime}=$
$\displaystyle\hskip
28.45274ptR\langle\varphi^{JM\pi}_{cL}|\varphi^{JM\pi}_{c^{\prime}L^{\prime}}\dfrac{g^{J\pi}_{c^{\prime}L^{\prime}}(R^{\prime})}{R^{\prime}}\rangle,$
(18)

and equivalent expressions hold for the kinetic energy and potential kernels.
Notice that similar terms shows up in the coupled-channel approach of transfer
reactions Cotanch (1975). We discuss the various contributions of (17) in the
next subsection.

For scattering states ($E>0$), a radial function $g^{J\pi}_{cL}(R)$ has the
asymptotic behavior at large $R$ values

$\displaystyle g^{J\pi}_{cL,\omega L_{\omega}}(R)\rightarrow{v_{c}}^{-1/2}$
$\displaystyle\times\bigl{[}I_{L}(k_{c}R)\delta_{c\omega}\delta_{LL_{\omega}}-U^{J\pi}_{cL,\omega
L_{\omega}}O_{L}(k_{c}R)\bigr{]},$ (19)

where $\omega$ is the entrance channel, $I_{L}$ and $O_{L}$ are the incoming
and outgoing Coulomb functions, $v_{c}$ and $k_{c}$ are the velocity and wave
number in channel $c$, and $\boldsymbol{U}^{J\pi}$ is the scattering matrix.
In the $R$-matrix formalism Descouvemont and Baye (2010), a channel radius $a$
separates the internal region, where all terms of the potentials contribute,
and the external region where only the monopole part of the Coulomb
interaction is present. In the internal region, the radial wave function is
expanded as

$\displaystyle g^{J\pi}_{cL,\omega
L_{\omega}}(R)=\sum_{i=1}^{N}f^{J\pi}_{cLi,\omega L_{\omega}}u_{i}(R),$ (20)

where the $N$ functions $u_{i}(R)$ represent the basis. The choice of the
basis functions is discussed in Subsec. II.D.

###  II.2 Non-local terms

The non-local potential (17) arises from exchange effects due to the operator
$P$ (9,10). The overlap and potential kernels in (17) are obtained from

$\displaystyle\begin{Bmatrix}{\cal
N}^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})\\\ {\cal
V}^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})\end{Bmatrix}={\cal
J}RR^{\prime}\iint\varphi^{JM\pi\ast}_{cL}(\Omega_{R},\boldsymbol{r})$
$\displaystyle\hskip 28.45274pt\times\begin{Bmatrix}1\\\
V_{Cv}+V_{CC}\end{Bmatrix}\varphi^{JM\pi}_{c^{\prime}L^{\prime}}(\Omega_{R^{\prime}},\boldsymbol{r}^{\prime})d\Omega_{R}\,d\Omega_{R^{\prime}},$
(21)

where

$\displaystyle{\cal J}=\gamma^{3}$ (22)

is the Jacobian from coordinates $(\boldsymbol{r},\boldsymbol{R})$ to
$(\boldsymbol{R},\boldsymbol{R}^{\prime})$. Coordinates
$(\boldsymbol{r},\boldsymbol{r}^{\prime},\boldsymbol{R}_{CC})$ are expressed
as

$\displaystyle\boldsymbol{r}=-\gamma(\alpha\boldsymbol{R}+\boldsymbol{R}^{\prime}),$
$\displaystyle\boldsymbol{r}^{\prime}=-\gamma(\boldsymbol{R}+\alpha\boldsymbol{R}^{\prime}),$
$\displaystyle\boldsymbol{R}_{CC}=\boldsymbol{r}-\boldsymbol{r}^{\prime}=\frac{1}{\alpha+1}(\boldsymbol{R}-\boldsymbol{R}^{\prime}).$
(23)

In the case of a nucleon transfer, $\alpha$ is close to 1, and the relative
coordinates $(\boldsymbol{r},\boldsymbol{r}^{\prime})$ are large, even for
relatively small values of $(\boldsymbol{R},\boldsymbol{R}^{\prime})$. This
means that the core-valence wave function needs to be accurately known up to
large distances. If the binding energy is small, the non-local potentials are
therefore very sensitive to the long-range part of bound-state wave function.

The potential term is quite similar to the matrix elements involved in
Distorted Wave Born Approximation (DWBA) calculations Satchler (1983), and its
calculation is explained in several references (see e.g. Refs. Satchler
(1983); Thompson (1988); Shubhchintak and Descouvemont (2019)). The
calculation of the overlap and potential kernels is based on the expansions

$\displaystyle\frac{u^{\ell
j}_{n}(r)\,u^{\ell^{\prime}j^{\prime}}_{n^{\prime}}(r^{\prime})}{r^{\ell+1}r^{\prime\ell^{\prime}+1}}=\sum_{K}N^{K}_{cc^{\prime}}(R,R^{\prime})P_{K}(\cos\theta_{R}),$
$\displaystyle\frac{u^{\ell
j}_{n}(r)}{r^{\ell+1}}\biggl{(}V_{Cv}(r^{\prime})+V_{CC}(R_{CC})\biggr{)}\frac{u^{\ell^{\prime}j^{\prime}}_{n^{\prime}}(r^{\prime})}{r^{\prime\ell^{\prime}+1}}=$
$\displaystyle\hskip
28.45274pt\sum_{K}V^{K}_{cc^{\prime}}(R,R^{\prime})P_{K}(\cos\theta_{R}),$
(24)

where $\theta_{R}$ is the angle between $\boldsymbol{R}$ and
$\boldsymbol{R}^{\prime}$, and $P_{K}(x)$ is a Legendre polynomial. The
components $N^{K}_{cc^{\prime}}(R,R^{\prime})$ and
$V^{K}_{cc^{\prime}}(R,R^{\prime})$ are obtained from numerical integrations.
Notice that these expansions only depend on the quantum numbers of the
projectile. They do not depend on the angular momenta $J,L,L^{\prime}$, and
are therefore performed once.

The derivation of the kinetic-energy kernels is more tedious. However, we show
in the Appendix, that they can be deduced from the overlap kernel as

$\displaystyle{\cal
T}^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})=-\frac{\hbar^{2}}{2\mu}\biggl{[}\frac{\partial^{2}}{\partial
R^{2}}-\frac{L(L+1)}{R^{2}}\biggr{]}{\cal
N}^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime}).$ (25)

The calculation of the non-local kernels from (21) and (24) requires some
angular-momentum algebra. We use the addition theorem

$\displaystyle
r^{\ell}Y_{\ell}^{m}(\Omega_{r})=\sum_{\lambda=0}^{\ell}C_{\ell}^{\lambda}r_{1}^{\lambda}r_{2}^{{\ell}-{\lambda}}\bigl{[}Y_{\lambda}(\Omega_{1})\otimes
Y_{{\ell}-{\lambda}}(\Omega_{2})\bigr{]}^{\ell m},$ (26)

where $\boldsymbol{r}=\boldsymbol{r}_{1}+\boldsymbol{r}_{2}$ and with

$\displaystyle
C_{\ell}^{\lambda}=\biggl{[}\frac{4\pi(2\ell+1)!}{(2\lambda+1)!(2\ell-2\lambda+1)!}\biggr{]}^{1/2}.$
(27)

When the valence particle $v$ has a spin zero, the overlap kernel can be
expanded as

$\displaystyle{\cal N}^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})={\cal
J}(-\gamma)^{\ell+\ell^{\prime}}\sum_{K}N^{K}_{c,c^{\prime}}(R,R^{\prime})$
$\displaystyle\times\sum_{\lambda_{1}\lambda_{2}}C_{\ell}^{\lambda_{1}}C_{\ell^{\prime}}^{\lambda_{2}}\alpha^{\lambda_{1}+\lambda_{2}}R^{\ell^{\prime}+\lambda_{1}-\lambda_{2}+1}R^{\prime\ell-\lambda_{1}+\lambda_{2}+1}$
$\displaystyle\times
F^{J\pi}_{cL,c^{\prime}L^{\prime}}(K,\lambda_{1},\lambda_{2}),$ (28)

with

$\displaystyle F^{J\pi}_{cL,c^{\prime}L^{\prime}}(K,\lambda_{1},\lambda_{2})=$
$\displaystyle\biggl{\langle}\biggl{[}Y_{L}(\Omega_{R})\otimes\bigl{[}Y_{\lambda_{1}}(\Omega_{R})\otimes
Y_{\ell-\lambda_{1}}(\Omega_{R^{\prime}})\bigr{]}^{\ell}\biggr{]}^{JM}|P_{K}(\cos\theta_{R})|$
$\displaystyle\biggl{[}Y_{L^{\prime}}(\Omega_{R^{\prime}})\otimes\bigl{[}Y_{\ell^{\prime}-\lambda_{2}}(\Omega_{R})\otimes
Y_{\lambda_{2}}(\Omega_{R^{\prime}})\bigr{]}^{\ell^{\prime}}\biggr{]}^{JM}\biggr{\rangle}.$
(29)

The analytical calculation of these coefficients requires some algebra to
modify the order of angular-momentum couplings, and involves $6j$
coefficients. When the valence particle has a spin, further angular-momentum
recoupling is necessary. A simple value is obtained for
$\ell=\ell^{\prime}=0$, where we have

$\displaystyle F^{J\pi}(K,0,0)=\frac{\delta_{KJ}}{4\pi(2J+1)}.$ (30)

###  II.3 Symmetry of the non-local kernels

Let us briefly discuss the symmetry properties of the non-local potential
(17). According to Eq. (9), we have, when the cores are bosons

$\displaystyle P\Psi^{JM\pi}=\Psi^{JM\pi},$ (31)

which means that the property

$\displaystyle[H,P]=0$ (32)

should be satisfied. This implies that both $V_{Cv}$ potentials in Eqs. (1)
and (12) are identical. For example, in the ${}^{13}{\rm C}+^{12}$C system,
the $n+^{12}$C potential (real) associated with the 13C ground state should be
identical to the $n+^{12}$C optical potential which describes the neutron-
target scattering. If this condition is fulfilled, the non-local potential
$W^{J\pi}$ is symmetric, and we have

$\displaystyle
W^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})=W^{J\pi}_{c^{\prime}L^{\prime},cL}(R^{\prime},R).$
(33)

This test is very strong since, individually, all terms of the r.h.s. of (17)
are not symmetric. In practical applications, however, it may seem more
physical to choose different potentials: one which binds the $C+v$ system, and
the other which is adapted to the $C+v$ scattering. In such a case, the
symmetry property (33) is approximately satisfied (see the discussion in Ref.
Imanishi and von Oertzen (1987)). This means that some properties of the
scattering matrix, such as the symmetry or the unitarity (for real
potentials), are not any more valid. We choose to restore the symmetry of the
non-local potential through

$\displaystyle
W^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})\rightarrow\frac{1}{2}\biggl{[}W^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})+W^{J\pi}_{c^{\prime}L^{\prime},cL}(R^{\prime},R)\biggr{]}.$
(34)

We will see in some examples that these effects, in practice, are small. The
main reason is that, in general, the contribution of the optical potential
$V_{Cv}$ is small compared to the optical potential between the cores.

###  II.4 Elastic cross sections

The elastic cross sections are obtained from the scattering matrices
$\boldsymbol{U}^{J\pi}$ [see Eq. (19)]. According to the scattering theory,
the elastic cross section between different particles is defined from the
scattering amplitude as

$\displaystyle\frac{d\sigma}{d\theta}=|f(\theta)|^{2},$ $\displaystyle
f(\theta)=f^{\rm N}(\theta)+f^{\rm C}(\theta),$ (35)

where $\theta$ is the scattering angle. In this definition, $f^{\rm
C}(\theta)$ is the Coulomb amplitude and $f^{\rm N}(\theta)$ is the nuclear
amplitude, defined by

$\displaystyle f^{\rm
N}(\theta)=\frac{1}{2ik}\sum_{J}(2J+1)P_{J}(\cos\theta)(U^{J}-1)e^{2i\sigma_{J}},$
(36)

where $\sigma_{J}$ is the Coulomb phase shift. For the sake of simplicity, we
consider single-channel systems with spin 0 nuclei. The generalization is
straightforward (see, for example, Ref. Satchler (1983)).

As explained in Sec. II.A, the scattering matrices are obtained from the
resolution of a Schrödinger equation involving a non-local potential. Let us
denote as $U^{J}_{0}$ and $g^{J}_{0}(R)$ the scattering matrix and wave
function obtained from the local term (16) only. These quantities are obtained
without symmetrization of the wave function (9). The integral definition of
the scattering matrix Canto and Hussein (2013) provides a relationship between
$U^{J}$ and $U^{J}_{0}$ as

$\displaystyle U^{J}=U^{J}_{0}+U^{J}_{\rm ex},$ $\displaystyle U^{J}_{\rm
ex}=-\frac{i}{\hbar}\iint
g^{J}_{0}(R)W^{J}(R,R^{\prime})g^{J}(R^{\prime})dRdR^{\prime}.$ (37)

Consequently the nuclear scattering amplitude can be written as

$\displaystyle f^{\rm N}(\theta)=f^{\rm N}_{0}(\theta)+f^{\rm N}_{\rm
ex}(\theta),$ (38)

where $f^{\rm N}_{0}(\theta)$ is obtained from (36) with the scattering
matrices $U^{J}_{0}$, and where the exchange amplitude $f^{\rm N}_{\rm
ex}(\theta)$ is the non-local contribution

$\displaystyle f^{\rm N}_{\rm
ex}(\theta)=\frac{1}{2ik}\sum_{J}(2J+1)P_{J}(\cos\theta)U^{J}_{\rm
ex}e^{2i\sigma_{J}}.$ (39)

A similar decomposition has been suggested in Refs. Phuc  _et al._ (2018,
2021, 2019). Since the calculation of the exchange amplitude $f^{\rm N}_{\rm
ex}(\theta)$ is based on a non-local potential, it is common in the literature
to assume that it can be simulated by the transfer of the valence particle
(see for example Ref. Phuc  _et al._ (2019) and references therein). In other
words, the exchange amplitude can be approximated as

$\displaystyle f^{\rm N}_{\rm ex}(\theta)\approx f_{\rm tr}(\pi-\theta),$ (40)

where $f_{\rm tr}(\pi-\theta)$ is the transfer amplitude, usually computed at
the DWBA approximation Thompson (1988). In this approximation, however, the
exact wave function $g^{J}(R)$ is replaced by the non-symmetrized wave
function $g^{J}_{0}(R)$ in Eq. (37), and the non-local potential
$W^{J}(R,R^{\prime})$ is approximated from an auxiliary potential.

Another difference is related to the spectroscopic factors. As expected in
transfer calculations, the DWBA amplitude $f_{\rm tr}(\pi-\theta)$ is
multiplied by the spectroscopic factor of the projectile. This multiplicative
factor, however, does not show up in the present approach. The two-body wave
function (6) is of course an approximation which can be improved, either by
introducing a spectroscopic factor or by including core excitations. The
introduction of a spectroscopic factor in (6), however, represents a global
renormalization of the expansion (7), and therefore of the symmetrized
definition (9). The scattering matrices deduced from (19) are therefore not
affected by a spectroscopic factor. Introducing core excitations in the
present approach is a challenge for future works, but is beyond the scope of
the present work.

###  II.5 Lagrange functions

As mentioned before, the scattering matrices are calculated with the
$R$-matrix method, which is based on a channel radius $a$ and on the choice of
basis functions $u_{i}$. As in previous works, we choose Lagrange functions
which permit fast and accurate calculations of the matrix elements, in
particular for non-local potentials (see Ref. Baye (2015) for detail).

The calculation of the $R$ matrix is based on matrix elements between basis
functions over the internal region. The main input is the matrix defined from

$\displaystyle C^{J\pi}_{cLi,c^{\prime}L^{\prime}j}=$ $\displaystyle\langle
u_{i}|(T_{R}+E_{c}-E)\delta_{cc^{\prime}}\delta_{LL^{\prime}}$
$\displaystyle+V^{J\pi}_{cL,c^{\prime}L^{\prime}}+W^{J\pi}_{cL,c^{\prime}L^{\prime}}|u_{j}\rangle.$
(41)

For example, matrix elements of the local and non-local potentials are given
by

$\displaystyle\langle
u_{i}|V|u_{j}\rangle=\int_{0}^{a}u_{i}(R)V(R)u_{j}(R)\,dR,$
$\displaystyle\langle
u_{i}|W|u_{j}\rangle=\int_{0}^{a}\int_{0}^{a}u_{i}(R)W(R,R^{\prime})u_{j}(R^{\prime})\,dRdR^{\prime}.$
(42)

These calculations are greatly simplified by using Lagrange functions for
$u_{i}$ which are defined by

$\displaystyle
u_{i}(R)=(-1)^{N+i}\frac{R}{R_{i}}\sqrt{R_{i}\biggl{(}1-\frac{R_{i}}{a}\biggr{)}}\,\frac{P_{N}(2R/a-1)}{R-R_{i}},$
(43)

where $R_{i}$ are the zeros of

$\displaystyle P_{N}(2R/a-1)=0.$ (44)

The normalization of (43) is chosen in such a way that the Lagrange condition

$\displaystyle u_{i}(R_{j})=\frac{1}{\sqrt{a\lambda_{i}}}\delta_{ij}$ (45)

is satisfied. In this equation, $\lambda_{i}$ is the weight of the Gauss-
Legendre quadrature associated with the $[0,1]$ interval.

With the choice of basis function (43), the calculation of the matrix elements
(42) is extremely simple if the Gauss approximation of order $N$ is used for
the quadratures. The matrix elements are given by

$\displaystyle\langle u_{i}|V|u_{j}\rangle=V(R_{i})\delta_{ij},$
$\displaystyle\langle
u_{i}|W|u_{j}\rangle=a\sqrt{\lambda_{i}\lambda_{j}}W(R_{i},R_{j}),$ (46)

and no numerical integral is required for the matrix elements. We refer to
Refs. Descouvemont and Baye (2010); Baye (2015) for details.

##  III Applications

###  III.1 The ${}^{13}{\rm C}+^{12}$C system

The ${}^{13}{\rm C}+^{12}$C system has been intensively studied experimentally
Liénard  _et al._ (1995); Voit  _et al._ (1988) as well as theoretically von
Oertzen (1970); Imanishi and von Oertzen (1987); Imanishi  _et al._ (1997). It
is known that non-local effects can be simulated by a parity-dependent optical
potential Buttle and Goldfarb (1966). Owing to the one-neutron exchange, the
potentials for even and odd partial waves are different Baye (1986). In the
“extreme” situation of identical nuclei, odd partial waves are strictly
forbidden.

We have determined the non-local potential (17) from ${}^{12}{\rm C}+^{12}$C
and $n+^{12}$C potentials. For ${}^{12}{\rm C}+^{12}$C, we take the optical
potential derived by Treu et al. Treu  _et al._ (1980), and defined (in MeV)
as

$\displaystyle V_{\rm CC}(r)=$
$\displaystyle-\frac{100}{1+\exp[(r-5.45)/0.48]}$
$\displaystyle-i\frac{15}{1+\exp[(r-5.77)/0.26]}$ (47)

where $r$ is expressed in fm. A Coulomb point-sphere potential of radius
$R_{C}=5.45$ fm is added. This optical potential is fitted on elastic-
scattering data around the Coulomb barrier. In all applications, we use the
integer masses with $\hbar^{2}/2m_{N}=20.736$ MeV fm2 ($m_{N}$ is the nucleon
mass).

The $n+^{12}$C potential is chosen as in Ref. von Oertzen and Imanishi (1984),
i.e.

$\displaystyle V_{nC}(\boldsymbol{r})=$
$\displaystyle-\frac{V_{0}}{1+\exp[(r-r_{0})/a_{0}]}$
$\displaystyle-(\boldsymbol{\ell}\cdot\boldsymbol{s})\frac{V_{\ell
s}}{r}\frac{d}{dr}\frac{1}{1+\exp[(r-r_{0})/a_{0}]},$ (48)

where $V_{0}=62.70$ MeV for $\ell$ even and $50.59$ MeV for $\ell$ odd, and
where $r_{0}=2.656$ fm and $a_{0}=0.705$ fm. The spin-orbit amplitude is
$V_{\ell s}=28.406$ MeV. This potential reproduces the experimental energies
of the first $1/2^{-},1/2^{+}$, and $5/2^{+}$ states in 13C. Between the
target and the projectile, the same $n+^{12}$C potential is adopted for each
partial wave, corresponding to the central part of (48). Once potentials (47)
and (48) are determined, the model does not contain any free parameter.

In Fig. 2, the ${}^{13}{\rm C}+^{12}$C elastic cross sections at
$\mbox{$E_{\rm c.m.}$}=7.8$ and 14.2 MeV are shown. In each case, we consider
four conditions: (1) when only the local potential is included, the backward
angle enhancement of the cross section is not reproduced; (2) the non-local
calculation involving the 13C ground state only (dashed line) reproduces
fairly well the data; (3) when the $1/2^{+}$ and $5/2^{+}$ excited states are
introduced (solid line), elastic scattering is not significantly modified; (4)
in the $n+^{12}$C potential between the target and the projectile (dotted
line), we have replaced the real potential (48) by the Koning-Delaroche
parametrization Koning and Delaroche (2003). Although some symmetry properties
are lost (see Sec. II.B), there is a weak influence on the ${}^{13}{\rm
C}+^{12}$C cross section.

Figure 2: ${}^{13}{\rm C}+^{12}$C elastic cross sections (divided by the
Rutherford cross section) at $\mbox{$E_{\rm c.m.}$}=7.8$ MeV (a) and 14.2 MeV
(b). The data are taken from Refs. Liénard  _et al._ (1995) (full dots) and
Voit  _et al._ (1988) (open dots). The solid lines are obtained with the
$1/2^{-},1/2^{+},5/2^{+}$ states of 13C. The dashed (red) lines correspond to
the $1/2^{-}$ ground state only. The dotted (blue) lines are obtained with
different $n+^{12}$C potentials in the entrance and exit channels (see text).
Panel (c) presents the amplitudes $|U^{J\pi}_{L}|$ of the scattering matrices
at $\mbox{$E_{\rm c.m.}$}=7.8$ MeV, and for $J=L-1/2$.

Figure 2(c) presents the amplitudes $|U^{J\pi}_{L}|$ of the scattering
matrices (elastic channel) at $\mbox{$E_{\rm c.m.}$}=7.8$ MeV. We choose here
$J=L-1/2$ but a similar behaviour is observed for $J=L+1/2$. Without the non-
local part of the potential, the variation is smooth. As expected, the non-
locality leads to a splitting between odd and even $L$ values. This property
gives rise to the backward angle enhancement of the cross section, and
justifies the use of parity-dependent optical potentials Buttle and Goldfarb
(1966) to simulate non-local effects. An advantage of the present method is
that it does not require any additional parameter. In addition, excited states
of the $n+^{12}$C system are included in a straightforward way.

In Fig. 3, we investigate the inelastic cross sections to the 13C($1/2^{+}$)
and 13C($5/2^{+}$) states, which have been measured in Ref. Fröhlich  _et al._
(1984) at energies around the Coulomb barrier. The effect of non-locality is
quite important. With the local potential only, the theoretical cross sections
are far below the data. The calculation involves the ${}^{13}{\rm
C(gs},1/2^{+},5/2^{+})+^{12}$C channels. Here the cross sections are more
sensitive to the choice of the $n+^{12}$C potential. Let us emphasize that
there is no fit of the cross sections. All inputs are kept identical as for
elastic scattering.

Figure 3: Inelastic ${}^{13}{\rm C}+^{12}$C cross sections to the $1/2^{+}$
(a) and $5/2^{+}$ (b) states of 13C at $\mbox{$E_{\rm c.m.}$}=9.88$ MeV. The
experimental data are taken from Ref. Fröhlich  _et al._ (1984). Solid
(dashed) lines are obtained with identical (different) $n+^{12}$C potentials
in the entrance and exit channels (see text).

###  III.2 The ${}^{13}{\rm N}+^{12}$C system

With the development of radioactive beams, the ${}^{13}{\rm N}+^{12}$C mirror
system has attracted much attention in the literature Liénard  _et al._
(1995); Imanishi  _et al._ (1997). Measurements have provided some information
about charge-symmetry and about the parity effect Liénard  _et al._ (1995).
Figure 4 shows the ${}^{13}{\rm N}+^{12}$C calculated cross sections. Only the
13N ground state is bound, and has been introduced in the calculation. With
respect to the ${}^{13}{\rm C}+^{12}$C system, the only difference is the
introduction of a Coulomb term for $p+^{12}$C (with $R_{C}=2.7$ fm). The
binding energy of 13N is $-1.90$ MeV, in fair agreement with experiment
($-1.94$ MeV). Here again the model reproduces remarkably well the
experimental data Liénard  _et al._ (1995).

Figure 4: ${}^{13}{\rm N}+^{12}$C elastic cross sections (divided by the
Rutherford cross section) at $\mbox{$E_{\rm c.m.}$}=7.8$ (a) and 9.6 (b) MeV.
The data are taken from Ref. Liénard  _et al._ (1995). Solid lines correspond
to the full calculation, and dashed lines to the local potential only.

###  III.3 The ${}^{16}{\rm O}+^{12}$C system

In the ${}^{16}{\rm O}+^{12}$C scattering, an $\alpha$ particle is exchanged
between the target and the projectile. This system has been intensively
investigated in the literature (see, for example, Refs. Phuc  _et al._ (2019,
2018, 2021) for recent works).

We consider the elastic-scattering data of Villari et al. Villari  _et al._
(1989) at the typical energy $\mbox{$E_{\rm c.m.}$}=23.14$ MeV, where a
backward-angle enhancement of the cross sections is observed. We use the same
${}^{12}{\rm C}+^{12}$C core-core potential (47) as in previous applications.
For the $\alpha+^{12}$C system, we adopt the potentials used in Ref. Oulebsir
_et al._ (2012), i.e. a Woods-Saxon potential with a range $R_{0}=4.15$ fm and
a diffuseness $a=0.55$ fm. In addition to the $0^{+}$ ground state, we include
the $1^{-},3^{-}$, and $2^{+}$ excited states. Since the $2^{+}$ state
presents a cluster structure, $R_{0}$ is chosen larger (4.5 fm) for this
state. The Coulomb potential has a point-sphere shape with a radius
$R_{C}=4.15$ fm. The depths of the potentials $V_{0}$ are adjusted to the
experimental binding energies, which provides $V_{0}=-43.25,-68.95,-41.3$, and
$-42.4$ MeV for the $0^{+},2^{+},1^{-}$, and $3^{-}$ states, respectively.
Between the target and the projectile, the same $\alpha+^{12}$C potential is
adopted, corresponding to the 16O ground state.

The elastic cross section is presented in Fig. 5(a) in different conditions.
The calculation with the local potential provides a fair description of the
data up to $\theta\approx 90^{\circ}$, but does not reproduce the enhancement
at large angles. With the non-local term, even if the oscillations at backward
angles are not exactly reproduced, the role of inelastic channels is obvious.
The present model, based on the exchange of an $\alpha$ particle during the
collision, cannot be expected to be perfect. Other channels are open, such as
the neutron or proton transfer, but are neglected here. There are usually
treated by phenomenological potentials involving additional parameters.

In the ${}^{16}{\rm O}+^{12}$C system, the role of the core-valence potential
in more important than in ${}^{13}{\rm C}+^{12}$C. To assess this sensitivity,
we have also used the $\alpha$ optical potential of Avrigeanu et al. Avrigeanu
_et al._ (2009) (dotted line in Fig. 5(a)). This is explained by the long
range of the $\alpha+^{12}$C potential, due to the Coulomb term. In this
system, using a consistent $\alpha+^{12}$C potential in the entrance and in
the exit channels seems more appropriate.

The amplitude of the scattering matrices $|U^{J\pi}|$ is displayed in Fig.
5(b). As for ${}^{13}{\rm C}+^{12}$C, the calculation with the non-local term
provides differences between even and odd partial waves. These differences,
however, are weaker than in ${}^{13}{\rm C}+^{12}$C. As observed in Ref. Phuc
_et al._ (2019), the even-odd effect is stronger around the grazing angular
momentum ($J\approx 18$).

Figure 5: ${}^{16}{\rm O}+^{12}$C elastic cross section (divided by the
Rutherford cross section) at $\mbox{$E_{\rm c.m.}$}=23.14$ MeV (a). The data
are taken from Ref. Villari  _et al._ (1989). The solid lines correspond to
the multichannel model. The dashed lines are obtained with the 16O ground
state only, and the dotted lines with different $\alpha+^{12}$C potentials
(see text). Panel (b) presents the amplitudes $|U^{J\pi}|$ of the scattering
matrices.

##  IV Discussion of the non-locality

###  IV.1 Local equivalent potentials

The effects of the non-locality can be simulated by an equivalent local
potential. For the sake of simplicity, we assume a single-channel problem. The
extension to multichannel systems is simple and does not modify the
conclusions. In a single-channel model, the radial Schrödinger equation reads

$\displaystyle\biggl{[}-\frac{\hbar^{2}}{2\mu}\biggl{(}\frac{d^{2}}{dr^{2}}-\frac{L(L+1)}{R^{2}}\biggr{)}+V^{J\pi}(R)-E\biggr{]}g^{J\pi}(R)$
$\displaystyle\hskip 28.45274pt+\int
W^{J\pi}(R,R^{\prime})g^{J\pi}(R^{\prime})dR^{\prime}=0,$ (49)

and can be replaced by

$\displaystyle\biggl{[}-\frac{\hbar^{2}}{2\mu}\biggl{(}\frac{d^{2}}{dr^{2}}-\frac{L(L+1)}{R^{2}}\biggr{)}+V^{J\pi}_{\rm
eq}(R)-E\biggr{]}g^{J\pi}(R)=0,$ (50)

where the equivalent potential $V^{J\pi}_{\rm eq}(R)$ is given by

$\displaystyle V^{J\pi}_{\rm eq}(R)=V^{J\pi}(R)+\frac{1}{g^{J\pi}(R)}\int
W^{J\pi}(R,R^{\prime})g^{J\pi}(R^{\prime})dR^{\prime}.$ (51)

This potential depends on the angular momentum, and presents singularities at
the nodes of the wave functions. Thompson et al. Thompson  _et al._ (1989)
have proposed to define a smooth, $J$-independent, effective potential by

$\displaystyle V^{\pi}_{\rm eff}(R)=\frac{\sum_{J}\omega^{J\pi}(R)V_{\rm
eq}^{J\pi}(R)}{\sum_{J}\omega^{J\pi}(R)},$ (52)

where the weight factors $\omega^{J\pi}(R)$ are given by

$\displaystyle\omega^{J\pi}(R)=(2J+1)(1-|U^{J\pi}|^{2})|g^{J\pi}(R)|^{2}.$
(53)

In this way, the influence of the nodes is reduced and the potential (53) does
not depend on $J$. However, the scattering matrices obtained with (52) are not
strictly identical to those obtained with (49) or (50). A test with the cross
sections must be performed to check the accuracy of the potential (52).

As we expect the equivalent local potentials to depend on parity, the
potential (52) is defined for each parity. From $V^{+}_{\rm eff}(R)$ and
$V^{-}_{\rm eff}(R)$, we determine central and parity-dependent potentials as

$\displaystyle V_{0}(R)=\frac{1}{2}\bigl{(}V^{+}_{\rm eff}(R)+V^{-}_{\rm
eff}(R)\bigr{)},$ $\displaystyle V_{\pi}(R)=\frac{1}{2}\bigl{(}V^{+}_{\rm
eff}(R)-V^{-}_{\rm eff}(R)\bigr{)}.$ (54)

The present work, based on rigorous non-local potentials, offers the
possibility to investigate the parity potential which, in general, is
phenomenological (see, for example, Ref. Liénard  _et al._ (1995)). Notice
that both $V_{0}(R)$ an $V_{\pi}(R)$ contain real and imaginary components. A
parity effect can be also deduced from techniques based on data inversion
Mackintosh (2019); Phuc  _et al._ (2019).

###  IV.2 Asymptotic form of the non-local kernels

Here we present some qualitative aspects regarding the non-local kernels, and
in particular the overlap kernel ${\cal
N}^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})$. Let us first discuss the
overlap functions $N^{K}_{cc^{\prime}}(R,R^{\prime})$ which show up in Eq.
(24). To simplify the presentation, we assume that the valence particle is a
neutron in a $s$ state. In that case, the core-valence wave function tends to

$\displaystyle u_{0}(r)\rightarrow C\exp(-k_{B}r),$ (55)

where $k_{B}$ is the wave number and $C$ the asymptotic normalization
coefficient (ANC).

To develop further, we use the expansion Buttle and Goldfarb (1966)

$\displaystyle\frac{\exp(-k|\boldsymbol{r}_{1}-\boldsymbol{r}_{2}|)}{k|\boldsymbol{r}_{1}-\boldsymbol{r}_{2}|}=$
$\displaystyle\hskip
28.45274pt\frac{2}{\pi}\sum_{\ell}(2\ell+1)i_{\ell}(kr_{<})k_{\ell}(kr_{>})P_{\ell}(\cos\theta_{r}),$
(56)

where $r_{<}=\min(r_{1},r_{2})$ and $r_{>}=\max(r_{1},r_{2})$, and where
$\theta_{r}$ is the angle between $\boldsymbol{r}_{1}$ and
$\boldsymbol{r}_{2}$. In this definition, $i_{\ell}(x)$ and $k_{\ell}(x)$ are
modified spherical Bessel functions Abramowitz and Stegun (1972). For large
arguments, they tend to

$\displaystyle i_{\ell}(x)\rightarrow\frac{\exp(x)}{2x},$ $\displaystyle
k_{\ell}(x)\rightarrow\pi\frac{\exp(-x)}{2x}.$ (57)

Using the expansion (56) for $u_{0}(r)$ and $u_{0}(r^{\prime})$, and using
relations (23), we find, for large $(R,R^{\prime})$ values

$\displaystyle\frac{u_{0}(r)u_{0}(r^{\prime})}{rr^{\prime}}\rightarrow\sum_{K}N^{K,{\rm
as}}(R,R^{\prime})P_{K}(\cos\theta_{R}).$ (58)

In the range $\alpha R^{\prime}<R<R^{\prime}/\alpha$, the asymptotic kernels
are defined by

$\displaystyle N^{K,{\rm
as}}(R,R^{\prime})=\frac{4}{\pi^{2}}k_{B}^{2}C^{2}(-1)^{K}\sum_{\ell_{1}\ell_{2}}(2\ell_{1}+1)(2\ell_{2}+1)$
$\displaystyle\times\langle\ell_{1}0\ell_{2}0|K0\rangle^{2}k_{\ell_{2}}(\gamma
k_{B}R)k_{\ell_{1}}(\gamma k_{B}R^{\prime})$ $\displaystyle\times
i_{\ell_{1}}(\alpha\gamma k_{B}R)i_{\ell_{2}}(\alpha\gamma k_{B}R^{\prime}),$
(59)

and Eq. (57) provides

$\displaystyle N^{K,{\rm as}}(R,R^{\prime})\sim\frac{1}{R^{2}R^{\prime
2}}\exp\biggl{(}-\frac{k_{B}}{1+\alpha}(R+R^{\prime})\biggr{)}.$ (60)

As $\alpha$ is in general close to 1, Eqs. (59,60) are valid for $R\approx
R^{\prime}$. Otherwise, a similar development gives

$\displaystyle N^{K,{\rm as}}(R,R^{\prime})\sim\frac{1}{R^{2}R^{\prime
2}}\exp\biggl{(}-\frac{k_{B}}{1-\alpha}|R-R^{\prime}|\biggr{)}.$ (61)

This shows that the non-locality overlap kernel presents an exponential
decrease, associated with the binding energy of the projectile.

When the angular momentum is not an $s$ wave, the expansion (56) can be
generalized Buttle and Goldfarb (1966), but this does not change the general
trend. In addition, if the transferred particle is charged, the asymptotic
behaviour takes the form

$\displaystyle u_{0}(r)\rightarrow C\frac{\exp(-k_{B}r)}{r^{\eta_{B}}},$ (62)

where $\eta_{B}$ is the Sommerfeld parameter. The faster decrease can be
simulated by using (55) with a (larger) effective wave number, which simulates
Coulomb effects Buttle and Goldfarb (1968). Consequently, Eq. (59) remains
qualitatively valid, even for charged transferred particles.

The non-local potentials (17) also involve kinetic-energy and nucleus-nucleus
potential terms. As discussed in Sec. II.B, the kinetic-energy kernel is
directly deduced from a second derivative of the overlap. The asymptotic
behaviour is therefore similar to (59). To determine the potential
contribution, we start from definition (24). The first term in the potential
$V_{Cv}(r^{\prime})$ is involved in standard DWBA calculations (post form)
Satchler (1983). For a transferred neutron, the nuclear contribution in
$V_{Cv}$ makes this term short-ranged. The asymptotic behaviour (60) is
therefore modified, with a smaller range. In contrast, the core-core
interaction $V_{CC}$ always presents a Coulomb term. One therefore expects
this interaction to be dominant at large distances.

The non-local potentials in Eq. (17) are obtained from (24),(55), and angular
matrix elements Satchler (1983); Shubhchintak and Descouvemont (2019). The
procedure is similar to the one followed in DWBA calculations. If the internal
angular momenta are taken as $\ell=\ell^{\prime}=0$, the calculation is
simple, as the sum over $K$ contains a single term $K=J$. We have

$\displaystyle{\cal N}^{J\pi}(R,R^{\prime})=$ $\displaystyle{\cal
J}RR^{\prime}N^{J}(R,R^{\prime})/(2J+1),$ $\displaystyle{\cal
V}^{J\pi}(R,R^{\prime})=$ $\displaystyle{\cal
J}RR^{\prime}\bigl{(}V^{J}_{Cv}(R,R^{\prime})+V^{J}_{CC}(R,R^{\prime})\bigr{)}$
$\displaystyle/(2J+1).$ (63)

Equation (59) contains a phase factor $(-1)^{J}$ at large distances, which
shows that the non-local potentials have opposite signs for even and odd
partial waves. This effect arises from the symmetrization of the wave
functions for the core exchange. As shown by von Oertzen von Oertzen (1970),
it can be simulated by a local, parity-dependent, potential. An application to
the analysis of the ${}^{13}{\rm C}+^{12}$C and of the ${}^{13}{\rm N}+^{12}$C
systems can be found in Ref. Liénard  _et al._ (1995).

###  IV.3 Application to ${}^{13}{\rm C}+^{12}$C

In Fig. 6, we present the effective potentials for the ${}^{13}{\rm C}+^{12}$C
system at $\mbox{$E_{\rm c.m.}$}=7.8$ MeV. For the sake of simplicity, we
discuss the single-channel calculation. Figure 6(a) shows the local potential
(dashed lines), and the parity-dependent effective potentials (solid lines)
obtained with (54). The elastic cross sections obtained with this effective
potential is very close to the original cross section (they are almost
indistinguishable in a figure). As expected, the potentials present small
oscillations due to the nodes in the scattering wave functions. We have
repeated the calculations with different numerical conditions (channel radius,
number of basis functions), and checked that the effective potentials are
quite stable. Figure 6(b) display the parity potential $V_{\pi}$ [see Eq.
(54)]. The parity effect is important in the real part as well as in the
imaginary part.

Figure 6: Local effective potentials for the ${}^{13}{\rm C}+^{12}$C system at
$\mbox{$E_{\rm c.m.}$}=7.8$ MeV. In panel (a), the dashed curve corresponds to
the local potential, and the solid curves to the local equivalent potentials
(52). Panel (b) displays the parity potential $V_{\pi}(R)$ (54).

In Fig. 7, we display the different contributions in the non-local kernel
$W^{J\pi}(R,R^{\prime})$ (17) for $R=R^{\prime}$ at $\mbox{$E_{\rm
c.m.}$}=7.8$ MeV. Four terms are present: the overlap, the kinetic energy, and
two contributions of the potential (the core-core and core-valence terms).
Partial waves $J=1/2^{+}\ (L=1)$ and $J=1/2^{-}\ (L=0)$ are shown in panels
(a) and (b). The change of sign between both parities is confirmed. At short
distances, the main contribution comes from $V_{CC}$, but the kinetic-energy
is dominant at large distances. For consistency, the overlap is multiplied by
$-\mbox{$E_{\rm c.m.}$}$, and represents a small contribution. The inset of
panel (a) confirms the asymptotic behaviour. All terms but $V_{Cv}$ have the
same exponential behaviour. The contribution associated with $V_{Cv}$ presents
a faster decrease owing to the absence of the Coulomb interaction.

Figure 7: Contributions of the overlap, kinetic energy and potentials
($V_{CC}$ and $V_{Cn}$) to the non-local kernel (17) for $R=R^{\prime}$, and
for $J=1/2^{+}$ (a) and $J=1/2^{-}$ (b) in the ${}^{13}{\rm C}+^{12}$C system
($\mbox{$E_{\rm c.m.}$}=7.8$ MeV). The overlap kernel is multiplied by
$-\mbox{$E_{\rm c.m.}$}$. For $V_{CC}$, the real part is displayed. The inset
focuses on the long-range part in a logarithmic scale.

In Fig. 8, we analyze the non-locality for $R\neq R^{\prime}$ and for
$J=1/2^{+}$. We plot the overlap and real-potential kernels as a function of
$(R^{\prime}-R)$ for various $R$ values. For small $R$ values, the shape of
the overlap kernel is close to a Gaussian, which justifies the Perey-Buck
approximation Perey and Buck (1962). However, for $R>1$ fm, the shape is more
complicated than a Gaussian factor. As expected the main effect of the non-
locality comes from $|R-R^{\prime}|\lesssim 1$ fm.

Figure 8: Non-local kernels ${\cal N}(R,R^{\prime})$ (a) and $W(R,R^{\prime})$
(real part, (b)) for the ${}^{13}{\rm C}+^{12}$C system ($J=1/2^{+}$), as a
function of $(R^{\prime}-R)$ [see Eq. (17)]. The curves are plotted for
different $R$ values.

###  IV.4 Application to ${}^{16}{\rm O}+^{12}$C

The local effective potentials (52) for ${}^{16}{\rm O}+^{12}$C at
$\mbox{$E_{\rm c.m.}$}=23.14$ MeV are presented in Fig. 9(a). Again, we limit
the discussion to the single-channel system. Since the depth of these
potentials is rather large, the differences with the exact local potential
(16) are small, indistinguishable at the scale of the figure. Figure 9(b)
shows the parity potential (54). As for ${}^{13}{\rm C}+^{12}$C, these
potentials present oscillations due to the nodes of the scattering wave
functions. Although the amplitude is small, the effect of the parity potential
extends to large distances. This explains that this potential provides a large
backward-angle enhancement of the cross section (see Fig. 5). The cross
sections provided by the effective potential (52) are very close to the
original cross sections of Fig. 5.

Figure 9: Local effective potentials for the ${}^{16}{\rm O}+^{12}$C system at
$\mbox{$E_{\rm c.m.}$}=23.14$ MeV. In panel (a), the local potential and the
equivalent potential (52) are superimposed. Panel (b) displays the parity
potential $V_{\pi}(R)$ (54).

The decomposition of the non-local kernel (17) in different terms is presented
in Fig. 10 for $J=0^{+}$ (a) and $J=1^{-}$ (b). The dominant contribution
comes from the core-core potential $V_{CC}$. There is a cancellation effect of
the overlap, kinetic-energy and $V_{C\alpha}$ contributions. At large
distances the decrease of the potentials (inset of Fig. 10(a) is much faster
than in ${}^{13}{\rm C}+^{12}$C (see Fig. 7). This is explained by Eq. (60)
since $k_{B}$ is much larger in 16O than in 13C (and coefficient $\alpha$ is
smaller).

Figure 10: Contributions of the overlap, kinetic energy and potential
($V_{CC}$ and $V_{C\alpha}$) to the non-local kernel (17) for $R=R^{\prime}$,
and for $J=0^{+}$ (a) and $J=1^{-}$ (b) in the ${}^{16}{\rm O}+^{12}$C system.
The real part of $V_{CC}$ is shown. The inset focuses on the long-range part
in a logarithmic scale.

The non-locality for $R\neq R^{\prime}$ is illustrated in Fig. 11 for
$J=0^{+}$. As for ${}^{13}{\rm C}+^{12}$C, the shape of the overlap kernel is
close to a Gaussian for small $R$ values, but deviates when $R$ increases.

Figure 11: Non-local kernels ${\cal N}(R,R^{\prime})$ (a) and
$W(R,R^{\prime})$ (real part, (b)) for the ${}^{16}{\rm O}+^{12}$C system
($J=0^{+}$), as a function of $(R^{\prime}-R)$ [see Eq. (17)]. The curves are
plotted for different $R$ values.

##  V Conclusion

Our main goal is an exploratory study of a rigorous method to treat exchange
effects in nucleus-nucleus scattering. Starting from a three-body model, we
have derived local and non-local kernels in a coupled-channel formalism. This
represents a natural extension of the three-body CDCC method. The coupled-
channel system, involving non-local potentials, is solved with the $R$-matrix
theory, associated with the Lagrange-mesh technique. This permits fast and
accurate calculations of the scattering matrices and of the cross sections.

We have applied the formalism to typical reactions: ${}^{13}{\rm C}+^{12}$C,
${}^{13}{\rm N}+^{12}$C and ${}^{16}{\rm O}+^{12}$C, which illustrate the
transfer of a neutron, of a proton, and of an $\alpha$ particle, respectively.
In each case we started from a core-core optical potential, taken from the
literature, and which fits ${}^{12}{\rm C}+^{12}$C elastic-scattering data. As
expected, exchange effects lead to a backward-angle enhancement of the elastic
cross sections. For the core-valence potential, we have two “natural” choices:
either it is fitted on the spectroscopic properties of the heavy particle, or
it is fitted on scattering properties. The first choice guarantees the
symmetry of the scattering matrix, since the same Hamiltonian is used in the
entrance and exit channels. With our examples, we have however shown that both
options provide similar cross sections. This can be explained by the weak
contribution of this potential to the non-local kernel.

We have shown in Sec. IV that the bound-state wave functions of the $C+v$
system need to be accurately determined up to large distances. If not, the
calculation of the various kernels is inaccurate, and the scattering matrices
are unstable. For this reason, CDCC calculations involving continuum states
raise numerical difficulties with pseudostates as they tend to zero at very
large distances only. The use of bins would be preferable. In our
applications, however, the binding energies of the projectile are large, and
no strong continuum effects are expected.

The present work opens the path to more ambitious calculations, such as the
$\alpha+^{8}$Be system, where exchange effects could affect the existing
calculations Ogata  _et al._ (2009). Also it could be extended to more
complicated systems, such as $d+^{11}$Be or ${}^{13}{\rm N}+^{13}$C which
require a four-body theory. Another application of the formalism deals with
Coupled Reaction Channel (CRC) calculations where similar non-local kernels
are involved Satchler (1983).

## Acknowledgments

This work was supported by the Fonds de la Recherche Scientifique - FNRS under
Grant Numbers 4.45.10.08 and J.0049.19. It benefited from computational
resources made available on the Tier-1 supercomputer of the Fédération
Wallonie-Bruxelles, infrastructure funded by the Walloon Region under the
grant agreement No. 1117545.

##  Appendix A Calculation of the non-local kinetic energy

The kinetic-energy kernel ${\cal T}^{J\pi}_{cL,c^{\prime}L^{\prime}}$ is
implicitly defined by

$\displaystyle\int{\cal
T}^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})g^{J\pi}_{c^{\prime}L^{\prime}}(R^{\prime})dR^{\prime}=$
$\displaystyle
R\langle\varphi^{JM\pi}_{cL}|T_{\boldsymbol{R}}|\varphi^{JM\pi}_{c^{\prime}L^{\prime}}\dfrac{g^{J\pi}_{c^{\prime}L^{\prime}}(R^{\prime})}{R^{\prime}}\rangle,$
(64)

where the integration is performed over $\boldsymbol{r}$ and $\Omega_{R}$. An
explicit expression for this kernel can be obtained by, first, writing the
kinetic-energy operator $T_{\boldsymbol{R}}$ in spherical coordinates,

$\displaystyle T_{\boldsymbol{R}}=-\dfrac{\hbar^{2}}{2\mu
R}\dfrac{d^{2}}{dR^{2}}R+\dfrac{\hat{L}^{2}_{\boldsymbol{R}}}{2\mu R^{2}},$
(65)

where $\hat{L}_{\boldsymbol{R}}$ is the orbital angular momentum associated
with the coordinate $\boldsymbol{R}$. The radial part of $T_{\boldsymbol{R}}$
can be moved outside the matrix element since the function
$\varphi^{JM\pi}_{cL}$ is independant of $R$ and since the matrix element does
not involve any integration over the coordinate $R$. Besides, since the
operator $\hat{L}^{2}_{\boldsymbol{R}}$ is Hermitian, it can be applied, with
a trivial effect, on the bra. We therefore obtain

$\displaystyle
R\langle\varphi^{JM\pi}_{cL}|T_{\boldsymbol{R}}|\varphi^{JM\pi}_{c^{\prime}L^{\prime}}\dfrac{g^{J\pi}_{c^{\prime}L^{\prime}}(R^{\prime})}{R^{\prime}}\rangle$
$\displaystyle=$
$\displaystyle-\dfrac{\hbar^{2}}{2\mu}\left[\dfrac{d^{2}}{dR^{2}}-\dfrac{L(L+1)}{R^{2}}\right]R\langle\varphi^{JM\pi}_{cL}|\varphi^{JM\pi}_{c^{\prime}L^{\prime}}\dfrac{g^{J\pi}_{c^{\prime}L^{\prime}}(R^{\prime})}{R^{\prime}}\rangle$
(66) $\displaystyle=$
$\displaystyle\int\left\\{-\dfrac{\hbar^{2}}{2\mu}\left[\dfrac{\partial^{2}}{\partial
R^{2}}-\dfrac{L(L+1)}{R^{2}}\right]{\cal
N}^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})\right\\}g^{J\pi}_{c^{\prime}L^{\prime}}(R^{\prime})dR^{\prime},$
(67)

where the definition (18) of the norm kernel and the Leibniz integral rule
have been used. By comparison of (64) and (67), we get the relation (25),
linking the kinetic-energy and norm kernels,

$\displaystyle{\cal
T}^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime})=-\frac{\hbar^{2}}{2\mu}\biggl{[}\frac{\partial^{2}}{\partial
R^{2}}-\frac{L(L+1)}{R^{2}}\biggr{]}{\cal
N}^{J\pi}_{cL,c^{\prime}L^{\prime}}(R,R^{\prime}).$ (68)

## References

  * Feshbach (1958) H. Feshbach, Ann. Rev. Nucl. Sci. 8, 49 (1958).
  * Dickhoff and Charity (2019) W. H. Dickhoff and R. J. Charity, Prog. Part. Nucl. Phys. 105, 252 (2019).
  * Rawitscher (1974) G. H. Rawitscher, Phys. Rev. C 9, 2210 (1974).
  * Kamimura  _et al._ (1986) M. Kamimura, M. Yahiro, Y. Iseri, S. Sakuragi, H. Kameyama, and M. Kawai, Prog. Theor. Phys. Suppl. 89, 1 (1986).
  * Austern  _et al._ (1987) N. Austern, Y. Iseri, M. Kamimura, M. Kawai, G. Rawitscher, and M. Yahiro, Phys. Rep. 154, 125 (1987).
  * Yahiro  _et al._ (2012) M. Yahiro, T. Matsumoto, K. Minomo, T. Sumi, and S. Watanabe, Prog. Theor. Phys. Supp. 196, 87 (2012).
  * Matsumoto  _et al._ (2004) T. Matsumoto, E. Hiyama, K. Ogata, Y. Iseri, M. Kamimura, S. Chiba, and M. Yahiro, Phys. Rev. C 70, 061601 (2004).
  * Descouvemont (2018) P. Descouvemont, Phys. Rev. C 97, 064607 (2018).
  * de Diego  _et al._ (2014) R. de Diego, J. M. Arias, J. A. Lay, and A. M. Moro, Phys. Rev. C 89, 064609 (2014).
  * Pesudo  _et al._ (2017) V. Pesudo, M. J. G. Borge, A. M. Moro, J. A. Lay, E. Nácher, J. Gómez-Camacho, O. Tengblad, L. Acosta, M. Alcorta, M. A. G. Alvarez, C. Andreoiu, P. C. Bender, R. Braid, M. Cubero, A. Di Pietro, J. P. Fernández-García, P. Figuera, M. Fisichella, B. R. Fulton, A. B. Garnsworthy, G. Hackman, U. Hager, O. S. Kirsebom, K. Kuhn, M. Lattuada, G. Marquínez-Durán, I. Martel, D. Miller, M. Moukaddam, P. D. O’Malley, A. Perea, M. M. Rajabali, A. M. Sánchez-Benítez, F. Sarazin, V. Scuderi, C. E. Svensson, C. Unsworth, and Z. M. Wang, Phys. Rev. Lett. 118, 152502 (2017).
  * Ogata  _et al._ (2009) K. Ogata, M. Kan, and M. Kamimura, Prog. Theor. Phys. 122, 1055 (2009).
  * Descouvemont (2017) P. Descouvemont, Phys. Lett. B 772, 1 (2017).
  * von Oertzen (1970) W. von Oertzen, Nucl. Phys. A 148, 529 (1970).
  * Imanishi and von Oertzen (1987) B. Imanishi and W. von Oertzen, Phys. Rep. 155, 29 (1987).
  * Buttle and Goldfarb (1966) P. Buttle and L. Goldfarb, Nucl. Phys. 78, 409 (1966).
  * Sparenberg  _et al._ (2000) J.-M. Sparenberg, D. Baye, and B. Imanishi, Phys. Rev. C 61, 054610 (2000).
  * Descouvemont and Baye (2010) P. Descouvemont and D. Baye, Rep. Prog. Phys. 73, 036301 (2010).
  * Baye (2015) D. Baye, Phys. Rep. 565, 1 (2015).
  * Dohet-Eraly (2017) J. Dohet-Eraly, Eur. Phys. J. Plus 132, 362 (2017).
  * Cotanch (1975) S. Cotanch, Phys. Lett. B 57, 123 (1975).
  * Satchler (1983) G. R. Satchler,  _Direct Nuclear Reactions_ (Oxford University Press, 1983).
  * Thompson (1988) I. J. Thompson, Comput. Phys. Rep. 7, 167 (1988).
  * Shubhchintak and Descouvemont (2019) Shubhchintak and P. Descouvemont, Phys. Rev. C 100, 034611 (2019).
  * Canto and Hussein (2013) L. F. Canto and M. S. Hussein,  _Scattering Theory of Molecules, Atoms and Nuclei_ (World Scientific Publishing, Singapore, 2013).
  * Phuc  _et al._ (2018) N. T. T. Phuc, N. H. Phuc, and D. T. Khoa, Phys. Rev. C 98, 024613 (2018).
  * Phuc  _et al._ (2021) N. H. Phuc, D. T. Khoa, and N. T. T. Phuc, Eur. Phys. J. A 57, 7 (2021).
  * Phuc  _et al._ (2019) N. T. T. Phuc, R. S. Mackintosh, N. H. Phuc, and D. T. Khoa, Phys. Rev. C 100, 054615 (2019).
  * Liénard  _et al._ (1995) E. Liénard, D. Baye, T. Delbar, P. Descouvemont, P. Duhamel, W. Galster, M. Kurokawa, P. Leleux, I. Licot, P. Lipnik, C. Michotte, T. Motobayashi, A. Ninane, J. Vanhorenbeeck, and J. Vervier, Phys. Rev. C 52, 775 (1995).
  * Voit  _et al._ (1988) H. Voit, N. Bischof, W. Tiereth, I. Weitzenfelder, W. von Oertzen, and B. Imanishi, Nucl. Phys. A 476, 491 (1988).
  * Imanishi  _et al._ (1997) B. Imanishi, V. Denisov, and T. Motobayashi, Phys. Rev. C 55, 1946 (1997).
  * Baye (1986) D. Baye, Nucl. Phys. A 460, 581 (1986).
  * Treu  _et al._ (1980) W. Treu, H. Fröhlich, W. Galster, P. Dück, and H. Voit, Phys. Rev. C 22, 2462 (1980).
  * von Oertzen and Imanishi (1984) W. von Oertzen and B. Imanishi, Nucl. Phys. A 424, 262 (1984).
  * Koning and Delaroche (2003) A. J. Koning and J. P. Delaroche, Nucl. Phys. A 713, 231 (2003).
  * Fröhlich  _et al._ (1984) H. Fröhlich, N. Bischof, W. Tiereth, H. Voit, W. von Oertzen, and B. Imanishi, Nucl. Phys. A 420, 124 (1984).
  * Villari  _et al._ (1989) A. C. C. Villari, A. Lépine-Szily, R. L. Filho, O. P. Filho, M. M. Obuti, J. M. O. Jr., and N. Added, Nucl. Phys. A 501, 605 (1989).
  * Oulebsir  _et al._ (2012) N. Oulebsir, F. Hammache, P. Roussel, M. G. Pellegriti, L. Audouin, D. Beaumel, A. Bouda, P. Descouvemont, S. Fortier, L. Gaudefroy, J. Kiener, A. Lefebvre-Schuhl, and V. Tatischeff, Phys. Rev. C 85, 035804 (2012).
  * Avrigeanu  _et al._ (2009) M. Avrigeanu, A. Obreja, F. Roman, V. Avrigeanu, and W. von Oertzen, At. Data Nucl. Data Tables 95, 501 (2009).
  * Thompson  _et al._ (1989) I. Thompson, M. Nagarajan, J. Lilley, and M. Smithson, Nucl. Phys. A 505, 84 (1989).
  * Mackintosh (2019) R. S. Mackintosh, Eur. Phys. J. A 55, 147 (2019).
  * Abramowitz and Stegun (1972) M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions_ (Dover, London, 1972).
  * Buttle and Goldfarb (1968) P. Buttle and L. Goldfarb, Nucl. Phys. A 115, 461 (1968).
  * Perey and Buck (1962) F. Perey and B. Buck, Nucl. Phys. 32, 353 (1962).


# Efficient Reduced Magnetic Vector Potential Formulation for the Magnetic
Field Simulation of Superconducting Magnets

Laura A. M. D’Angelo1,2, Dominik Moll1, and Herbert De Gersem1,2 Manuscript
received Month Day, Year; revised Month Day, Year. Corresponding author: Laura
A. M. D’Angelo (email: dangelo@temf.tu-darmstadt.de). 1Institute for
Accelerator Science and Electromagnetic Fields, Technische Universität
Darmstadt, Darmstadt, Germany  2Graduate School Computational Engineering,
Technische Universität Darmstadt, Darmstadt, Germany

###### Abstract

The major advantage of reduced magnetic vector potential formulations (RMVPs)
is that complicated coil structures do not need to be resolved by a
computational mesh. Instead, they are modeled by thin wires, whose source
field is included into the simulation model along Biot-Savart’s law. Such an
approach has already been successfully employed in ROXIE for the simulation of
superconducting Large Hadron Collider magnets at CERN. This work presents an
updated RMVP approach, which significantly outperforms the original method.
The updated formulation is postulated, implemented, verified, compared to the
original formulation, and applied for the simulation of a quadrupole magnet.
The promising results of this work encourage further investigation towards an
updated simulation framework for next-generation accelerator magnets.

###### Index Terms:

Accelerator magnets, Biot-Savart law, finite element analysis, superconducting
coils

##  I Introduction

Without doubt, high-temperature superconducting (HTS) technology will lift
next-generation synchrotrons beyond today’s technological frontier [1]. With
this step forward, however, magnet design is confronted with new challenges
regarding the design of large, high-field and high-quality magnet systems.
Computer-aided design and numerical field simulation generally play a crucial
role in designing and optimizing accelerator magnet systems, and will do even
more so regarding the next-generation HTS magnet systems. For the past
decades, the simulation software ROXIE [2] proved to be an indispensable
workhorse for designing the low-temperature superconducting (LTS) magnets of
the Large Hadron Collider (LHC). ROXIE combines a hybrid finite-element (FE)
boundary-element method with a reduced magnetic vector potential (RMVP)
formulation [3], leading to fast and accurate simulations [4]. Herein, the
coils are modeled as thin wires, and their excitation is included as a source
magnetic field, which is calculated by Biot-Savart’s law. The major advantage
of this formulation is that these wires do not need to be resolved by a
computational mesh. Especially superconducting accelerator magnets, which
typically contain hundreds or thousands of coil windings, greatly benefit from
this approach.

Nonetheless, ROXIE and commercial out-of-the-box simulation tools struggle
with the multi-scale nature that is imposed by HTS tapes, resulting into
excessive computation times [5]. The goal of this work is to improve the
computational efficiency of the RMVP approach and to contribute towards
suitable simulation tools for future HTS magnet design campaigns. To do so,
ROXIE’s original formulation is first adapted to a pure FE model in Sec. II,
and then reformulated in order to diminish the number of Biot-Savart integrals
to be calculated to quantify the source magnetic field. This leads to a multi-
step calculation procedure, which is presented in Sec. III and demonstrated
for a case study with an eccentric line current in an iron tube. In Sec. IV,
the updated RMVP formulation is analyzed regarding accuracy and performance.
Herein, the proposed procedure proves to be clearly superior to ROXIE’s
original formulation. Sec. V finally showcases the updated RMVP formulation by
applying it to the two-dimensional (2D) nonlinear magnetostatic simulation of
the LHC’s LTS MQXA quadrupole magnet [6]. All simulations have been carried
out using the freely available FE solver GetDP [7].

##  II Original RMVP Formulation

In this section, the original RMVP formulation from [3] is recapitulated. The
original physical problem that has to be solved is the magnetostatic problem
with a homogeneous Dirichlet boundary condition, reading

$\displaystyle\nabla\times\left(\nu\nabla\times\vec{A}_{\mathrm{}}\right)$
$\displaystyle=\vec{J}$ $\displaystyle\text{in $V$},$ (1a)
$\displaystyle\vec{n}\times\vec{A}_{\mathrm{}}$ $\displaystyle=0$
$\displaystyle\text{on $\partial V$}.$ (1b)

Herein, $\vec{A}_{\mathrm{}}$ is the sought-for magnetic vector potential
(MVP), $\nu$ is the reluctivity and $\vec{J}$ is the current density, which
represents the excitation in this problem. The computational domain
$V=V_{\mathrm{a}}\cup V_{\mathrm{i}}$ consists of the coil and air domain
$V_{\mathrm{a}}$ and the iron domain $V_{\mathrm{i}}$, and $\partial V$ is its
boundary. Classically, the current excitation is modeled in the right-hand
side $\vec{J}$, e.g. by using winding functions [8]. This procedure requires
the explicit discretization of the individual wires (or at least half-turns)
in the FE mesh.

In contrast, the RMVP method represents the wires by one-dimensional curves,
which are not necessarily taken into account in the FE mesh. This benefits the
meshing workload in the overall simulation process. The magnetic vector
potential (MVP) is decomposed into

$\vec{A}_{\mathrm{}}=\vec{A}_{\mathrm{s}}+\vec{A}_{\mathrm{r}},$ (2)

where $\vec{A}_{\mathrm{s}}$ is called the source MVP, and
$\vec{A}_{\mathrm{r}}$ the reduced MVP. $\vec{A}_{\mathrm{s}}$ is obtained by
evaluating Biot-Savart’s law [9]

$\vec{A}_{\mathrm{s}}=\frac{\mu_{0}}{4\pi}\int\limits_{\mathcal{L}^{\prime}}\frac{I\,\mathrm{d}\vec{s}^{\prime}}{|\vec{r}-\vec{r}^{\prime}|}$
(3)

for _all_ spatial coordinates $\vec{r}\in V$. The source domain
$\mathcal{L}^{\prime}$ represents the line, on which the line current $I$ is
located. In a three-dimensional (3D) model, this would be an arbitrarily
complicated curve loop in space; while in a two-dimensional (2D) setting,
$\mathcal{L}^{\prime}$ is reduced to a point, which represents a line current
going in or out of plane. Furthermore, in 2D, the MVP is assumed to have only
a $z$-component, $\vec{A}_{\mathrm{}}=A_{z}(x,y)\vec{e}_{z}$, and Biot-
Savart’s law (3) becomes

$A_{z}=\frac{\mu_{0}}{2\pi}\int\limits_{\mathcal{L}^{\prime}}I\ln(\left|\vec{r}-\vec{r}^{\prime}\right|^{-1})\,\mathrm{d}r^{\prime}.$
(4)

Multiple sources are taken into account by superposition.

Eventually, the discrete source MVP

$\vec{A}_{\mathrm{s}}(\vec{r})\approx\sum\limits_{j=1}^{N_{\text{edge}}}\overset{\mathchoice{}{}{\mbox{\raisebox{-2.15277pt}[0.0pt][0.0pt]{$\scriptstyle\smallfrown$}}}{\mbox{\raisebox{-1.50696pt}[0.0pt][0.0pt]{$\scriptscriptstyle\smallfrown$}}}}{a}_{\mathrm{s},j}\,\vec{w}_{j}(\vec{r})$
(5)

can be computed on the mesh edges $j=1,\dots,N_{\text{edge}}$ in two ways: One
can utilize the partition-of-unity property and calculate the discrete
coefficients
$\overset{\mathchoice{}{}{\mbox{\raisebox{-2.15277pt}[0.0pt][0.0pt]{$\scriptstyle\smallfrown$}}}{\mbox{\raisebox{-1.50696pt}[0.0pt][0.0pt]{$\scriptscriptstyle\smallfrown$}}}}{a}_{\mathrm{s},j}$
per edge $e_{j}$ directly by weighting (3) with the $j$-th edge function
$\vec{w}_{j}$ and integrating over that edge $e_{j}$,

$\overset{\mathchoice{}{}{\mbox{\raisebox{-2.15277pt}[0.0pt][0.0pt]{$\scriptstyle\smallfrown$}}}{\mbox{\raisebox{-1.50696pt}[0.0pt][0.0pt]{$\scriptscriptstyle\smallfrown$}}}}{a}_{\mathrm{s},j}=\frac{\mu_{0}}{4\pi}\int\limits_{e_{j}}\int\limits_{\mathcal{L}^{\prime}}\frac{I\,\mathrm{d}\vec{s}^{\prime}}{|\vec{r}-\vec{r}^{\prime}|}\cdot\vec{w}_{j}\,\mathrm{d}s.$
(6)

Alternatively, one performs a weak $L^{2}$-projection of the Biot-Savart
integral onto $\vec{A}_{\mathrm{s}}$,

$(\vec{A}_{\mathrm{s}},\vec{A}_{\mathrm{s}}^{\prime})_{V}=\left(\frac{\mu_{0}}{4\pi}\int\limits_{\mathcal{L}^{\prime}}\frac{I\,\mathrm{d}\vec{s}^{\prime}}{|\vec{r}-\vec{r}^{\prime}|},\vec{A}_{\mathrm{s}}^{\prime}\right)_{V}$
(7)

with test functions $\vec{A}_{\mathrm{s}}^{\prime}\in H(\text{curl};V)$ in the
Hilbert space [10]

$H(\text{curl};V):=\\{\vec{A}_{\mathrm{}}\in L^{2}(V)\leavevmode\nobreak\
:\leavevmode\nobreak\ \nabla\times\vec{A}_{\mathrm{}}\in L^{2}(V)\\}.$ (8)

This work uses the built-in $L^{2}$-projection of GetDP.

The reduced MVP $\vec{A}_{\mathrm{r}}$ is computed by solving the boundary
value problem (BVP)

$\displaystyle\nabla\times(\nu\nabla\times\vec{A}_{\mathrm{r}})$
$\displaystyle=-\nabla\times(\nu\nabla\times\vec{A}_{\mathrm{s}})$
$\displaystyle\text{in $V_{\mathrm{i}}$},$ (9a)
$\displaystyle\nabla\times(\nu_{0}\nabla\times\vec{A}_{\mathrm{r}})$
$\displaystyle=0$ $\displaystyle\text{in $V_{\mathrm{a}}$},$ (9b)
$\displaystyle\vec{n}\times\vec{A}_{\mathrm{r}}$
$\displaystyle=-\vec{n}\times\vec{A}_{\mathrm{s}}$ $\displaystyle\text{on
$\partial V$}.$ (9c)

Using a Ritz-Galerkin approach, the weak formulation is obtained as: Find
$\vec{A}_{\mathrm{r}}\in H_{\mathrm{r}}(\text{curl};V)$, s.t.

$(\nu\nabla\times\vec{A}_{\mathrm{r}},\nabla\times\vec{A}_{\mathrm{r}}^{\prime})_{V}=-(\nu\nabla\times\vec{A}_{\mathrm{s}},\nabla\times\vec{A}_{\mathrm{r}}^{\prime})_{V_{\mathrm{i}}}$
(10)

$\forall\vec{A}_{\mathrm{r}}^{\prime}\in H_{\mathrm{r}}(\text{curl};V)$, where
$\vec{A}_{\mathrm{r}}^{\prime}$ is a test function, and

$H_{\mathrm{r}}(\text{curl};V)=\\{\vec{A}_{\mathrm{}}\in
H(\text{curl};V):\gamma_{\partial
V}(\vec{A}_{\mathrm{}})=-\vec{n}\times\vec{A}_{\mathrm{s}}\\}$ (11)

is chosen in order to fulfill (9c). Herein,

$\gamma_{\mathcal{B}}(\vec{A}_{\mathrm{}})=\vec{n}\times\vec{A}_{\mathrm{}}|_{\mathcal{B}}$
(12)

is the tangential trace operator w.r.t. a boundary $\mathcal{B}$ [10].

The weak formulation (10) is solved by a FE method employing standard edge
shape functions. Lastly, the total MVP $\vec{A}_{\mathrm{}}$ is composed of
$\vec{A}_{\mathrm{s}}$ and $\vec{A}_{\mathrm{r}}$ following (2). Note that the
Biot-Savart integral (3) must be evaluated in the whole domain $V$, whether
one is actually interested in the solution of the whole domain or of only a
small sub-domain.

##  III Updated RMVP Formulation

###  III-A Ansatz and Biot-Savart sub-formulation

The domain $V$ is decomposed into a non-permeable sub-domain $V_{\mathrm{a}}$
(consisting of e.g. air) and a source-free sub-domain $V_{\mathrm{i}}$
(containing e.g. iron). Herein, $V_{\mathrm{a}}$ is fully embedded into
$V_{\mathrm{i}}$, and $\Gamma$ denotes the interface between those domains.
The MVP $\vec{A}_{\mathrm{}}$ is then decomposed into

$\vec{A}_{\mathrm{}}=\left\\{\begin{array}[]{ll}\vec{A}_{\mathrm{g}}+\vec{A}_{\mathrm{s}}+\vec{A}_{\mathrm{m}}&\text{in
$V_{\mathrm{a}}$},\\\ \vec{A}_{\mathrm{g}}&\text{in
$V_{\mathrm{i}}$},\end{array}\right.$ (13)

where $\vec{A}_{\mathrm{s}}$ is called the source MVP, $\vec{A}_{\mathrm{m}}$
the image MVP, and $\vec{A}_{\mathrm{g}}$ the reaction MVP.
$\vec{A}_{\mathrm{s}}$ is obtained by evaluating Biot-Savart’s law (3) for all
spatial coordinates $\vec{r}\in V_{\mathrm{a}}$ that are of interest, but at
least for $\vec{r}\in\Gamma$. This is the major advantage to the original RMVP
approach [3], where $\vec{A}_{\mathrm{s}}$ needed to be calculated in the
whole domain $V$.

###  III-B Image sub-formulation

The image MVP $\vec{A}_{\mathrm{m}}$ is the solution of the BVP

$\displaystyle\nabla\times(\nu_{0}\nabla\times\vec{A}_{\mathrm{m}})$
$\displaystyle=0$ $\displaystyle\text{in $V_{\mathrm{a}}$},$ (14a)
$\displaystyle\vec{n}\times\vec{A}_{\mathrm{m}}+\vec{n}\times\vec{A}_{\mathrm{s}}$
$\displaystyle=0$ $\displaystyle\text{on $\Gamma$}.$ (14b)

The MVP $\vec{A}_{\mathrm{s}}+\vec{A}_{\mathrm{m}}$ is the equivalent of a
Green’s function obeying a homogeneous Dirichlet boundary condition at
$\Gamma$ [9]. Using a Ritz-Galerkin approach and the abbreviation
$\vec{n}\times\vec{H}_{\mathrm{m}}=\vec{n}\times\vec{H}_{\mathrm{m}}$ for the
tangential component of the image magnetic field strength, the weak
formulation is obtained as:
Find $\vec{A}_{\mathrm{m}}\in
H(\mathrm{curl};V_{\mathrm{a}}),\vec{n}\times\vec{H}_{\mathrm{m}}\in
H^{-1/2}(\mathrm{curl};\Gamma)$, s.t.

$\displaystyle(\nu_{0}\nabla\times\vec{A}_{\mathrm{m}},\nabla\times\vec{A}_{\mathrm{m}}^{\prime})_{V_{\mathrm{a}}}+(\vec{n}\times\vec{H}_{\mathrm{m}},\vec{A}_{\mathrm{m}}^{\prime})_{\Gamma}$
$\displaystyle=0,$ (15a)
$\displaystyle(\vec{A}_{\mathrm{m}},\vec{n}\times\vec{H}_{\mathrm{m}}^{\prime})_{\Gamma}+(\vec{A}_{\mathrm{s}},\vec{n}\times\vec{H}_{\mathrm{m}}^{\prime})_{\Gamma}$
$\displaystyle=0$ (15b)

$\forall\vec{A}_{\mathrm{m}}^{\prime}\in
H(\mathrm{curl};V_{\mathrm{a}}),\forall\vec{n}\times\vec{H}_{\mathrm{m}}^{\prime}\in
H^{-1/2}(\mathrm{curl};\Gamma)$, where $\vec{A}_{\mathrm{m}}$ is the discrete
image MVP, $\vec{A}_{\mathrm{m}}^{\prime}$ and
$\vec{n}\times\vec{H}_{\mathrm{m}}^{\prime}$ are corresponding test functions,
and $H^{-1/2}(\mathrm{curl};\Gamma)$ is a trace space [10]. Here, the boundary
condition (14b) is weakly imposed by (15b), yielding a saddle-point problem
[11]. In this way, the quantity $\vec{n}\times\vec{H}_{\mathrm{m}}$, which
will be needed for the calculation of $\vec{A}_{\mathrm{g}}$, is already at
hand. The weak sub-formulation (15) is eventually solved by a FE method
employing standard edge shape functions.

###  III-C Reaction sub-formulation

The reaction MVP $\vec{A}_{\mathrm{g}}$ is obtained by solving the BVP

$\displaystyle\nabla\times(\nu\nabla\times\vec{A}_{\mathrm{g}})$
$\displaystyle=\vec{J}_{\mathrm{g}}$ $\displaystyle\text{in $V$},$ (16a)
$\displaystyle\vec{n}\times\vec{A}_{\mathrm{g}}$ $\displaystyle=0$ on
$\partial V$ (16b)

with $\vec{J}_{\mathrm{g}}=\vec{K}_{\mathrm{g}}\,\delta_{\Gamma}$, where

$\vec{K}_{\mathrm{g}}=\vec{n}\times(\vec{H}_{\mathrm{s}}+\vec{H}_{\mathrm{m}})$
(17)

denotes the surface current density in A/m, and $\delta_{\Gamma}$ the delta
distribution function defined by

$\int\limits_{V}f\,\delta_{\Gamma}\,\mathrm{d}V=\int\limits_{\Gamma}f\,\mathrm{d}S\quad\forall
f.$ (18)

Here, $\vec{H}_{\mathrm{s}}=\nu_{0}\nabla\times\vec{A}_{\mathrm{s}}$ and
$\vec{H}_{\mathrm{m}}=\nu_{0}\nabla\times\vec{A}_{\mathrm{m}}$ are the source
and image magnetic field strengths, respectively. Visually, the current
excitation in $V_{\mathrm{a}}$ has been shifted onto the interface surface
$\Gamma$. The weak formulation of (16) reads:
Find $\vec{A}_{\mathrm{g}}\in H_{0}(\mathrm{curl};V)$ s.t.

$\displaystyle(\nu\nabla\times\vec{A}_{\mathrm{g}},\nabla\times\vec{A}_{\mathrm{g}}^{\prime})_{V}=$
(19)
$\displaystyle(\vec{n}\times\vec{H}_{\mathrm{s}},\vec{A}_{\mathrm{g}}^{\prime})_{\Gamma}$
$\displaystyle+(\vec{n}\times\vec{H}_{\mathrm{m}},\vec{A}_{\mathrm{g}}^{\prime})_{\Gamma}$

$\forall\vec{A}_{\mathrm{g}}^{\prime}\in H_{0}(\mathrm{curl};V)$, where
$\vec{A}_{\mathrm{g}}^{\prime}$ is a test function, and

$H_{0}(\mathrm{curl};V)=\\{\vec{A}_{\mathrm{}}\in
H(\mathrm{curl};V)\leavevmode\nobreak\ :\leavevmode\nobreak\ \gamma_{\partial
V}(\vec{A}_{\mathrm{}})=0\\}$ (20)

is chosen in order to fulfill (16b). The weak sub-formulation (19) is solved
by a FE method employing standard edge shape functions. Finally, the total MVP
$\vec{A}_{\mathrm{}}$ is composed of $\vec{A}_{\mathrm{g}}$,
$\vec{A}_{\mathrm{m}}$ and $\vec{A}_{\mathrm{s}}$ following (13). Note that
$\vec{A}_{\mathrm{g}}$ corresponds to the full MVP in the domain
$V_{\mathrm{i}}$, which facilitates the treatment of nonlinearities and eddy
currents.

$I$$\Gamma$$V_{\mathrm{a}}$$V_{\mathrm{i}}$ Figure 1: 2D case study model: An
infinitely long eccentric line current $I$ in an air domain $V_{\mathrm{a}}$
surrounded by an infinitely long iron tube $V_{\mathrm{i}}$.

(a) Equipotential lines of the source MVP $\vec{A}_{\mathrm{s}}$ in the
eccentric wire case study calculated by (3). The source current is depicted as
a green line. For the sake of visualization, $\vec{A}_{\mathrm{s}}$ is
calculated in the whole domain instead of only at the interface boundary.

(b) Equipotential lines of the image MVP $\vec{A}_{\mathrm{m}}$ in the
eccentric wire case study calculated by (15). The computation is only done in
the sub-domain $V_{\mathrm{a}}$.

(c) The green lines visualize the surface current density
$\vec{J}_{\mathrm{g}}$ at $\Gamma$ determined by (17). The equipotential lines
of the reaction MVP $\vec{A}_{\mathrm{g}}$ in the eccentric wire case study
calculated by (19).

(d) Equipotential lines of the total MVP
$\vec{A}=\vec{A}_{\mathrm{s}}+\vec{A}_{\mathrm{m}}+\vec{A}_{\mathrm{g}}$ in
the eccentric wire case study. The source current is depicted as a green line.

Figure 2: The MVP components obtained in the sub-formulations of the RMVP
formulation and the resulting total MVP for the case study of the eccentric
wire in an iron ring.

###  III-D Case study: Eccentric line current in an iron tube

Consider as a case study an infinitely long eccentric line current $I$ in air
surrounded by an infinitely long iron tube. This model’s 2D cross-section is
shown in Fig. 1. $V_{\mathrm{a}}$ is defined as the air and coil region within
the iron tube, while the iron tube itself is chosen as $V_{\mathrm{i}}$. Thus,
the interface $\Gamma$ is the circular boundary between those two domains.

The RMVP formulation is applied to compute the MVP and magnetic flux density
caused by the direct current $I$. Fig. 2(a) shows the source MVP
$\vec{A}_{\mathrm{s}}$, which is obtained by evaluating the Biot-Savart
integral (3). It represents the MVP as if the eccentric line current was
located in free space. The corresponding image MVP $\vec{A}_{\mathrm{m}}$
calculated by (15) is seen in Fig. 2(b). Fig. 2(c) visualizes the surface
current density $\vec{K}_{\mathrm{g}}$ determined by (17), from which the
reaction MVP $\vec{A}_{\mathrm{g}}$ calculated using (19) is shown as flux
lines in the same figure. The superposition of these three sub-solutions leads
to the total MVP $\vec{A}$ as shown in Fig. 2(d).

##  IV Numerical Studies

###  IV-A Benchmark model: Racetrack coil

The RMVP formulation is implemented in the freely available open-source FE
solver GetDP [7] and employed to carry out a 2D linear magnetostatic
simulation of a racetrack coil in an iron yoke. Fig. 3(a) shows the geometry,
which consists of two winding groups (hatched rectangles) containing wires
with rectangular cross-section and embedded in an air domain $V_{\mathrm{a}}$
(white), which is surrounded by an iron yoke $V_{\mathrm{i}}$ (gray). $\Gamma$
is the interface between $V_{\mathrm{a}}$ and $V_{\mathrm{i}}$. For the iron
yoke, a linear permeability of $\mu_{\mathrm{i}}=4000\mu_{0}$ has been chosen.
Since the line currents come together with singularities, the magnetic energy
strives for infinity. Therefore, the magnetic energy considered below is
evaluated in a sub-domain, which is the dashed circular domain
$V_{\text{eval}}$ within $V_{\mathrm{a}}$.

The configuration of each winding group is shown in Fig. 3(b). In the 2D
reference model, each half-turn of the coil is modeled as a surface (black
rectangles in Fig. 3(b)) with a given surface current density. In the 2D RMVP
setting, each half-turn is discretized by a set of points (see Fig. 3(b), gray
and purple dots), which represent line currents in or out of plane, depending
on the current orientation. In the simplest case, a half-turn is modeled by a
single line current located in the cross-section center as seen in Fig. 3(b)
(purple dots).

$\Gamma$$V_{\mathrm{i}}$$V_{\mathrm{a}}$$V_{\text{eval}}$ (a) Two coil
windings (hatched rectangles, see Fig. 3(b) for a detailed sketch) in air
(white, denoted as $V_{\mathrm{a}}$) surrounded by an iron yoke (gray, denoted
as $V_{\mathrm{i}}$). The interface boundary between iron and air is denoted
as $\Gamma$. The dashed circular domain represents the evaluation domain for
the magnetic energy.

(b) Detailed view of the coil winding configuration for $N_{x}=N_{y}=3$. In
the volumetric case (1), each winding is modeled as a rectangle with a given
volumetric current density. In the RMVP case (19), each half-turn is
represented by a set of points (purple and gray dots) representing line
currents in or out of plane, depending on the current orientation. In the
simplest case, each half-turn is modeled by a single line current located in
the cross-section center (purple dots).

Figure 3: Racetrack coil model used for the purpose of numerical studies.
Figure 4: Magnetic flux density magnitude in the racetrack coil obtained by
the updated RMVP formulation.

The results of the RMVP approach are verified against reference values
obtained by a conventional 2D FE simulation modeling the half-turns as
surfaces and solving (1) while using a very fine mesh. The resulting magnetic
flux density magnitude of the proposed method is shown in Fig. 4, where the
maximal value has been capped to that of the reference solution. Particularly
high fields are obtained at the four inner corners of the iron yoke due to the
sharp geometry as well as the linear material properties and, of course, at
the line currents due to their singular nature.

###  IV-B Convergence Analysis

$10^{-2}$$10^{-1}$$10^{-2}$$10^{-1}$$10^{0}$$\mathcal{O}(h)$Mesh length
$h$Relative error$L^{2}(V_{\text{eval}})$$L^{2}(V)$ Figure 5: Linear
convergence of the $L^{2}$-error in the evaluation domain (solid, purple)
w.r.t. the mesh length for the RMVP formulation. The $L^{2}$-error in the
whole domain (dashed, orange) does not converge due to the singularities.

For the convergence analysis, the $L^{2}$-error w.r.t. the volumetric
reference solution is considered. The line currents have to be excluded from
the considered domain, as otherwise the $L^{2}$-error would not converge due
to the singularities. Therefore, a sub-domain $V_{\text{eval}}\subset V$ is
chosen, such that $\mathcal{L}^{\prime}\nsubseteq V_{\text{eval}}$ as seen in
Fig. 3(a) (circular region).

For lowest order FEs, one would normally expect a quadratic convergence of the
$L^{2}$-error. However, this only holds if, among other criteria, the right-
hand side of the BVP is in $L^{2}(V)$ [12]. In the reaction sub-problem (16),
the right-hand side is a Dirac source term, making the problem not regular.
Although classical convergence results are therefore invalid here, it has been
shown that the $L_{2}$-error of 2D elliptic problems with Dirac source terms
converges linearly [13]. Indeed, this is observed for the
$L_{2}(V_{\text{eval}})$-error in Fig. 5 (solid purple line).

Varying the size of $V_{\text{eval}}$, and thereby the remoteness of $\partial
V_{\text{eval}}$ to the line currents, neither improves nor worsens the linear
convergence behavior, as long as the $V_{\text{eval}}$ excludes the line
currents or their immediate neighborhood. Then, the magnetic energy gets
highly overestimated and does not converge at all (see Fig. 5, dashed orange
line).

###  IV-C Performance comparison to the original formulation

The runtimes of the original and updated RMVP formulations are compared for
the 2D magnetostatic linear simulation of the benchmark model with more than
$100,000$ degrees of freedom on a standard workstation. Both the total
runtimes and the runtimes of only the Biot-Savart integral evaluation are
measured. Additionally, the updated RMVP simulation is carried out for the
optimal case that only the source MVP on $\Gamma$ has to be computed as well
as the worst case, in which the source MVP in $V_{\mathrm{a}}$ is wanted. The
latter would be the case if one is interested in the magnetic field in the
whole aperture of a magnet. In practice, magnet designers are often interested
in field values at particular points, e.g. in a circular curve for a multipole
coefficient analysis to investigate the field quality of the magnet [14].
Therefore, the average runtime is expected to range between those two extreme
cases.

Figure 6 shows the measured total runtimes (purple) and the portions of the
Biot-Savart law evaluation therein (orange). The following two observations
are made, from which two important conclusions are drawn:

  1. 1.

The Biot-Savart integral computation heavily dominates the total runtime of
both RMVP formulations. Therefore, the performance of the RMVP procedures can
be easily and significantly improved by parallelizing the integral evaluations
and employing fast-multipole methods [15].

  2. 2.

The updated RMVP formulation is by far computationally superior to the
original formulation. This holds true even for the worst case. It is even more
obvious when used for the calculation of the magnetic field at particular
points, along certain curves or in small-sized regions in the magnet’s
aperture. Then, the updated method impressively outshines the original
procedure.

Hence, Fig. 6 illustrates the performance gain of the updated RMVP formulation
with respect to the standard formulation and at the same time, parallelization
of the Biot-Savart solver as a straightforward and promising measure for
further improvement.

$0$$500$$1{,}000$$1{,}500$$2{,}000$$2{,}500$$3{,}000$$3{,}500$$4{,}000$$4{,}500$Updated
(on $\Gamma$)Updated (on $V_{\mathrm{a}}$)Original (on $V_{\mathrm{a}}\cup
V_{\mathrm{i}}$)$3581\,$s$3563\,$s$467\,$s$461\,$s$12\,$s$7\,$sTime (s) Figure
6: Runtime comparison between the original and updated RMVP. Both the total
times (purple) and the proportions of the Biot-Savart integral evaluations
alone (orange) are depicted. For the updated RMVP case, the worst and best
case with computation of $\vec{A}_{\text{s}}$ in $V_{\mathrm{a}}$ and $\Gamma$
has been considered, respectively.

##  V Simulation of the MQXA Quadrupole

Figure 7: Magnetic flux density magnitude in the MQXA quadrupole computed with
the RMVP formulation. Figure 8: Relative difference between the magnetic flux
density magnitude computed by the RMVP formulation and of the reference
simulation performed by a conventional 2D FE solver.

The proposed RMVP is employed to carry out a 2D magnetostatic nonlinear
simulation of the MQXA low-beta quadrupole of the LHC [6] utilizing GetDP and
FiQuS [16]. Each winding is approximated by one line current in its center,
resulting into 488 line currents in total. Fig. 7 shows the computed magnetic
flux density in the magnet. The results are compared to a reference simulation
performed by a conventional 2D FE method taking the windings into account by
surface current densities. Fig. 8 shows the relative difference
$\epsilon_{\mathrm{r}}$ of the magnetic flux density magnitude obtained by the
RMVP to that of the reference simulation. As expected, high discrepancies
($\epsilon_{\mathrm{r}}\geq 1$) occur in the neighborhood of the line
currents, which represent singularities (clearly visible in Fig. 8 as white
dots). Outside the immediate neighborhood of the singularities, a good
approximation is achieved with relative differences in the order of
$\epsilon_{\mathrm{r}}=10^{-3},\dots,10^{-1}$. An even better alignment
between the RMVP and reference simulation is expected if the longish winding
shape was resolved by multiple line currents instead of just one at its
center.

##  VI Conclusion

This work proposed an updated RMVP formulation for accurate and fast magnetic
field simulations of superconducting accelerator magnets. The formulation was
postulated and verified against a benchmark model. A runtime comparison showed
that the proposed method clearly outperforms the original formulation.
However, because of the Dirac source term occurring in one of the sub-
problems, the $L^{2}$-error only shows a linear convergence for lowest order
FEs instead of a quadratic one. This issue is well understood in the
scientific computing community, and different techniques exist to improve the
convergence order [17]. Finally, the updated RMVP procedure was embedded in
FiQuS, guaranteeing the greatest possible applicability in the accelerator
magnet engineer community, and successfully employed to carry out a 2D
nonlinear magnetostatic simulation of a quadrupole magnet.

The promising results of this work encourage an expansion of that method
towards an updated simulation framework for HTS accelerator magnets, including
physical formulations suitable for HTS magnets [18] and magnetization models
for HTS tapes [19].

## Acknowledgment

This work has been supported by CHART (http://chart.ch) in the context of the
MagNum project, by the German BMBF project BMBF-05P18RDRB1, by the Graduate
School Computational Engineering at TU Darmstadt, and by the DFG Research
Training Group 2128 "Accelerator Science and Technology for Energy Recovering
Linacs".

The authors thank Nicolas Marsic, Erik Schnaubelt, Andrea Vitrano and Bernhard
Auchmann for the fruitful discussions and their helpful advice.

## References

  * [1] Maria Durante et al. “Realization and First Test Results of the EuCARD 5.4-T REBCO Dipole Magnet” In _IEEE Trans. Appl. Supercond._ 28.3, 2018, pp. 1–5
  * [2] “ROXIE: Routine for the Optimization of Magnet X-Sections, Inverse Field Calculation and Coil End Design – First International ROXIE Users Meeting and Workshop”, 1998
  * [3] Christian Paul “Numerical calculation of the saturation induced field errors in the main dipole of CERN’s future Large Hadron Collider”, 1997
  * [4] S. Kurz, S. Russenschuck and N. Siegel “Accurate calculation of fringe fields in the LHC main dipoles” In _IEEE Trans. Appl. Supercond._ 10.1, 2000, pp. 85–88
  * [5] Francesco Grilli “Numerical Modeling of HTS Applications” In _IEEE Trans. Appl. Supercond._ 26.3, 2016
  * [6] R. Ostojic “The LHC insertion magnets” In _IEEE Trans. Appl. Supercond._ 12.1, 2002, pp. 196–201
  * [7] P. Dular, C. Geuzaine, F. Henrotte and W. Legros “A general environment for the treatment of discrete problems and its application to the finite element method” In _IEEE Trans. Magn._ 34.5, 1998, pp. 3395–3398
  * [8] Sebastian Schöps, Herbert De Gersem and Thomas Weiland “Winding functions in transient magnetoquasistatic field-circuit coupled simulations” In _COMPEL_ 32.6, 2013, pp. 2062–2082
  * [9] John David Jackson “Classical Electrodynamics” Hoboken, NJ, USA: John Wiley & Sons, 1998
  * [10] Peter Monk “Finite Element Methods for Maxwell’s Equations”, Numerical Mathematics and Scientific Computation New York, NY, USA: Oxford University Press, 2003
  * [11] Franco Brezzi and Michel Fortin “Mixed and Hybrid Finite Element Methods” New York, NY: Springer, 1991
  * [12] Alexandre Ern and Jean-Luc Guermond “Analysis of the edge finite element approximation of the Maxwell equations with low regularity solutions” In _Comput. Math. Appl._ 75.3, 2018, pp. 918–932
  * [13] Ridgway Scott “Finite Element Convergence for Singular Data” In _Numer. Math._ 21, 1973, pp. 317–327
  * [14] Stephan Russenschuck “Field Computation for Accelerator Magnets” Weinheim: Wiley-VCH, 2010
  * [15] Friedemann Groh et al. “Calculating Exciting Fields Using the Fast Multipole Method and an Integral Transformation to the Coil Surface” In _IEEE Trans. Magn._ 41.5, 2005, pp. 1384–1387
  * [16] Andrea Vitrano et al. “An Open-Source Finite Element Quench Simulation Tool for Superconducting Magnets” In _IEEE Trans. Appl. Supercond._ 33.5, 2023, pp. 1–6
  * [17] Silvia Bertoluzza, Astrid Decoene, Loïc Lacouture and Sébastien Martin “Local Error Estimates of the Finite Element Method for an Elliptic Problem with a Dirac Source Term” In _Numer. Methods Partial Differential Eq._ 34, 2018, pp. 97–120
  * [18] L. Bortot et al. “A Coupled A-H Formulation for Magneto-Thermal Transients in High-Temperature Superconducting Magnets” In _IEEE Trans. Appl. Supercond._ 30.5, 2020
  * [19] J. van Nugteren et al. “Measurement and Numerical Evaluation of AC Losses in a ReBCO Roebel Cable at 4.5 K” In _IEEE Trans. Appl. Supercond._ 26.3, 2016, pp. 1–7