# The 6Li$(p,\gamma)^{7}$Be reaction rate in the light of the new LUNA data

S. B. Dubovichenko  Fesenkov Astrophysical Institute ”NCSRT” ASA MDASI RK,
050020 Almaty, Kazakhstan  al-Farabi Kazakh National University, 050040
Almaty, Kazakhstan  A. S. Tkachenko  Fesenkov Astrophysical Institute ”NCSRT”
ASA MDASI RK, 050020 Almaty, Kazakhstan  R. Ya. Kezerashvili  New York City
College of Technology, City University of New York, Brooklyn, 11201 New York,
USA  Graduate School and University Center, City University of New York, 10016
New York, USA  N. A. Burkova  al-Farabi Kazakh National University, 050040
Almaty, Kazakhstan  A. V. Dzhazairov-Kakhramanov  Fesenkov Astrophysical
Institute ”NCSRT” ASA MDASI RK, 050020 Almaty, Kazakhstan

###### Abstract

We present new calculations of the astrophysical $S-$factor and reaction rate
for the 6Li$(p,\gamma)^{7}$Be reaction at energies of 10 keV to 5 MeV in the
framework of a modified potential cluster model with forbidden states,
including low lying resonances. The astrophysical $S(E)-$factor is compared
with the available experimental data and calculations done within different
models. The results for the $S-$factor are in good agreement with the data set
(for $E<0.3$ MeV) and calculations (for $E<0.6$ MeV) of LUNA collaboration
(Phys. Rev. C 102, 052802, 2020). The recommended extrapolated zero value
$S(0)$ turned out to be 101 eV $\cdot$ b. Using the theoretical total cross-
sections, the 6Li$(p,\gamma)^{7}$Be capture reaction rate is calculated at
temperatures ranging from 0.01 to 10 $T_{9}$ and compared with NACRE and NACRE
II. Analytical expressions for the $S-$factor and reaction rate are given, and
the effect of low-lying resonances on the reaction rate is estimated. We
suggest to update the NACRE and NACRE II databases in light of the new LUNA
data and present calculations.

low and astrophysical energies, $p^{6}$Li system, thermonuclear reaction rate,
potential cluster model

††preprint: APS/123-QED

##  I Introduction

The radiative 6Li$(p,\gamma)^{7}$Be capture reaction is of great interest in
nuclear astrophysics [1, 2]. Since 1955 the 6Li$(p,\gamma)^{7}$Be reaction at
low energies has been studied by several experimental groups [3, 4, 5, 6, 7,
8, 9, 10, 11]. Measurements of the astrophysical $S-$factor of this reaction
were limited to the energy range of 35 keV to 1.2 MeV. The astrophysical
$S-$factor and reaction rate were studied in the framework of different
theoretical approaches and methods [12, 13, 14, 15, 16, 17, 18, 11]. A
detailed review of the theoretical and experimental current status is given in
Ref. [9].

In 2020, new experimental data were obtained in the Laboratory for Underground
Nuclear Astrophysics (LUNA) [9], and it excluded the possibility of resonance
mentioned in [8]. It seems challenging to consider this reaction in the
astrophysical energy range, for which experimental data are available. In
particular, our interest is the re-examination of $S-$factor and has two foci:
i. to consider 6Li$(p,\gamma)^{7}$Be within the framework of the modified
potential cluster model with the classification of bound and scattering states
according to Young’s orbital diagrams [28]; ii. to describe the $S-$factor
using all available experimental data and to obtain the reaction rate. While
assessing the reliability of MPCM, it is reasonable to extend the energy
interval up to 5 MeV to estimate the role of resonances in this energy range.

Ten years ago, we studied this reaction, but limited ourselves to an energy
range up to 1 MeV. Moreover, we did not take into account resonances, nor did
we consider the reaction rate [14, 15].

In this paper, we investigate the energy dependence of the astrophysical
$S-$factor of $\ {}^{6}$Li$(p,\gamma)^{7}$Be reaction at energies of 10 keV to
5 MeV. By considering several resonances, including a wide resonance at
$E_{x}=9.9$ MeV, the reaction rate in the temperature range of 0.01 to 10
$T_{9}$ is calculated. We demonstrate that it is possible to correctly convey
the available experimental data based on potentials that are consistent with
the energies of bound states and their asymptotic constants. For the
scattering potentials, we use the parameters consistent with the resonance
spectrum of the final nucleus.

The results obtained for the reaction rate are approximated by curves of a
particular type to simplify their use in applied research. These results apply
to problems in nuclear astrophysics related to light atomic nuclei and ultra-
low energies.

This work is organized as follows. Secs. II and III present the theoretical
framework and fundamentals of the modified potential cluster model (MPCM),
constructing principles of discrete states potentials and the description of
the $p^{6}$Li channel in the continuous spectrum in the MPCM. Cross-section
for the radiative capture processes and discussion of wave functions
asymptotics are given in Sec. IV. Classification of cluster states and bound
states potentials are given in Sec. V, and Sec. VI contains the description of
the $p^{6}$Li channel in the continuous spectrum. Sec. VII is devoted to the
astrophysical $S-$factor and the 6Li$(p,\gamma)^{7}$Be reaction rate.
Appendices A and B include the description of the finite-difference method
that we are using in present calculations and a table of numerical values of
the $p^{6}$Li reaction rate in the temperature range of 0.001 $T_{9}$ to 10
$T_{9}$, respectively. We outline conclusions in Sec. VIII.

##  II Theoretical framework

Charged-particle induced reactions represent one of the main inputs in stellar
evolution. There are a number of theoretical methods used for description of
nuclear reactions at stellar energies that are based on fundamental principles
of quantum mechanics [20].

Since its first application in 1963 [21], the potential model approach has a
special place among models for description of low energy reactions. However,
over the course of over a half century, this model has been significantly
modified and improved. Below we present the fundamentals of the modified
potential cluster model (MPCM), where the Young diagrams are used for the
classification of orbital states and construction of potentials ([28] and
Refs. herein).

The basic features of the MPCM approach are as follows:

  1. 1.

The MPCM is a two-particle model that accounts for the internal
characteristics of clusters: their sizes, charges, masses, quadrupole and
magnetic momenta, which are used to calculate the reaction total cross-
sections or other characteristics of the final nucleus.

  2. 2.

The classification of cluster states is performed according to Young’s orbital
diagrams, leading to the concept of forbidden states in some partial waves
[28]. The Pauli principle is implemented via exclusion of the forbidden states
(FSs), manifesting in proper node behavior of the radial wave function (WF).
Forbidden states that lead to low-lying bound states are not physically
realized due to the orthogonality of corresponding functions and allowed state
functions.

  3. 3.

The Gaussian type inter-cluster interaction potentials are constructed, taking
into account these forbidden states in certain partial waves. For each partial
wave with specified quantum numbers, the potential is constructed with two
parameters, assuming it depends explicitly on Young’s orbital diagrams.

  4. 4.

Potentials of the bound states (BSs) are constructed based on asymptotic
constants (AC) and binding energies. Potentials of the scattering processes
are constructed based on the spectra of the final nucleus or the scattering
phase shifts of the particles of the input channel. Parameters of the
potentials are fixed or variable within the AC error intervals and vary within
the energy or width errors of resonant or excited states.

  5. 5.

The radial WFs of the allowed states of the continuous and discrete spectra
are tailored appropriately using correct asymptotics.

In light of the new experimental LUNA data [9], we reexamine the
6Li$(p,\gamma)^{7}$Be reaction $S-$factor and reaction rate within the
framework of the MPCM.

### Classification of cluster states

The total wave functions (WF) have the form of an antisymmetrized product of
completely antisymmetric internal wave functions of clusters
$\Psi(1,\dots,A_{1})=\Psi(\mathbf{R_{1}})$ and
$\Psi(A_{1}+1,\dots,A)=\Psi(\mathbf{\ R_{2}})$, multiplied by the
corresponding wave function $\Phi(\mathbf{R})$ of relative motion [22, 23, 24]

$\Psi=\hat{A}\\{\Psi(\mathbf{R_{1}})\Psi(\mathbf{R_{2}})\Phi(\mathbf{R})\\}.$
(1)

In Eq. (1) $\hat{A}$ is the antisymmetrization operator that permutes nucleons
from the clusters $A_{1}$ and $A_{2}$, $\mathbf{\ R_{1}}$ and $\mathbf{R_{2}}$
are the center-of-mass radius-vectors of the clusters, and
$\mathbf{R=R}_{1}\mathbf{-R}_{2}$ is the relative motion coordinate.

The wave functions (1) are characterized by specific quantum numbers,
including $JLS$ — total momentum, orbital quantum momentum, and spin,
respectively — and Young’s diagrams $\\{f\\}$, which determine the orbital
part of WF permutation symmetry of the relative motion of the clusters.

In the general case, the possible Young’s orbital diagram $\\{f\\}_{L}$ of
some nucleus $A(\\{f\\})$, consisting of two parts
$A_{1}(\\{f_{1}\\})+A_{2}(\\{f_{2}\\})$, is the direct outer product of
Young’s orbital diagrams $\\{f\\}_{L}=\\{f_{1}\\}_{L}\times\\{f_{2}\\}_{L}$
and is determined by Littlewood’s theorem [23, 24]. According to Elliot’s
theorem, each Young’s diagram is associated with a certain orbital angular
momentum or their combination.

Spin-isospin diagrams are the direct inner product of the spin and isospin
Young diagrams of a nucleus consisting of $A$ nucleons
$\\{f\\}_{ST}=\\{f\\}_{S}\otimes\\{f\\}_{T}$. For a system with no more than
eight particles such diagrams are provided in Table C of Ref. [25]. A detailed
procedure for defining the corresponding momenta can be found in the classical
monograph [26]. Let us note that in Ref. [26] the definition for inner and
outer products is reverses.

The total Young’s diagram of the nucleus is defined as the direct inner
product of the orbital and spin-isospin diagram
$\\{f\\}=\\{f\\}_{L}\otimes\\{f\\}_{ST}$. The total wave function of the
system under antisymmetrization does not vanish identically, only if it
contains an antisymmetric component $\\{1^{N}\\}$, where $N$ is the number of
nucleons. In this case the conjugates $\\{f\\}_{L}$ and $\\{f\\}_{ST}$ are
multiplied. Therefore, the diagrams $\\{f\\}_{L}$ conjugated to $\\{f\\}_{ST}$
are allowed in this channel. All other orbital symmetries are forbidden since
they lead to zero total wave function of the particle system after
antisymmetrization.

##  III The potentials construction within the MPCM

Let us describe in more detail the procedure for constructing the intercluster
partial potentials. Below we define the criteria and outline the sequence for
finding parameters for the potentials and indicating their errors as well as
ambiguities.

###  III.1 Discrete states

For the bound states of two clusters, the interaction potentials within the
framework of the MPCM are constructed based on the requirement imposed to
describe the main observable characteristics of such a nucleus. In this case,
the potential parameters are fixed. It should be noted that this requirement
is an idealized scenario that exists in the nucleus since it assumes that the
ground state (GS) is a two-body single channel with probability closed to
unity. First, we find the parameters of the bound state potentials. For the GS
with a given number of bound allowed and forbidden states in the partial wave,
these parameters are fixed unambiguously in terms of the binding energy, the
radius of the nucleus, and the AC. When constructing the partial interaction
potentials in the MPCM, it is assumed that interactions depend not only on the
orbital angular momentum $L$ but also on the total spin $S$ and the total
angular momentum $J$ of the system and also depend on Young’s orbital
diagrams. As in earlier work [28], we use Gaussian interaction potentials,
which depend on the quantum numbers $JLS$, and Young’s diagrams $\\{f\\}_{L}$.
Therefore, for different $JLS$, we have different values of the parameters of
the partial potentials.

The accuracy of determining the parameters of the BS potential is connected
directly with the accuracy of the AC. The potential does not contain any other
ambiguities, since according to Young’s diagrams, the classification of states
makes it possible to unambiguously fix the number of bound forbidden and
allowed states in a given partial wave. The number of bound states ultimately
determines the depth of the potential, while the width depends entirely on the
value of the AC. If one fixes two parameters of the potential using two
particular quantities — the binding energy and the AC — the error of the
binding energy is seen to be much less than that of the AC.

It should be noted that any calculations of the charge radius reflect the
errors of the underlying model. In any model, the magnitude of such a radius
depends on the integral of the model wave functions, thereby compounding
sources of error. At the same time, the values of AC are determined from the
asymptotic behavior of the model WFs at one point and contain significantly
less error. The potentials of the BSs are constructed to obtain the best
agreement with the values of the AC extracted independently from the
experimental data. For more details, see Ref. [27].

###  III.2 Continuum states

For the potentials of the continuous spectrum, the intercluster potential of
the nonresonant scattering process for a given number of allowed and forbidden
BSs in the considered partial wave is also constructed quite unambiguously
based on the scattering phase shifts. The accuracy of the potential parameters
sometimes as high as 20–30%, is associated with the precision of the extracted
scattering phase shifts from experimental data. For the 6Li$(p,\gamma)^{7}$Be
reaction, the potential is unambiguous since the classification, according to
Young’s diagrams, makes it possible to fix the number of bound states. This
completely determines the potential depth, and its width is determined by the
shape of the scattering phase shifts.

When constructing the nonresonant scattering potential based on the data for
the nuclear spectra, it is difficult to estimate the accuracy of the
parameters even for a given number of BSs. However, one can expect that it
will not exceed the error discussed above. This potential should lead to a
scattering phase shift close to zero or rise to a smoothly decreasing phase
shift at low energies, since there are no resonance levels in the spectra of
the nucleus.

In resonance scattering, when a relatively narrow resonance is present in the
partial wave at low energies for a given number of BSs, the potential is
constructed completely unambiguously. The accuracy of determining the
parameters of the interaction potentials is determined by the following
factors. The depth of the potential depends on the resonance energy $E_{x}$
and the number of BSs. The width is determined by the accuracy of the
experimental values of the level width $\Gamma$.

The error of the parameters, approximately 5–10%, usually does not exceed the
error of the energy level width. This also applies to the construction of the
partial potential from the resonant scattering phase shifts and the
determination of its parameters from the spectral resonance of the nucleus
[19, 28].

##  IV Cross-section and WFs asymptotics

To calculate the total cross-sections of radiative capture processes, we use
the well-known formula for the transitions of $NJ$ multipolarity [19, 28]

$\begin{gathered}\sigma(NJ,J_{f})=\frac{8\pi Ke^{2}}{\hbar^{2}\
k^{3}}\frac{\mu}{(2S_{1}+1)(2S_{2}+1)}\frac{J+1}{J[(2J+1)!!]^{2}}A_{J}^{2}(NJ,K)\times\\\
\times\sum\limits_{L_{i},J_{i}}P_{J}^{2}(NJ,J_{f},J_{i})I_{J}^{2}(k,J_{f},J_{i})\end{gathered}$
(2)

where the matrix elements of orbital $EJ(L)$-transitions have the following
form $(S=S_{i}=S_{f})$

$P_{J}^{2}(EJ,J_{f},J_{i})=\delta_{S_{i}S_{f}}[(2J+1)(2L_{i}+1)(2J_{i}+1)(2J_{f}+1)](L_{i}0J0|L_{f}0)^{2}\begin{Bmatrix}L_{i}&S&J_{i}\\\
J_{f}&J&L_{f}\end{Bmatrix},$ (3)
$A_{J}(EJ,K)=K^{J}\mu^{J}\left(\frac{Z_{1}}{m_{1}^{J}}+(-1)^{J}\frac{Z_{2}}{m_{2}^{J}}\right),\
\ \ \ I_{J}(k,J_{f},J_{i})=\langle\chi_{f}|r^{J}|\chi_{i}\rangle,$ (4)

and the matrix elements of the magnetic $M1(S)$-transition are written as
follows $(S=S_{i}=S_{f},L=L_{i}=L_{f})$

$P_{1}^{2}(M1,J_{f},J_{i})=\delta_{S_{i}S_{f}}\delta_{L_{i}L_{f}}[S(S+1)(2S+1)(2J_{i}+1)(2J_{f}+1)]\begin{Bmatrix}S&L&J_{i}\\\
J_{f}&1&S\end{Bmatrix},$ (5) $A_{1}(M1,K)=\frac{\hbar
K}{m_{0}c}\sqrt{3}\left(\muup_{1}\frac{m_{2}}{m_{1}+m_{2}}-\muup_{2}\frac{m_{1}}{m_{1}+m_{2}}\right),\
\ \ \ I_{J}(k,J_{f},J_{i})=\langle\chi_{f}|r^{J-1}|\chi_{i}\rangle,\ J=1.$ (6)

In Eqs. (2)–(6) $K=E_{\gamma}/\hbar c$ is the wave number of the emitted
photon with energy $E_{\gamma}$, $m_{1}$, $m_{2}$ and $\muup_{1}$, $\muup_{2}$
are the masses and magnetic momenta of the clusters, respectively, and $\mu$
is the reduced mass of the system. Namely, in the present calculations for the
reaction 6Li$(p,\gamma)^{7}$Be, $m_{1}\equiv m_{p}=1.00727646677$ amu,
$m_{2}\equiv m_{{}^{6}\text{Li}}=6.0151232$ amu [29] and
$\muup_{1}\equiv\muup_{p}=2.792847\muup_{0},$
$\muup_{2}\equiv\muup_{{}^{6}\text{Li}}=0.822\muup_{0}$ [30, 31], where
$\muup_{0}$ is nuclear magneton, $\hbar^{2}/m_{0}=41.4686$ MeV$\cdot$fm2,
where $m_{0}=931.494$ MeV is the atomic mass unit (amu).

The point-like Coulomb potential is of the form
$V_{\text{Coul}}\text{(MeV)}=1.439975Z_{1}Z_{2}/r$, where $r$ is the relative
distance between the particles of the channel in fm, and $Z_{1}$ and $Z_{2}$
are the charges in units of the elementary charge. The Coulomb parameter
$\eta=\mu Z_{1}Z_{2}e^{2}/k\hbar^{2}$ is represented in the form
$\eta=3.44476\cdot 10^{-2}\ Z_{1}Z_{2}\mu/k$, where $k$ is the wavenumber in
fm-1 and is determined by the energy $E_{c.m.}$ of the interacting particles,
$k^{2}=2\mu E_{c.m.}/\hbar^{2}$.

In our calculations, we use the dimensionless AC denoted as $C_{\text{w}}$
[32]

$\chi_{L}(r)=\sqrt{2k_{0}}\ C_{\text{w}}\ W_{-\eta L+1/2}(2k_{0}r).$ (7)

A dimensional asymptotic constant $C$ is related to the asymptotic
normalization coefficient (ANC) $A_{NC}$ by the expression [27]

$A_{NC}=\sqrt{S_{F}}C,$ (8)

where $S_{F}$ is the spectroscopic factor and $C$ is the dimensional AC that
can be represented using the asymptotics of the WF

$\chi_{L}(r)\xrightarrow[r\rightarrow R]{}CW_{-\eta
L+1/2}\left(2k_{0}r\right).$ (9)

In Eq. (9) $R$ is the large distance where the nuclear potential vanishes and
$\chi_{L}(r)$ is the wave function of the bound state obtained from the
solution of the radial Schr$\ddot{\text{o}}$dinger equation and normalized to
unity. The Whittaker function $W_{-\eta L+1/2}$ of the bound state determines
the asymptotic behaviour of the WF. The wave number $k_{0}$ is related to the
channel binding energy $E_{b}$.

For a continuous spectrum, the function $\chi_{i}$ found numerically is
matched to asymptotics $u_{L}(R)$ of the form

$N_{L}u_{L}(r)\xrightarrow[r\rightarrow
R]{}F_{L}(kr)+\tan(\deltaup_{S,L}^{J})G_{L}(kr).$ (10)

Here $F_{L}$ and $G_{L}$ are Coulomb regular and irregular functions [33].
They are the solutions of the Schr$\ddot{\text{o}}$dinger equation with the
Coulomb potential. $\deltaup_{S,L}^{J}$ are the scattering phase shifts
depending on the $JLS$ momenta of the system and $N_{L}$ is the normalizing
constant of the numerical radial function $u_{L}(R)$ for the continuum.

##  V Cluster states classification and the BS potentials

Consider the classification of the BSs of the $p^{6}$Li system according to
Young’s diagrams. In this case there is only one Young’s orbital diagram for
the 6Li nucleus $\\{42\\}$. It is believed that the system’s potentials are
dependent on the diagrams or combinations of these diagrams in various states.
Thus, if the orbital diagram $\\{42\\}$ allowed in the 2H4He-cluster channel
is accepted for the 6Li nucleus, then the $p^{6}$Li system with spin $S=1/2$
contains a forbidden level with diagram $\\{52\\}$ and orbital momenta of
$L=0,2$, and the allowed states with configurations $\\{43\\}$ for $L=1,3$ and
$\\{421\\}$ for $L=1,2$. Hence, the $p^{6}$Li potentials must have a forbidden
state related to $\\{52\\}$ in the $S$ wave. The allowed bound state
corresponds to the $P$ wave with the two Young’s diagrams $\\{43\\}$ and
$\\{421\\}$. In the quartet spin channel $S=3/2$ of the system, only one
diagram $\\{421\\}$ is allowed for $L=1,2$ [28]. Since there are two allowed
diagrams $\\{43\\}$ and $\\{421\\}$ in the doublet spin state of the $p^{6}$Li
system, the scattering states turn out to be mixed in orbital symmetries. At
the same time, only one allowed diagram $\\{43\\}$ usually corresponds to the
doublet ground state of the 7Be nucleus in the $p^{6}$Li channel with
$J^{\pi}=3/2^{-}$ and $L=1$.

Here, the $p^{6}$Li system is completely analogous to the $p^{2}$H channel in
the 3He nucleus. In the latter case the doublet state is also mixed according
to Young’s diagrams $\\{3\\}$ and $\\{21\\}$ [24]. Therefore, the potentials
constructed based on the elastic scattering phase shifts of the $p^{6}$Li
system or the level spectra of the 7Be nucleus cannot be used to describe the
GS of the 7Be nucleus in the $p^{6}$Li channel. Pure in orbital symmetry with
Young’s diagram $\\{43\\}$, the ${}^{2}P_{3/2}$ potential of the ground state
of 7Be reproduces the binding energy of the GS of the nucleus consistent with
the $p^{6}$Li system and its asymptotic constant.

The scattering potentials are constructed based on the spectra of the 7Be
nucleus from Ref. [30]; it has no major difference from the newer compilation
[34].

The orbital state’s classification of the $p^{6}$Li system is shown in Table
1.

Table 1: The classification of the orbital states of the $p^{6}$Li system [14, 19]. The following notations are used: $S$ and $L$ are spin and orbital angular momentum of the system, respectively, $\\{f\\}_{S}$, $\\{f\\}_{T}$ and $\\{f\\}_{ST}$ for isospin $T=1/2$, and $\\{f\\}_{L}$ are spin, isospin and spin-isospin, and possible orbital Young’s diagrams, and $\\{f\\}_{\text{AS}}$, $\\{f\\}_{\text{FS}}$ are Young’s diagrams of allowed and forbidden orbital states. $S$ | $\\{f\\}_{S}$ | $\\{f\\}_{T}$ | $\\{f\\}_{ST}$ | $\\{f\\}_{L}$ | $L$ | $\\{f\\}_{\text{AS}}$ | $\\{f\\}_{\text{FS}}$
---|---|---|---|---|---|---|---
$1/2$ | {43} | {43} | {7}+{61}+{52}+{511}+{43}+ | {52} | $0,2$ | — | {52}
+{421}+{331}+{4111}+ | {43} | $1,3$ | {43} | —
+{322}+{3211}+{2221} | {421} | $1,2$ | {421} | —
$3/2$ | {52} | {43} | {61}+{52}+{511}+ | {52} | $0,2$ | — | {52}
+{43}+2{431}+{331}+ | {43} | $1,3$ | — | {43}
+{322}+{3211} | {421} | $1,2$ | {421} | —

We use a Gaussian potential that depends on the momenta of the system [28] and
Young’s diagrams

$V(r,JLS,\\{f_{L}\\})=-V_{0}(JLS,\\{f_{L}\\})\exp(-\alpha_{JLS,\\{f_{L}\\}}r^{2}).$
(11)

In Eq. (11) $V_{0}$ is the potential depth and $\alpha$ is related to the
potential width. The choice of parameters for the bound and scattering states
is discussed in detail in Sec. VI.

A compilation of the AC data for the ground and first excited state in the
$p^{6}$Li channel of 7Be is presented in Table 2, with
$C_{\text{w}}^{2}=A_{NC}^{2}/2k_{0}S_{F}$. Here $\sqrt{2k_{0}}=0.983$ fm-1/2
for the GS and $\sqrt{2k_{0}}=0.963$ fm-1/2 for the FES.

Table 2: AC data for the ground and the first excited states in $p^{6}$Li channel of 7Be. BS |  Reference |  $A_{NC}$, fm-1/2 |  $S_{F}$ |  $C_{\text{w}}$
---|---|---|---|---
GS |  Nollett and Wiringa [35], 2011 |  2.85(3) |  1 |  2.90(3)
|  Huang et al. [13], 2010 |  2.01 |  0.66 – 1 |  2.28(24)
|  Timofeyuk [36], 2013 |  1.80 |  0.46 – 0.87 |  2.32(37)
|  Burtebayev et al. [37], 2013 |  1.77(8) |  0.55 – 0.81 |  2.23(31)
|  Gnech and Marcucci [18], 2019 |  2.654 |  1.003 |  2.65
|  Kiss et al. [11], 2021 |  2.19(9) |  0.98(30) |  2.35(43)
FES |  Huang et al. [13], 2010 |  1.91 |  0.66 – 1.02 |  2.20(24)
|  Timofeyuk [36], 2013 |  1.91 |  0.62 – 1.21 |  2.17(36)
|  Burtebayev et al. [37], 2013 |  1.95(9) |  0.85 – 1.03 |  2.10(20)
|  Gnech and Marcucci [18], 2019 |  2.528 |  1.131 |  2.53
|  Kiss et al. [11], 2021 |  2.18(6) |  1.08(32) |  2.26(39)

Summarizing the data in Table 2, we conclude that all cited data on the AC are
overlapped. While constructing the corresponding GS and FES potentials, we
used the average values indicated in bold font in Table 2. Just at the end of
our calculation story, a publication by Kiss et al. [11] appeared, and it
happened that these latest experimental results turned to be within the
defined intervals for ANC given above.

We use GS and FES in the form of only doublet ${}^{2}P_{3/2}$ and
${}^{2}P_{1/2}$ states, but we take the experimental data on ANC from [37]
since it is assumed that these states result in the observed ANC values. We do
not consider these states as a mix of doublet and quartet states, for
instance, ${}^{2+4}P_{3/2}$, is a prime example, as the quartet channel is not
allowed for the orbital Young’s diagram $\\{43\\}$ in Table 1.

All AC values are used here as a framework to obtain the parameters of the
$p^{6}$Li interaction BSs potentials. These potentials correspond to the
lower, upper, and average values of AC and accurately reproduce the binding
energies [30] of the bound states. The parameters of the potentials are
presented in Table 3.

Table 3: Parameters of $p^{6}$Li system bound state potentials and bound states characteristics. $E_{x}$, $E_{b}$, and $V_{0}$ are provided in MeV. The $R_{\text{ch}}$ and $R_{\text{m}}$ are given in fm. $\\#$ |  BS |  $E_{x}$ |  $J^{\pi}$ |  $E_{b}$ |  ${}^{2S+1}L_{J}$ |  $V_{0}$ |  $\alpha$, fm-2 |  $C_{\text{w}}$ |  $R_{\text{ch}}$ |  $R_{\text{m}}$
---|---|---|---|---|---|---|---|---|---|---
1 |  GS |  0 |  3/2- |  –5.60580 |  ${}^{2}P_{3/2}$ |  $100.750920$ |  $0.25$ |  $1.75(1)$ |  $2.49$ |  $2.51$
2 |  GS |  0 |  3/2- |  –5.60580 |  ${}^{2}P_{3/2}$ |  $74.504070$ |  $0.17$ |  $2.26(1)$ |  $2.58$ |  $2.58$
3 |  GS |  0 |  3/2- |  –5.60580 |  ${}^{2}P_{3/2}$ |  $60.998575$ |  $0.13$ |  $2.74(1)$ |  $2.64$ |  $2.61$
4 |  FES |  0.4291 |  1/2- |  –5.17670 |  ${}^{2}P_{1/2}$ |  $99.473500$ |  $0.25$ |  $1.68(1)$ |  $2.52$ |  $2.54$
5 |  FES |  0.4291 |  1/2- |  –5.17670 |  ${}^{2}P_{1/2}$ |  $73.333835$ |  $0.17$ |  $2.16(1)$ |  $2.59$ |  $2.59$
6 |  FES |  0.4291 |  1/2- |  –5.17670 |  ${}^{2}P_{1/2}$ |  $59.898120$ |  $0.13$ |  $2.61(1)$ |  $2.65$ |  $2.62$

##  VI $\boldmath p^{6}\text{Li}$ channel in the continuous spectrum

We usually assume that the scattering potentials can lead to the FSs [28], and
if the FSs are absent, the potential depth can be set to zero. The present
case refers to the scattering $P$ potentials without the FS, while the $S$ and
$D$ potentials have the bound forbidden state and, even at zero phase shifts,
must have a nonzero depth.

In the presence of one BS, the phase shift starts at $180^{\circ}$ [24], as
shown in Fig. 1 for the $S$ scattering phase shifts. Furthermore, we consider
$E2$ transitions from resonant $F$ scattering states with phase shifts shown
in Fig. 2. The parameters of potentials for all scattering processes for
transitions to GS and FES are given in Tables 4 and 5, respectively.

Table 4: The spectrum of 7Be levels [30] and scattering states in the $p^{6}$Li channel for the capture to the ${}^{2}P_{3/2}$ GS at a binding energy of 5.6058 MeV, along with $P^{2}_{J}$ from expressions (3) and (5). $E_{x}$, $E_{\text{res}}$, $\Gamma_{c.m.}$ and $V_{0}$ are provided in MeV. # |  $E_{x}$, | $J^{\pi}$ |  $E_{\text{res}}$, |  $\Gamma_{c.m.}$, | GS transition: | $P^{2}_{J}$ | $V_{0}$ |  $\alpha$, | $E_{\text{res}}$ | $\Gamma_{c.m.}$
---|---|---|---|---|---|---|---|---|---|---
exp. | exp. | exp. | $\left[{}^{2S+1}L_{J}\right]_{i}\rightarrow\left[{}^{2S+1}L_{J}\right]_{f}$ |  fm-2 | theory | theory
1 | No res. | $5/2^{+}$ | — | — | $E1:\ ^{2}D_{5/2}\rightarrow\ ^{2}P_{3/2}$ | $36/5$ | $58.0$ | $0.4$ | — | —
2 | No res. | $3/2^{+}$ | — | — | $E1:\ ^{2}D_{3/2}\rightarrow\ ^{2}P_{3/2}$ | $4/5$ | $58.0$ | $0.4$ | — | —
3 | No res. | $1/2^{+}$ | — | — | $E1:\ ^{2}D_{5/2}\rightarrow\ ^{2}P_{3/2}$ | $4$ | $58.0$ | $0.4$ | — | —
4 | No res. | $1/2^{-}$ | — | — | $M1:\ ^{2}P_{1/2}\rightarrow\ ^{2}P_{3/2}$ | $4/3$ | $0.0$ | $1.0$ | — | —
5 | $7.2(1)$ | $5/2^{-}$ | $1.59(10)$ | $0.40(5)$ | $E2:\ ^{2}F_{5/2}\rightarrow\ ^{2}P_{3/2}$ | $12/7$ | $111.60$ | $0.1$ | $1.60(1)$ | $0.62(1)$
6 | $9.29(31)$ | $7/2^{-}$ | $3.68(31)$ | $1.93(96)$ | $E2:\ ^{2}F_{7/2}\rightarrow\ ^{2}P_{3/2}$ | $72/7$ | $44.34$ | $0.05$ | $3.68(1)$ | $1.50(1)$
7 | $9.9$ | $3/2^{-}$ | $4.3$ | $1.8$ | $M1:\ ^{2}P_{3/2}\rightarrow\ ^{2}P_{3/2}$ | $5/3$ | $432.0$ | $1.5$ | $4.30(1)$ | $1.80(2)$
Table 5: The spectrum of 7Be energy levels [30] and scattering states in the $p^{6}$Li channel for the proton capture to the ${}^{2}P_{1/2}$ FES at a binding energy of 5.1767 MeV, along with $P^{2}_{J}$ from expressions (3) and (5). $E_{\text{x}}$, $E_{\text{res}}$, $\Gamma_{c.m.}$ and $V_{0}$ are provided in MeV. # |  $E_{\text{x}}$, | $J^{\pi}$ |  $E_{\text{res}}$, |  $\Gamma_{c.m.}$, | FES transition: | $P^{2}_{J}$ | $V_{0}$ |  $\alpha$, | $E_{\text{res}}$ | $\Gamma_{c.m.}$
---|---|---|---|---|---|---|---|---|---|---
exp. | exp. | exp. | $\left[{}^{2S+1}L_{J}\right]_{i}\rightarrow\left[{}^{2S+1}L_{J}\right]_{f}$ | fm-2 | theory | theory
1 | No res. | $3/2^{+}$ | — | — | $E1:\ ^{2}D_{3/2}\rightarrow\ ^{2}P_{1/2}$ | $4$ | $58.0$ | $0.4$ | — | —
2 | No res. | $1/2^{+}$ | — | — | $E1:\ ^{2}S_{1/2}\rightarrow\ ^{2}P_{1/2}$ | $2$ | $58.0$ | $0.4$ | — | —
3 | No res. | $1/2^{-}$ | — | — | $M1:\ ^{2}P_{1/2}\rightarrow\ ^{2}P_{1/2}$ | $1/6$ | $0.0$ | $1.0$ | — | —
4 | $7.2(1)$ | $5/2^{-}$ | $1.59(10)$ | $0.40(5)$ | $E2:\ ^{2}F_{5/2}\rightarrow\ ^{2}P_{1/2}$ | $6$ | $111.6$ | $0.1$ | $1.60(1)$ | $0.62(1)$
5 | $9.9$ | $3/2^{-}$ | $4.3$ | $1.8$ | $M1:\ ^{2}P_{3/2}\rightarrow\ ^{2}P_{1/2}$ | $4/3$ | $432.0$ | $1.5$ | $4.30(1)$ | $1.80(2)$
Figure 1: Doublet and quartet $S$ phase shifts of elastic $p^{6}$Li scattering
at low energies. The ${}^{2}S$ and ${}^{4}S$ phase shifts are taken from Ref.
[38] and shown by $\medblackcircle$ and $\medblacktriangleup$, respectively.
Results from [14, 15] obtained according to [38] are shown by the dashed-
dotted curves. Results of the present work are shown by the solid curve.

In addition, we consider the resonance at an excitation energy of 9.9 MeV [30]
according to Fig. 3 (4.3 MeV above the threshold) in a ${}^{2}P_{3/2}$
scattering state of width 1.8 MeV in c.m. Considering such an $M1$ transition
to the ${}^{2}P_{3/2}$ ground state or an $M1$ transition from the
${}^{2}P_{1/2}$ scattering state to the ${}^{2}P_{1/2}$ FES are possible due
to the presence of different Young’s diagrams in the bound and scattering
states. Recall that the BSs have the diagram $\\{43\\}$, and the scattering
states are mixed according to the two diagrams $\\{43\\}+\\{421\\}$ [15].

Table 4 shows possible transitions to the 7Be nucleus GS from various
$p^{6}$Li scattering states with ${}^{2S+1}L_{J}$. The possible transitions to
the FES from different scattering states with ${}^{2S+1}L_{J}$ are shown in
Table 5. The resonance energies and widths are obtained with the corresponding
parameters of the scattering potentials. For the $P_{1/2}$ scattering wave,
zero-depth potentials are used since the scattering $P$ waves do not contain
forbidden BSs. For the ${}^{2}D$ wave potentials, ${}^{2}S$ wave parameters
are used for $L=2$.

Figure 2: $P$ and $F$ phase shifts of elastic $p+^{6}$Li scattering obtained
for scattering potentials with the parameters from Tables 4 and 5. Figure 3:
Schematics of the energy spectrum of 7Be. The energy are given in MeV and
figure is not drown to scale [30]. a the above-threshold resonance at 6.73 MeV
with a width of 1.2 MeV refers to the 4He3He channel that is not considered in
the present work.

Resonant phase shifts of elastic scattering $p+^{6}$Li are shown in Fig. 2.
The above-threshold resonance at 6.73 MeV with a width of 1.2 MeV indicated in
Fig. 3 refers to the 4He3He channel [30] and is not considered in our previous
works [14, 15]. Note again that in the MPCM we used, Young’s orbital diagram
{43} is forbidden in the quartet state, as shown in Table 1, and this
particular diagram corresponds to the GS of the 7Li nucleus. Therefore, in the
GS there is only a doublet ${}^{2}P_{3/2}$ state (without impurity of
${}^{4}P_{3/2}$), which is allowed for the diagram {43}. Thus, our model [14,
15] predicted the absence of resonance at 6.73 MeV with $J^{\pi}=5/2$ in the
nucleon channel, or, in other words, the impossibility of the $M1$ transition
from this resonance to the GS. This has been confirmed by the new LUNA results
[9] and, indirectly, by the data of [30]. The width of the resonance peak at
9.29 MeV is taken from Table 7.10 of [30], although another state,
${}^{2}P_{1/2}$, is indicated therein. At 9.27 MeV, the given moment is $7/2$,
as per Table 7.7 in Ref. [30], so we infer the presence of an $F$ state.
However, this resonance leads to a minimal increase in cross-sections at the
$E2$ transition. It is negligible against the background of the $E1$ resonance
at 4.3 MeV with a transition from the ${}^{2}P_{3/2}$ scattering state for the
potential parameters #7 (Table 4) or #5 (Table 5), respectively.

In Fig. 1 are shown the doublet and quartet $S$ phase shifts of elastic
$p^{6}$Li scattering at low energies. The ${}^{2}S$ potential from Table 4
with a depth of 58 MeV has one FS and allows one to describe the phase shifts
of [14] up to 1 MeV, shown in Fig. 1 by solid circles. Moreover, it gives
phase shifts below 2 MeV that coincide with the phase shifts obtained with the
potential from Ref. [14], with a depth of 126 MeV and a width of 0.15 fm-2.
This early potential has two FSs and does not agree with our new
classification from Table 1. To compare the results, we construct a new
potential (depth 58 MeV) that gives the most overlap in phase shifts from
prior work [14]. The phase shift of the new potential is given in Fig. 1 by a
solid curve, while dash-dotted curves refers to results from Ref. [14].

##  VII Astrophysical S-factor and reaction rate

A special feature of cross-sections of nuclear reactions with charged
particles at low and ultra-low energies is an extreme reduction by several
orders of a cross-section magnitude due to the decrease in transmission
probability through the Coulomb barrier. For practical purposes, the
astrophysical S-factor is introduced as

$S(E)=\dfrac{\sigma(E)}{P(E)}E,\ \ \ P(E)=e^{-2\pi\eta},$ (12)

where the factor $P(E)$ reflects the permeability of the Coulomb barrier. The
advantage of the astrophysical $S-$factor is that it shows a smooth energy
dependence at low energies.

Following the excellent manuscript by Christian Iliadis [39], we would like to
provide a brief discussion of the use of the $S-$factor in conventional
calculation schemes in order to clarify the current approach. The definition
of the $S-$factor (12) allows to write the following expression for the
reaction rate:

$N_{A}\langle\sigma\nu\rangle=\left(\dfrac{8}{\pi\mu}\right)^{1/2}N_{A}\left(\kappa
T_{9}\right)^{-3/2}\int\limits_{0}^{\infty}e^{-2\pi\eta}S(E)e^{-E/\kappa
T_{9}}dE.$ (13)

In Eq.(13) $\kappa$ is the Boltzmann constant and $N_{A}$ is Avogadro’s
number. A notable effort has been expended to bring the integral in (13) to
analytical form. This is possible only if the expansion of the $S(E)-$factor
in the $E$ series, given by

$S(E)\approx S(0)+S^{\prime}(0)E+S^{\prime\prime}E^{2}.$ (14)

is valid at low energies (c.f. Ref. [39], Section 3.2). Expression (14)
explains the active interest to determine the value $S(0)$ as a key one for
the calculation of the reaction rate in the form of (13) with its further
analytical parameterizations.

Table 6 presents the available experimental data of the astrophysical
$S-$factor for the 6Li$(p,\gamma)^{7}$Be reaction, as well as the extrapolated
$S(0)$ values. Results of previous theoretical calculations and the present
work are presented in Table 7. Our calculations are analyzed and interpreted
based on the experimental data [9] as one of the newest and most accurate.

In the present work we introduce some corrections to [14, 15] that allow us to
extend the energy interval for the cross-sections and corresponding
$S-$factors. It is also worth mentioning that in Ref. [18], the authors used
the calculation scheme based on [15]. The value $S(0)=95.0$ eV$\cdot$b is
obtained which is consistent with the LUNA experimental data [9]. However, the
astrophysical reaction rate is missing in Ref. [18].

Table 6: Experimental data on the astrophysical $S-$factor of6Li$(p,\gamma)^{7}$Be. $S(E)$ and $S(0)$ are given in eV$\cdot$b and $S(0)$ are extrapolated values. $E$, keV | $S(E)$ | $S(0)$ | Reference |  Method/project
---|---|---|---|---
$135$ | $51(15)$ | — | Switkowski et al. [4], 1979 |  $\gamma$-ray Ge(Li) spectrometers for proton bombarding energies 200 – 1200 keV
$340$ | $43(3)$ | — | Ostojic et al. [5], 1983 |  Direct radiative capture
$40$ | $65$ | $65$ | Cecil et al. [40], 1992 |  Thick-target $\gamma$-ray-to-charged-particle branching ratio measurements
$35$ | $40(14)$ | — | Bruss [6], 1993 |
— | — | $79(18)$ | Prior et al. [10], 2004 |  Polarized proton beams, TUNL
$250$ | $95(10)$ | — | He et al. [8], 2013 |  320 keV platform with highly charged ions
$60$ | $92(6)$ | $95(9)$ | Piatti et al. [9], 2020 |  LUNA collaboration
Table 7: Theoretical calculations of the astrophysical $S-$factor of 6Li$(p,\gamma)^{7}$Be. $S(0)$, eV$\cdot$b | Reference | Model
---|---|---
$106$ | Barker [41], 1980 | Direct-capture potential model
$105$ | Arai et al. [12], 2002 | Four-cluster microscopic model
$95,5$ | Huang et al. [13], 2010 | Single-particle model
$114$ | Dubovichenko et al. [14, 15], 2010, 2011 | Modified potential cluster model
$73^{+56}_{-11}$ | Xu et al. [16], 2013 | Direct-capture potential model
$88,34$ | Dong et al. [17], 2017 | Gamow shell model
$103,9$ | Gnech and Marcucci [18], 2019 | Potential cluster model
$96,5\pm 5,7$ | Kiss et al. [11], 2021 | Modified two-body potential method
$92\pm 12$ | Kiss et al. [11], 2021 | Modified two-body potential method
$98,3$ | Present work | Modified potential cluster model

The results of the present calculations of $S-$factors along with available
experimental data are shown in Figs. 4, 5, and 6. Fig. 4 shows the
astrophysical $S(E)-$factor of 6Li$(p,\gamma_{0})^{7}$Be capture to the GS of
the 7Be nucleus in an energy range up to 5 MeV. The solid red curve 2 and the
two dashed curves, blue 1 and green 3, show the calculation for all
transitions to the GS given in Table 4. Parameters $V_{0}$ (only integer
values) and $\alpha$ of the corresponding potentials as well as $C_{\text{w}}$
from Table 3 are indicated in the figures, and the parameters of the GS
potentials are taken from Table 4. The solid red curve 2 is the result for the
potential with the set of parameters #2 from Table 3, leading to the average
value of AC.

Figure 4: Astrophysical $S-$factor of the 6Li$(p,\gamma_{0})^{7}$Be capture to
GS. Experimental data are taken from $\medblacktriangleup$ – [4],
$\medblackcircle$ – [8], $\medblacksquare$ – [9], $\medblackdiamond$ – [7],
$\medblackstar$ – [42], $\medsquare$ – [3], $\medtriangledown$ – [5].
Parameters of the continuum potential are given in Table 4. A band shows the
sensitivity to changes in $C_{\text{w}}$.

In Fig. 5, similar curves show the results for transitions and potentials from
Table 5 to FES. FES potentials have three sets of parameters from Table 3. The
solid red curve 2 shows the results for capture with set # 5 from Table 3,
allowing us to determine the average value of AC. This result is in good
agreement with the experimental data [6] presented in Fig. 5. The two dashed
curves 1 and 3 almost completely cover the interval or band of cross-section
errors of the capture to the FES.

Figure 5: Astrophysical $S-$factor of the 6Li$(p,\gamma_{1})^{7}$Be capture to
FES. Experimental data $\medcircle$ are taken from [6]. Parameters of the
continuum potential are given in Table 5. A band shows the sensitivity to
changes in $C_{\text{w}}$.

In Fig. 6, similar curves show the astrophysical $S-$factor for the total
cross-sections corresponding to the transition to GS and FES. The two dashed
curves 1 and 3 show the range of $S-$factor values due to ambiguities in the
AC of the GS and FES. For the scattering potentials, the parameters from
Tables 4 and 5 are used.

Figure 6: Astrophysical $S-$factor of the 6Li$(p,\gamma_{0+1})^{7}$Be capture
to GS and FES. Experimental data for capture are from $\medblacktriangleup$ –
[4], $\medblackcircle$ – [8], $\medblacksquare$ – [9], $\medblackdiamond$ –
[7], $\medsquare$ – [3], $\medtriangledown$ – [5], $\medblackstar$ – [42],
$\medcircle$ – [6]. The parameters of GS and FES potentials are listed in the
figure. A band shows the sensitivity to changes in $C_{\text{w}}$.

The best agreement of the $S-$factor with experimental data is achieved for
the values of $C_{\text{w}}=2.74$ for the GS and $C_{\text{w}}=1.68$ for the
FES. We recommend these values as the most reliable benchmarks for future
experimental studies.

Fig. 6 shows that almost all experimental data lie between the solid red curve
2 and the green dashed curve 3. If we use the GS and FES potentials set of
parameters #3 and #4 from Table 3, the result is shown in Fig. 6 by the black
curve 4. In this case the experimental data [9] are reproduced entirely, and
$S-$factor at 10 keV is found to be 101 eV$\cdot$b. For the scattering
potentials, the data from Table 4 and Table 5 are used.

Due to the uncertainty of the $S-$factor that arises from the uncertainty of
the AC, it is desirable to select other options for the potentials of the GS
and FES to correctly describe the LUNA data [9]. This can be the subject of
future work if more accurate data is compiled for the AC of the ground and
first excited states of the 7Be nucleus.

The approximation of the $S-$factor shown by the black curve 4 in Fig. 6 has
an analytical form

$S(E)=S_{0}+S_{1}E+S_{2}E^{2}$ (15)

with parameters $S_{0}=98.31$ eV$\cdot$b, $S_{1}=-187.18$
MeV${}^{-1}\cdot$eV$\cdot$b and $S_{2}=442.51$ MeV${}^{-2}\cdot$eV$\cdot$b.
This approximation leads to $\chi^{2}=2.4\cdot 10^{-4}$ with the error of 5%
in the energy range of 30 to 100 keV. This shows that $S(0)=98.3$ eV$\cdot$b
and $S(30)=93.1$ eV$\cdot$b.

New experimental data from LUNA [9] can be approximated to the first order

$S(E)=S_{0}+S_{1}E$ (16)

with parameters $S_{0}=91.952$ eV$\cdot$b and $S_{1}=-75.471$
MeV${}^{-1}\cdot$eV$\cdot$b, leading to $\chi^{2}=0.6$ and $S(0)=92$
eV$\cdot$b.

To compare the calculated $S-$factor at zero energy (10 keV), we present the
known results for the total $S(0)$: $79(18)$ eV$\cdot$b [10], 105 eV$\cdot$b
(at 10 keV) [12] and 106 eV$\cdot$b [41].The $S-$factor for transitions to the
ground state in [40], 39 eV$\cdot$b is specified, and for the transition to
the first excited state, the $S-$factor value is equal to 26 eV$\cdot$b, the
total $S-$factor is 65 eV$\cdot$b. In our previous works [14, 15], a value of
114 eV$\cdot$b was obtained. The summary for $S-$factor experimental and
theoretical values are presented in Tables 6 and 7.

To sum up, our astrophysical $S-$factor is given in Fig. 7 with a solid red
curve, together with experimental data and theoretical calculations. The
$R-$matrix fit of the data from LUNA collaboration [9] and Switkowski et al.
[4] is represented with the solid blue curve. A solid green curve was obtained
by Kiss et al. [11] using the weighted means of the ANCs from the analysis of
the 6Li(3He,$d$)7Be transfer reaction within the modified two-body potential
method (MTBPM). In addition, [11] contains the results for the $S-$factor of
the 6Li$(p,\gamma)^{7}$Be reaction calculated within the MTBPM, using the
values of ANCs obtained from the analysis of the experimental astrophysical
$S-$factors of the 6Li$(p,\gamma)^{7}$Be reaction [9]. These results are given
in Fig. 7 with the solid black curve.

Figure 7: Comparison of 6Li$(p,\gamma)^{7}$Be reaction astrophysical
$S-$factors. Experimental data for capture to GS and FES are from
$\medblacktriangleup$ – [4], $\medblackcircle$ – [8], $\medblacksquare$ – [9],
$\medblackdiamond$ – [7], $\medsquare$ – [3], $\medtriangledown$ – [5],
$\medblackstar$ – [42]. Results of calculations: red curve – present work;
blue curve – Ref. [9]; black and green curves – Ref. [11].

To calculate the 6Li$(p,\gamma)^{7}$Be capture reaction rate in units of cm
3mol-1s-1, we used the expression [43] analogous to Eq. (13), but substituting
the corresponding constants values

$N_{A}\langle\sigma\nu\rangle=3.7313\cdot
10^{4}\mu^{-1/2}T_{9}^{-3/2}\int\limits_{0}^{\infty}\sigma(E)E\exp(-11.605E/T_{9})dE.$
(17a) In Eq. (17a) $E$ is given in MeV, the total cross-section $\sigma(E)$ is
taken in $\muup$b, and $\mu$ is the reduced mass in amu and $T_{9}=10^{9}$ K
[43]. Using real integration limits $E_{min}$ and $E_{max}$ Eq. (17a) becomes
$N_{A}\langle\sigma\nu\rangle=3.7313\cdot
10^{4}\mu^{-1/2}T_{9}^{-3/2}\int\limits_{E_{min}}^{E_{max}}\sigma(E)E\exp(-11.605E/T_{9})dE.$
(17b)

It is important to stress this fact as the choice of $E_{max}$ in Eq. (17b)
may have a significant impact on the final result for the reaction rate. The
reaction rate (17b) is calculated based on cross-sections, displayed in the
form of $S-$factors (12) in Figs. 4, 5 and 6 within the energy range
$E_{min}=1$ keV to $E_{max}=5$ MeV. The results of these calculations are
plotted in Fig. 8.

Figure 8: Total 6Li$(p,\gamma)^{7}$Be capture reaction rate. Curves indicate
the different sums of the capture rates to the GS and FES. Curves are
designated as in Fig. 6. The yellow band is taken from Ref. [16].

As noted above, curve 4 in Fig. 6 is in best agreement with all experimental
data of the $S-$factor. Therefore, the corresponding reaction rate, also
marked as curve 4 in Fig. 7 is most recommended description of the reaction
rate. Curve 4 can be approximated by a function of the form [44]

$\begin{gathered}N_{A}\langle\sigma\nu\rangle=a_{1}/T_{9}^{1/3}\exp\left(-a_{2}/T_{9}^{2/3}\right)\left(1+a_{3}T_{9}^{1/3}+a_{4}T_{9}^{2/3}+a_{5}T_{9}+a_{6}T_{9}^{4/3}+a_{7}T_{9}^{5/3}\right)+\\\
+a_{8}T_{9}^{2/3}\exp\left(-a_{9}/T_{9}^{1/3}\right).\end{gathered}$ (18)

The parameters of approximation (18) with an average value of $\chi^{2}=0.014$
and the error of 5% are given in Table 8.

Table 8: The reaction rate approximation parameters. $i$ | $a_{i}$ | $i$ | $a_{i}$ | $i$ | $a_{i}$
---|---|---|---|---|---
$1$ | $0.00319$ | $4$ | $-544516.9$ | $7$ | $85197.34$
$2$ | $4.16292$ | $5$ | $16401.8$ | $8$ | $924167.3$
$3$ | $3721.884$ | $6$ | $-8044.932$ | $9$ | $8.36494$

Comparing our results with the reaction rate presented in NACRE [43] and NACRE
II [16], we follow the format of Fig. 4 from [9]. We added our results for the
reaction rate, normalized to the reaction rate from NACRE [43]. The comparison
is shown in Figs. 10 and 10.

Figure 9: Comparison of the astrophysical reaction rates in the range 0.001 to
1 $T_{9}$ from [43, 16, 9] to present work, normalized to the NACRE rate [43].
Dotted curves represent the uncertainty of the NACRE [43] rate, while shaded
areas represent the uncertainties from LUNA [9] (pale red) and NACRE II [16]
(yellow).

Figure 10: Comparison of the astrophysical reaction rates in the range 0.001
to 10 $T_{9}$ from [43, 16, 9] to present work, normalized to the NACRE rate
[43]. Dotted curves represent the uncertainty of the NACRE [43] rate, while
shaded areas represent the uncertainties from LUNA [9] (pale red) and NACRE II
[16] (yellow).

As stated in Ref. [9], the LUNA ”thermonuclear reaction rate is 9% lower than
NACRE [43] and 33% higher than reported in NACRE II [16] at 2 MK, and the
reaction rate uncertainty has been significantly reduced”. Fig. 10 shows that
the deviation between the adopted reaction rate obtained in [9] and the
present calculations in the range of 0.01 to 1 $T_{9}$ does not exceed 5%.
Therefore, the present calculations confirm the above conclusion by Piatti et
al. [9].

Fig. 10 shows two of our results — the blue curve 1 is the reaction rate
calculated on the basis of the total cross sections in the energy range from 1
keV to 5 MeV, and the green one 2 shows the rate at the upper range of the
integration limit $E_{max}=0.6$ MeV, and this value corresponded to the upper
limit of energy in work [9]. The difference between these curves illustrates
the importance of the resonance region contribution for the cross sections in
the range of 0.6 MeV to 5 MeV at temperatures above 1 $T_{9}$.

##  VIII Conclusion

We present the results of calculations and analyses of the $S-$factor and
astrophysical reaction rate for the 6Li$(p,\gamma)^{7}$Be reaction in the
framework of MPCM. It is demonstrated that the MPCM approach has only one
ambiguity arising from the accuracy of the experimentally determined
asymptotic constants. This effect manifests as bands in Figs. 4 – 6 for the
astrophysical $S-$factor. Precise LUNA experimental data played a role of the
criterion, in reducing the ANC ambiguity with theoretical simulations.

Comparing the $R-$matrix method which is constrained by the parameterization
of the experimental cross-sections data, MPCM enables to implement
calculations in wider energy ranges. We extended the energy interval for the
total cross-sections and $S-$factors up to 5 MeV, including resonances in the
continuum. The numerical signature of this extension is seen in Fig. 10 for
the reaction rate.

It was also shown in the present work that MPCM had predicted the absence of
resonance at 6.73 MeV in the nucleon channel [14, 15], which was confirmed by
the LUNA results [9], as well as, indirectly, by the data from [30].

We suggest that the NACRE [43] and NACRE II [16] databases should be updated
in light of LUNA data [9] and present calculations.

## Dedication

We dedicate this paper to the memory of our colleague, Dr. Albert Dzhazairov-
Kakhramanov, who recently passed away from COVID19.

## Acknowledgments

This work was supported by the grant of the Ministry of Education and Science
of the Republic of Kazakhstan #AP08855556 ”Study of additional thermonuclear
reactions flowing in the process of controlled thermonuclear fusion on lithium
isotopes” through the V.G. Fesenkov Astrophysical Institute of the ”National
Center of Space Research and Technology” of the Aerospace committee of the
Ministry of Digital Development, Innovations and Aerospace Industry of the
Republic of Kazakhstan.

## Appendix A

A solution of a two-body problem for a discrete energy spectrum with a given
potential requires finding the binding energy of the system and the wave
function of the state. This problem can be solved using the Variational Method
and the Finite Difference Method (FDM) [45]. If both methods are used for the
same system of particles, it is possible to control the correctness of the
search for the binding energy and WF of the state. We already have used such
an approach for $p^{2}$H and $p^{3}$H systems in [46, 28] and demonstrated
that the FDM provides more precise description of the systems. Below we
present the FDM approach.

The calculation of the binding energy of a two-cluster system by the FDM
relies on the representation of the Schr$\ddot{\text{o}}$dinger equation in
finite differences [47]. The radial equation for the central potential [45]

$u^{\prime\prime}_{L}(r)+\left[k^{2}-V(r)\right]u_{L}(r)=0$ (19)

with some boundary condition for $k^{2}<0$, ($k^{2}=2\mu E/\hbar$) takes the
form of a Sturm-Liouville type boundary value problem. Recasting the second
derivative in finite difference form, we obtain

$\begin{gathered}u^{\prime\prime}=\left[u_{n+1}-2u_{n}+u_{n-1}\right]/h^{2},\
\ \ u_{n}=u(r_{n})\end{gathered}$ (20)

and (19) becomes a closed system of linear algebraic equations. Thus, for a
certain $k_{0}$, $D_{N}(k)=0$

$D_{N}(k)=\begin{pmatrix}\thetaup_{1}&1&0&.&.&.&0\\\
\alphaup_{2}&\thetaup_{2}&1&0&.&.&0\\\ 0&\alphaup_{3}&\thetaup_{3}&1&0&.&0\\\
.&.&.&.&.&.&.\\\ .&.&.&.&.&.&.\\\ 0&.&0&0&\alphaup_{N-1}&\thetaup_{N-1}&1\\\
0&.&0&0&0&\alphaup_{N}&\thetaup_{N}\end{pmatrix}=0.$ (21)

Eq. (21) allows one to determine the binding energy $E_{b}$ of a system of two
particles. The elements of the tridiagonal determinant (21) are defined as
follows:

$\displaystyle\alphaup_{n}=1,$
$\displaystyle\thetaup_{n}=k^{2}h^{2}-2-V_{n}h^{2},\ \ \ n=1,2,\dots,N-1,$
(22) $\displaystyle\alphaup_{N}=2,$
$\displaystyle\thetaup_{N}=k^{2}h^{2}-2-V_{N}h^{2}+2hf(\eta,L,Z_{N}),$
$\displaystyle Z_{n}=2kr_{n},$ $\displaystyle
f(k,\eta,L,Z_{n})=-k-\dfrac{2k\eta}{Z_{n}}-\dfrac{2k(L-\eta)}{Z_{n}^{2}}.$

Here $\eta$ is the Coulomb parameter, $k=|\sqrt{k^{2}}|$ is the wave number
expressed in fm-1 and determined by the energy of interacting particles in the
input channel, and $V_{n}=V(r_{n})$ is the interaction potential of clusters
at the point $r_{n}=nh$ from the interval of zero to $R$. The number of
equations $N$ or the dimension of the determinant, which usually turns out to
be in the range $100\ 000-1\ 000\ 000$ [45], $h=\Delta r/N$ is the step of the
finite difference grid and $\Delta r$ is the solution interval of the system
(usually from zero to $r_{N}=R$).

By writing $f(k,\eta,L,Z_{n})$ in the form given in Eq. (22) it is possible to
take the Coulomb interaction into account [33]. The form of the logarithmic
derivative of the WF in the external region can be obtained from the integral
representation of the Whittaker function [33]

$f(k,\eta,L,Z)=-k-\dfrac{2k\eta}{Z}-\dfrac{2k(L-\eta)}{Z^{2}}S(\eta,L,Z),$
(23)

where

$S(\eta,L,Z)=\dfrac{\int\limits_{0}^{\infty}t^{L+\eta+1}(1+t/Z)^{L-\eta-1}e^{-t}dt}{\int\limits_{0}^{\infty}t^{L+\eta}(1+t/Z)^{L-\eta}e^{-t}dt}.$
(24)

Calculations show that the value $S(\eta,L,Z)$ does not exceed 1.05, and its
effect on the binding energy of a two-particle system is negligible [45]. When
$f(k,0,0,Z)=-k$ in Eq. (23), the binding energy search process is noticeably
accelerated.

The calculation of the band determinant $D_{N}(k)$ for a given $k$ is carried
out using recurrent formulas of the form [47]

$\displaystyle D_{-1}=0,$ $\displaystyle
D_{n}=\thetaup_{n}D_{n-1}-\alphaup_{n}D_{n-2},$ (25) $\displaystyle D_{0}=1,$
$\displaystyle n=1,\dots,N.$

Any energy $E$ or wave number $k$ that leads to zero determinant

$D_{N}(k_{0})=0$ (26)

is an eigenenergy of the system $E_{b}$ or $k_{0}$, and the wave function at
this energy, determined by recurrent process below, is an eigenfunction of the
problem.

Methods for determining the zero of some functional of one variable $k$ are
well known [48]. The number $N_{D}$ of determinant values is determined
automatically from the accuracy condition of the binding energy value. The
latter one is usually set to the level $\epsilon\approx 10^{-5}$–$10^{-9}$MeV,
and $r_{N}=R$ is fixed on the range 20–30 fm [45].

After determining the eigenenergy $E_{b}$, the WF of this state is sought. To
find the shape of the eigenfunctions of bound states, the recurrent procedure

$\displaystyle u_{0}=0,$ $\displaystyle u_{n}=\thetaup_{n-1}u_{n-1}+u_{n-2},$
(27) $\displaystyle u_{1}=\text{const},$ $\displaystyle n=2,\dots,N.$

is carried out, where $u_{1}$ is an arbitrary number, usually fixed on the
range 0.01–0.1 [48].

For bound states, the determined WF is normalized to unity. Comparing it to
Whittaker asymptotics, one can find an asymptotic constant denoted by
$C_{\text{w}}$ (see Sec. IV).

The WF search area $R$ is usually of 20 to 30 fm, and the number of steps
$N_{WF}$ for the desired WF is fixed between 10 000 and 50 000. Only in the
case of a very low binding energy (0.1–0.2 MeV) the WF search area increased
to 100–200 fm or more.

The recurrence relation (27) is also used to search for WFs in the case of a
continuous spectrum of eigenvalues at predetermined positive energy
$(k^{2}>0)$ of interacting particles[45]. However, the WF must now be matched
with asymptotics of the form

$N_{L}u_{L}(r)\xrightarrow[r\rightarrow
R]{}F_{L}(kr)+\tan\left(\deltaup_{S,L}^{J}\right)G_{L}(kr).$ (28)

Matching the numerical solution $u_{L}(R)$ of Eq. (19) for two points at large
distances ($R$ on the order of 10–20 fm) with asymptotics (28), it is possible
to calculate the scattering phase shifts for each value of the momenta $JLS$
for a given energy of interacting particles, as well as the normalization of
the WF for scattering processes [45]. To calculate the WF, one can also use
the Numerov method [49]. When the number of steps exceeds 10 000, both methods
yield the same results within the typical required accuracy. Such results can
be compared by calculating the values of AC or charge radii for the BS or the
matrix elements for the scattering processes [45].

## Appendix B

Table 9: The astrophysical 6Li$(p,\gamma)^{7}$Be reaction rate in the range of 0.001 to 10 $T_{9}$ $T_{9}$ | Rate | $T_{9}$ | Rate | $T_{9}$ | Rate | $T_{9}$ | Rate
---|---|---|---|---|---|---|---
0.001 | $3.20\times 10^{-29}$ | 0.035 | $6.90\times 10^{-5}$ | 0.19 | $1.35\times 10^{0}$ | 2.25 | $7.80\times 10^{2}$
0.002 | $7.07\times 10^{-22}$ | 0.040 | $1.92\times 10^{-4}$ | 0.20 | $1.67\times 10^{0}$ | 2.50 | $9.07\times 10^{2}$
0.003 | $2.53\times 10^{-18}$ | 0.045 | $4.57\times 10^{-4}$ | 0.25 | $3.98\times 10^{0}$ | 2.75 | $1.03\times 10^{3}$
0.004 | $4.35\times 10^{-16}$ | 0.050 | $9.60\times 10^{-4}$ | 0.30 | $7.64\times 10^{0}$ | 3.00 | $1.16\times 10^{3}$
0.005 | $1.68\times 10^{-14}$ | 0.055 | $1.83\times 10^{-3}$ | 0.35 | $1.28\times 10^{1}$ | 3.25 | $1.29\times 10^{3}$
0.006 | $2.71\times 10^{-13}$ | 0.060 | $3.25\times 10^{-3}$ | 0.40 | $1.94\times 10^{1}$ | 3.50 | $1.42\times 10^{3}$
0.007 | $2.49\times 10^{-12}$ | 0.065 | $5.41\times 10^{-3}$ | 0.45 | $2.75\times 10^{1}$ | 3.75 | $1.56\times 10^{3}$
0.008 | $1.54\times 10^{-11}$ | 0.070 | $8.55\times 10^{-3}$ | 0.50 | $3.71\times 10^{1}$ | 4.0 | $1.69\times 10^{3}$
0.009 | $7.20\times 10^{-11}$ | 0.075 | $1.30\times 10^{-2}$ | 0.55 | $4.81\times 10^{1}$ | 4.5 | $1.95\times 10^{3}$
0.010 | $2.70\times 10^{-10}$ | 0.080 | $1.89\times 10^{-2}$ | 0.60 | $6.03\times 10^{1}$ | 5.0 | $2.22\times 10^{3}$
0.011 | $8.58\times 10^{-10}$ | 0.085 | $2.68\times 10^{-2}$ | 0.65 | $7.37\times 10^{1}$ | 5.5 | $2.50\times 10^{3}$
0.012 | $2.38\times 10^{-9}$ | 0.090 | $3.69\times 10^{-2}$ | 0.70 | $8.83\times 10^{1}$ | 6.0 | $2.77\times 10^{3}$
0.013 | $5.92\times 10^{-9}$ | 0.095 | $4.97\times 10^{-2}$ | 0.75 | $1.04\times 10^{2}$ | 6.5 | $3.05\times 10^{3}$
0.014 | $1.35\times 10^{-8}$ | 0.10 | $6.55\times 10^{-2}$ | 0.80 | $1.20\times 10^{2}$ | 7.0 | $3.32\times 10^{3}$
0.015 | $2.83\times 10^{-8}$ | 0.11 | $1.08\times 10^{-1}$ | 0.85 | $1.38\times 10^{2}$ | 7.5 | $3.59\times 10^{3}$
0.016 | $5.60\times 10^{-8}$ | 0.12 | $1.67\times 10^{-1}$ | 0.90 | $1.56\times 10^{2}$ | 8.0 | $3.86\times 10^{3}$
0.017 | $1.05\times 10^{-7}$ | 0.13 | $2.48\times 10^{-1}$ | 0.95 | $1.74\times 10^{2}$ | 8.5 | $4.13\times 10^{3}$
0.018 | $1.86\times 10^{-7}$ | 0.14 | $3.52\times 10^{-1}$ | 1.00 | $1.94\times 10^{2}$ | 9.0 | $4.38\times 10^{3}$
0.019 | $3.18\times 10^{-7}$ | 0.15 | $4.84\times 10^{-1}$ | 1.25 | $2.99\times 10^{2}$ | 9.5 | $4.63\times 10^{3}$
0.020 | $5.23\times 10^{-7}$ | 0.16 | $6.48\times 10^{-1}$ | 1.50 | $4.13\times 10^{2}$ | 10 | $4.88\times 10^{3}$
0.025 | $4.13\times 10^{-6}$ | 0.17 | $8.45\times 10^{-1}$ | 1.75 | $5.32\times 10^{2}$ |  |
0.030 | $1.98\times 10^{-5}$ | 0.18 | $1.08\times 10^{0}$ | 2.00 | $6.55\times 10^{2}$ |  |

## References

  * [1] C. A. Barnes, D. D. Clayton, and D. N. Schramm, eds., _Essays in nuclear astrophysics: presented to William A. Fowler, on the occasion of his seventieth birthday_ (Cambridge University Press, New York, 1982) p. 562.
  * [2] R. N. Boyd, C. R. Brune, G. M. Fuller, and C. J. Smith, Phys. Rev. D 82, 105005 (2010), arXiv:1008.0848.
  * [3] S. Bashkin and R. R. Carlson, Phys. Rev. 97, 1245 (1955).
  * [4] Z. E. Switkowski, J. C. Heggie, D. L. Kennedy, D. G. Sargood, F. C. Barker, and R. H. Spear, Nucl. Phys. A 331, 50 (1979).
  * [5] R. Ostojic, K. Subotic, and B. Stepancic, Il Nuovo Cimento A 76, 73 (1983).
  * [6] R. Bruss, in _Nuclei in the Cosmos: Proceedings of the Second International Symposium on Nuclear Astrophysics_ , edited by F. Kappeler and K. Wisshak (IOP Publishing Ltd, Karlsruhe, Germany, 1993) p. 648.
  * [7] T. Paradellis (1999), unpublished, quoted by [12] as Ref. [25].
  * [8] J. J. He, S. Z. Chen, C. E. Rolfs, S.W. Xu, J. Hu, X.W. Ma, M.Wiescher, R. J. DeBoer, T. Kajino, M. Kusakabe, L. Y. Zhang, S. Q. Hou, X. Q. Yu, N. T. Zhang, G. Lian, Y. H. Zhang, X. H. Zhou, H. S. Xu, G. Q. Xiao, and W. L. Zhan, Phys. Let. B 725, 287 (2013).
  * [9] D. Piatti, T. Chillery, R. Depalo, M. Aliotta, D. Bemmerer, A. Best, A. Boeltzig, C. Broggini, C. G. Bruno, A. Caciolli, F. Cavanna, G. F. Ciani, P. Corvisiero, L. Csedreki, T. Davinson, A. Di Leva, Z. Elekes, F. Ferraro, E. M. Fiore, A. Formicola, Z. Fulop, G. Gervino, A. Gnech, A. Guglielmetti, C. Gustavino, G. Gyurky, G. Imbriani, M. Junker, I. Kochanek, M. Lugaro, L. E. Marcucci, P. Marigo, E. Masha, R. Menegazzo, V. Mossa, F. R. Pantaleo, V. Paticchio, R. Perrino, P. Prati, L. Schiavulli, K. Stockel, O. Straniero, T. Szucs, M. P. Takacs, and S. Zavatarelli, Phys. Rev. C 102, 052802 (2020).
  * [10] R. M. Prior, M. C. Spraker, A. M. Amthor, K. J. Keeter, S. O. Nelson, A. Sabourov, K. Sabourov, A. Tonchev, M. Ahmed, J. H. Kelley, D. R. Tilley, H. R. Weller, and H. M. Hofmann, Phys. Rev. C 70, 10.1103 (2004).
  * [11] G. G. Kiss, M. La Cognata, R. Yarmukhamedov, K. I. Tursunmakhatov, I. Wiedenhover, L. T. Baby, S. Cherubini, A. Cvetinovic, G. D’Agata, P. Figuera, G. L. Guardo, M. Gulino, S. Hayakawa, I. Indelicato, L. Lamia, M. Lattuada, F. Mudo, S. Palmerini, R. G. Pizzone, G. G. Rapisarda, S. Romano, M. L. Sergi, R. Spart‘a, C. Spitaleri, O. Trippella, A. Tumino, M. Anastasiou, S. A. Kuvin, N. Rijal, B. Schmidt, S. B. Igamov, S. B. Sakuta, Z. Fulop, G. Gyurky, T. Szucs, Z. Halasz, E. Somorjai, Z. Hons, J. Mrazek, R. E. Tribble, and A. M. Mukhamedzhanov, Phys. Rev. C 104, 015807 (2021).
  * [12] K. Arai, D. Baye, and P. Descouvemont, Nucl. Phys. A 699, 963 (2002).
  * [13] J. T. Huang, C. A. Bertulani, and V. Guimaraes, Atomic Data and Nuclear Data Tables 96, 824 (2010), arXiv:0810.3867.
  * [14] S. B. Dubovichenko, N. Burtebaev, D. M. Zazulin, and A. S. Amar, Russian Physics Journal 53, 743 (2010).
  * [15] S. B. Dubovichenko, N. Burtebaev, D. M. Zazulin, Z. K. Kerimkulov, and A. S. Amar, Physics of Atomic Nuclei 74, 984 (2011).
  * [16] Y. Xu, K. Takahashi, S. Goriely, M. Arnould, M. Ohta, and H. Utsunomiya, Nucl. Phys. A 918, 61 (2013), arXiv:1310.7099.
  * [17] G. X. Dong, N. Michel, K. Fossez, M. Płoszajczak, Y. Jaganathen, and R. M. Betan, Journal of Physics G 44, 1 (2017), arXiv:arXiv:1601.06660v1.
  * [18] A. Gnech and L. E. Marcucci, Nucl. Phys. A 987, 1 (2019).
  * [19] S. B. Dubovichenko and A. V. Dzhazairov-Kakhramanov, Nucl. Phys. A 941, 335 (2015).
  * [20] P. Descouvemont, Frontiers in Astronomy and Space Sciences 7, 9 (2020).
  * [21] T. A. Tombrello and P. D. Parker, Phys. Rev. 131, 2582 (1963).
  * [22] K. Wildermuth and Y. C. Tang, _A Unified Theory of the Nucleus_ (Vieweg+Teubner Verlag, 1977).
  * [23] V. I. Kukulin, V. G. Neudatchin, I. T. Obukhovski, and Y. F. Smirnov, in _Clusters as Subsystems in Light Nuclei_ (Vieweg+Teubner Verlag, Wiesbaden, 1983) pp. 1–155.
  * [24] V. G. Neudatchin, V. I. Kukulin, V. N. Pomerantsev, and A. A. Sakharuk, Phys. Rev. C 45, 1512 (1992).
  * [25] C. Itzykson and M. Nauenberg, Reviews of Modern Physics 38, 95 (1966).
  * [26] A. Bohr and B. R. Mottelson, in _Nuclear Structure_ (World Scientific Publishing Company, 1998) pp. 137–307.
  * [27] A. M. Mukhamedzhanov and R. E. Tribble, Phys. Rev. C 59, 3418 (1999).
  * [28] S. B. Dubovichenko, Thermonuclear processes in Stars and Universe, 2nd ed. (Scholar’s Press, Saarbrucken, 2015) p. 332.
  * [29] Physical Measurement Laboratory; http://physics.nist.gov/cgi-bin/cuu/Value?mud%7Csearch_for =atomnuc!
  * [30] D. R. Tilley, C. M. Cheves, J. L. Godwin, G. M. Hale, H. M. Hofmann, J. H. Kelley, C. G. Sheu, and H. R. Weller, Energy levels of light nuclei A = 5, 6, 7 (2002).
  * [31] V. Varlamov, B. Ishkhanov, and S. Komarov, Atomic Nuclei Map (2015); http://cdfe.sinp.msu.ru/services/ground/NuclChart_release.html
  * [32] G. R. Plattner and R. D. Viollier, Nucl. Phys. A 365, 8 (1981).
  * [33] M. Abramowitz and I. A. Stegun, eds., Handbook of mathematical function, 10th ed. (U.S. Government Printing Office, Washington, 1972) p. 1046.
  * [34] Sukhoruchkin S.I., Supplement to I/25 A-F (Springer Berlin Heidelberg, 2016).
  * [35] K. M. Nollett and R. B. Wiringa, Phys. Rev. C 83, 10.1103 (2011).
  * [36] N. K. Timofeyuk, Phys. Rev. C 88, 044315 (2013).
  * [37] N. Burtebayev, J. Burtebayeva, N. Glushchenko, Z. Kerimkulov, A. Amar, M. Nassurlla, S. Sakuta, S. Artemov, S. Igamov, A. Karakhodzhaev, K. Rusek, and S. Kliczewski, Nucl. Phys. A 909, 20 (2013).
  * [38] M. Skill, R. Baumann, G. Keil, N. Kniest, E. Pfaff, M. Preiss, G. Reiter, G. Clausnitzer, M. Haller, and W. Kretschmer, Nucl. Phys. A 581, 93 (1995).
  * [39] C. Iliadis, Nuclear physics of stars, 2nd ed. (Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2015) p. 672.
  * [40] F. E. Cecil, D. Ferg, H. Liu, J. C. Scorby, J. A. McNeil, and P. D. Kunz, Nucl. Phys. A 539, 75 (1992).
  * [41] F. Barker, Australian Journal of Physics 33, 159 (1980).
  * [42] A. Amar and N. Burtebayev, Journal of Nuclear Sciences 1, 15 (2014).
  * [43] C. Angulo, M. Arnould, M. Rayet, P. Descouvemont, D. Baye, C. Leclercq-Willain, A. Coc, S. Barhoumi, P. Aguer, C. Rolfs, R. Kunz, J.W. Hammer, A. Mayer, T. Paradellis, S. Kossionides, C. Chronidou, K. Spyrou, S. Degl’Innocenti, G. Fiorentini, B. Ricci, S. Zavatarelli, C. Providencia, H. Wolters, J. Soares, C. Grama, J. Rahighi, A. Shotter, and M. Lamehi Rachti, Nucl. Phys. A 656, 3 (1999).
  * [44] G. R. Caughlan and W. A. Fowler, Thermonuclear reaction rates V (1988).
  * [45] S. B. Dubovichenko, Methods for calculating nuclear characteristics. Nuclear and thermonuclear processes., 2nd ed. (Lambert Academic Publishing, Saarbrucken, 2012) p. 425.
  * [46] S. B. Dubovichenko, A. V. Dzhazairov-Kakhramanov, and N. V. Afanasyeva, Nucl. Phys. A 963, 52 (2017).
  * [47] G. Marchuk and V. Kolesov, Application of Numerical Methods to Neutron Cross-Section Calculations (Atomizdat, Moscow, 1970) p. 304\.
  * [48] T. Korn and G. A. Korn, Mathematical Handbook for Scientists and Engineers (McGraw-Hill, New-York, 1968) p. 832.
  * [49] E. Hairer, S. P. Nørsett, and G. Wanner, Solving Ordinary Differential Equations I (Springer-Verlag Berlin Heidelberg, 1993) p. 528.


# Collective dynamics of evanescently coupled excitable lasers with saturable
absorber

M. Lamperti1 and A. M. Perego2
1,∗ Politecnico di Milano, Department of Physics and IFN-CNR, Via G. Previati
1/C, 23900, Lecco, Italy
2 Aston Institute of Photonic Technologies, Aston University, Birmingham B4
7ET, Aston Street, UK
<EMAIL_ADDRESS>
###### Abstract

We present a numerical study of the collective dynamics in a population of
coupled excitable lasers with saturable absorber. At variance with previous
studies where real-valued (lossy) coupling was considered, we focus here on
the purely imaginary coupling (evanescent wave coupling). We show that
evanescently coupled excitable lasers synchronize in a more efficient way
compared to the lossy coupled ones. Furthermore we show that out-of-diagonal
disorder-induced localization of excitability takes place for imaginary
coupling too, but it can be frustrated by nonvanishing linewidth enhancement
factor.

##  1 Introduction

The last decade has been characterized by the rise and spread of data storage
and analytics over many industrial and commercial fields, due to the
decreasing cost of data collection, storage and transmission. While it is
clear that large datasets are needed to obtain meaningful statistical
information about complex systems, the instruments used for such information
distillation have not evolved with the same speed as the data infrastructure,
due on one hand to our still rough understanding of large, complex systems,
and on the other hand to the exponentially increasing computing resources
needed to apply more and more complex models for data interpretation. A
promising framework for advanced data processing tasks is constituted by
neural networks [1], which currently see widespread usage to perform actions
that are traditionally a premise of human beings, such as pattern recognition,
natural language processing (information extraction from texts and
translation), and identification of objects in images. Although specialized
hardware and software has been developed to increment the efficiency of
implementing such framework with the currently most advanced computing
paradigma, based on a Von Neumann architecture executed on CMOS-based
integrated circuits, the level of complexity and efficiency of the human brain
is far from reach by many orders of magnitude.

The maturity and apparent limits of the current technology have pushed
scientists from many disciplines to propose new computing paradigms that could
enable a leap towards the higher processing power required for the
aforementioned purposes; the most promising results take inspiration from what
is considered to be the most advanced naturally evolved computer: the human
brain. The resulting field of neuromorphic computing aims at taking advantage
from systems that exhibit naturally the characteristics that make biological
neural networks so powerful and efficient, mapping the process paradigm to its
underlying dynamics rather than abstracting away from it, as it is currently
done with software implementation of neural networks on serialized, digital
hardware [2].

The two key ingredients for implementing neuromorphic computing are neuron-
like behavior and a large network of interconnections [1]. The first
ingredient is provided by excitable systems: excitability is an ubiquitous
process in nature mostly known in biology [3, 4] in the context of the
neuronal cell activity and can be defined as the generation of spike-like
behaviour in one or more system dynamical variables in response to external
perturbations whose magnitude exceeds a given threshold. In such regime the
system response does not depend on the perturbation strength and each spike
generation is followed by a refractory time during which the system remains
silent; after such refractory time the emission of another spike can take
place again. Excitability needs to be provided by a suitable neuron-like
system, which in turn will constitute the building block to perform simple
computations in a highly parallel fashion, making the system fast and
efficient. What is needed to perform complex computations is a large network
of interconnects that weight the outputs of previous neurons and route them to
other computation elements, i.e. other neuron-like blocks. A network of
interconnects is thus the key ingredient to transform a bunch of excitable
systems into a highly-efficient computer able to solve complex tasks.

Currently, the golden standard for electronics-based computing, CMOS
integrated circuits, is unable to provide these two elements while keeping the
power efficiency high [5]. As we approach the few-atom transistor, with bigger
and bigger challenges without any significant decrease of power consumption,
new directions have been probed in search for another platform that could
provide high integration, low power consumption and scalability. Thanks to
recent advances in optoelectronics integration, photonics technologies
constitute one of the most promising platforms for neuromorphic computing:
neuromorphic photonics is gaining momentum as a research field where emerging
photonics technologies are used to mimic neuronal dynamics and/or to perform
computational tasks based on brain inspired strategies [6]. Excitability in
photonics has been reported in a variety of different lasers and amplifiers
systems [7, 8, 9, 10, 11, 12, 13, 14, 15]. Possibly the first historically
studied example of a photonics system exhibiting excitability is the
semiconductor laser with saturable absorber [7]. Excitability in such laser is
enabled by a phase space portrait exhibiting a limit cycle close to a saddle
node bifurcation. In the regime where the absorption is the slow dynamical
variable, if the below, but close-to-threshold laser is perturbed strongly
enough, then the stimulated emission process builds up producing the emission
of giant light pulse. Such emission depletes the gain and is followed by a
refractory time after which the laser is ready to be excited again. Such
dynamics can be associated to a so called type III excitability [8].

Engineering optically excitable elements connectivities and studying their
collective properties and dynamics are both crucial tasks towards achieving a
general understanding of neuromorphic photonics systems. Focussing our
attention to the excitable laser with saturable absorber, it is important to
stress that coupling engineering of excitable semiconductor lasers with
saturable absorber has been shown to allow pattern recognition [14], neuronal
circuits design [16] and coincidence detection devices [17]. As far as the
collective dynamics is concerned, we have shown theoretically in two recent
works [18, 19] that temporal and intensity synchronization, array enhanced
coherence resonance, and even disorder-induced localization of excitability
can take place in arrays of excitable lasers with nearest neighbour real-
valued (lossy) coupling. In this work we extend our previous studies on
synchronization and disorder-induced localization of excitability to the case
of an ensemble of excitable lasers with saturable absorber coupled via a
nearest neighbour purely imaginary coupling coefficient. Imaginary valued
couplings describe the physical evanescent waves interaction between adjacent
laser cavities. Such interaction is most likely the relevant coupling
mechanism for micropillars lasers, where collective excitable dynamics could
be observed [20].

##  2 The model

We consider here a population of $n$ lasers with nearest neighbour coupling.
The following normalized Yamada model describes the $i$-th laser dynamics:

$\displaystyle\dot{F}_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{2}[G_{i}(1-i\alpha)-Q_{i}(1-i\beta)-1]F_{i}+\sigma_{i}$
$\displaystyle-$ $\displaystyle
i(K_{i,i+1}+K_{i,i-1})F_{i}+iK_{i+1,i}F_{i+1}+iK_{i-1,i}F_{i-1},$
$\displaystyle\dot{G}_{i}$ $\displaystyle=$
$\displaystyle\gamma(A-G_{i}-I_{i}G_{i}),$ $\displaystyle\dot{Q}_{i}$
$\displaystyle=$ $\displaystyle\gamma(B-Q_{i}-aQ_{i}I_{i}).$ (1)

Here $F_{i}$ denotes the electric field strength, $I_{i}=|F_{i}|^{2}$ its
intensity, $G_{i}$ and $Q_{i}$ gain and absorption respectively. $A$ is the
pump parameter, $B$ the background absorption, $a$ the differential absorption
relative to the differential gain, $\gamma$ is the absorber and gain decay
rate, $\alpha$ and $\beta$ denote the linewidth enhancement factors for the
gain and the absorber respectively; these parameters do not have subscript $i$
since they have been considered identical for all lasers. $\sigma_{i}$
describes a delta correlated Gaussian noise term of strength $D$ with
$\langle\sigma_{i}(t_{1})\sigma_{j}(t_{2})\rangle=\sqrt{2D}\delta(t_{1}-t_{2})\delta_{ij}$
providing the perturbations needed for excitable behaviour. The dot denotes
temporal derivative and the time variable has been normalized to the uncoupled
laser photon lifetime. $K_{i,j}$ denotes the nearest neighbour coupling
strength describing light coupling from the $i$-th to the $j$-th laser where
the reciprocity condition $K_{i,j}=K_{j,i}$ has been imposed. The identical
coupling case $K_{i,i+1}=K$ results in an effective discrete Laplace
diffraction operator in the array.
Periodic conditions at the array boundary have been applied. Across the whole
paper we have set $A=6.5$, $B=5.8$, $a=1.8$, $\gamma=10^{-3}$ while for
$\alpha$ and $\beta$ different values have been used and specified across the
paper.

##  3 Synchronization

We have first studied the effect of coupling on the synchronization of the
firing events in the lasers array as a function of the number of lasers
constituting the array itself and of the values of the linewidth enhancement
factors. To this purpose we have considered independent additive noise sources
to be present in all lasers and varied their amplitude, $D$, from 0 to 0.15.
The coupling strength, $K=K_{0}$, has been considered identical for all lasers
and has been varied from 0 to 1. For each pair $(D,K)$ we run one simulation
for a time duration $T=100000$. The pulse synchronization can be understood
both in space and time: in the first case it refers to neighbor lasers
emitting a pulse at the same time while in the second case it characterizes a
regular (i.e. equally spaced) emission of pulses from the single laser. To
give a qualitative flavour of temporally synchronous pulses (spikes) emission
of coupled lasers we have plotted in Fig. 1 the temporal traces corresponding
to coupled and uncoupled laser (left and right column, respectively).

Figure 1: In panels a) to e) the temporal traces of five different uncoupled
lasers have been plotted showing uncorrelated lasers firing. In panels f) to
j) the time traces of five randomly-picked up lasers in a chain of 50 coupled
elements have been depicted. Spatial synchronization is clearly evident in
presence of non vanishing coupling. Parameters used are $D=0.01$,
$\alpha=\beta=0$ for both the coupled and uncoupled case; $K_{0}=0$ in a) to
e) and $K=0.05$ in f) to j).

We have characterized the spatial synchronization by using an indicator ($S$)
that quantifies the phase slip between the firing events of nearest-neighbors
lasers [18, 21]. The phase of the $i$-th laser reads:

$\displaystyle\phi_{i}(t)=\frac{t-\tau_{k}}{\tau_{k+1}-\tau_{k}}+2k\pi$ (2)

where $\tau_{k}$ is the time of the $k$-th firing event, i.e. the position in
time of the $k$-th emitted pulse. We define then

$\displaystyle s_{i}=\sin\left(\frac{\phi_{i}-\phi_{i+1}}{2}\right)^{2}$ (3)

which characterizes the phase synchronization between the $i$-th and $i+1$-th
lasers. The average both across all the array elements and along temporal
duration of the time trace gives the $S$ indicator, which provides a measure
of the spatial synchronization of the whole system:

$\displaystyle
S=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}\left(\frac{1}{n}\sum_{i=1}^{n}s_{i}\right)dt.$
(4)

The maximum synchronization occurs for $S=0$ while for the completely non
synchronized state $S=0.5$.

We have furthermore characterized the regular emission of pulses of the
individual laser (temporal synchronization) by means of the single laser pulse
jitter defined as $J=\sigma_{T}/\langle T\rangle$, where $\langle T\rangle$ is
the average time interval between two consecutive pulses and $\sigma_{T}$ is
the standard deviation. A value of $J$ close to 0 indicates regular pulses
emission, while values close or greater than 1 indicate poor regularity.

Figure 2: The single-laser temporal synchronization index, $J$, is plotted in
the $(D,K)$ space for different numbers of lasers (rows - 4, 20 and 150 from
top to bottom) and different values of the line enhancement factors (columns -
0, 1, 3 and 5 from left to right).

The panels in Fig. 2 and in Fig. 3 show the $J$ and $S$ synchronization
indicators, respectively, for different number of lasers constituting the
chain and different values of the line enhancement factors, versus the noise
and the coupling strength. The behavior observed here is qualitatively
different from the case of real (dissipative) coupling, showing no array-
enhanced coherence resonance [18]. Instead of enlarging the region of the
$(D,K)$ parameter space characterized by high synchronization, an increase in
the number of coupled lasers improves the temporal synchronization in a
limited region of the parameters space (the lower-left one) while worsening
the synchronization in all the remaining region. In other words, enlarging the
array creates a sharper separation between the regions of high and low
synchronization, without changing its area.
On the other hand, increasing the number of lasers in the array has the effect
of reducing the region of highest spatial synchronization, i.e. low value of
the $S$ indicator. This effect is well visible in the part of the $(D,K)$
space characterized by $D>0.05$, and manifests itself quickly with the
increase of lasers forming the chain (compare first and second row of Fig. 3).
It is worth noting that even in the worst synchronization cases, the $S$
indicator is smaller than $10^{-5}$, indicating good spatial synchronization
in all the $(D,K)$ space.

Figure 3: The logarithm of the spatial synchronization, $\log_{10}(S)$, is
plotted in the $(D,K)$ space for different numbers of lasers (rows - 4, 20 and
150 from top to bottom) and different values of the line enhancement factors
(columns - 0, 1, 3 and 5 from left to right).

The second control parameter we considered, the line enhancement factor,
differently from the number of lasers modifies the shape and size of the
synchronization region. For a small number of lasers (4, first row in Figs.
2-3), its progressive increase from 0 to 5 causes an enlargement of the region
of high synchronization towards larger values of the noise, both in time and
space. For a longer chain (20 and 150 lasers, two bottom rows), high
synchronization occurs in a more and more localized region towards small
values of $K$ and $D$; even in a strong coupling regime high temporal
synchronization of spikes emitted by the same laser results more difficult to
achieve.

As compared to the real coupling case [18] we can observe that a much smaller
value of the coupling strength is needed to achieve comparable
synchronization. Imaginary coupling is hence a much more efficient way to
achieve synchronization of excitability.

##  4 Disorder-induced localization

In our previous study for real-valued couplings we pointed out how randomness
in the laser array coupling strength is strikingly able to induce spatial
exponential localization of excitable behavior in a given area of the array
[19]. Such disorder-induced localization does not occur due to a trivial
breaking of the lasers chain but due to a dynamical process entailing both
scattering and dissipation. We want to stress that the localization we are
discussing here, although mediated by disorder, can not be assimilated
directly to the celebrated _Anderson localization_ [22]. The latter, although
it has inspired our work, is occurring in ideally conservative systems while
in our system dissipation plays of course a paramount important role. Still
inspired by the conspicuous literature on _Anderson localization_ we can
borrow some useful terminology. In the traditional studies of _Anderson
localization_ of electronic transport in solid state physics, the distinction
between diagonal and out-of-diagonal disorder is made when referring to
disorder on the individual atomic sites, or in the hopping terms between
neighbour atomic sites respectively [23]. In our system we identify clearly
the former as corresponding to randomness in some internal laser parameter
(e.g. pump parameter), while the latter to the coupling between the lasers. In
this paper we will consider out-of-diagonal disorder defining the coupling as:
$K_{i,i+1}=K_{0}+\rho_{i,i+1}$, being $K_{0}$ a value common to all lasers and
$\rho_{i,i+1}$ a random number constant in time drawn from a uniform
distribution in the interval $[-r,+r]$ which describes light coupling from the
$i$-th laser to nearest neighbour on the right. We recall also that the
reciprocity constraint implies $\rho_{i,i+1}=\rho_{i+1,i}$.

Figure 4: In a) the intensity dynamics in absence disorder: an excitability
wave propagates through the array. In b) an example of the temporal dynamics
in presence of disorder is shown: excitability is spatially localized.
Parameters used are: $\alpha=\beta=0$, $K_{0}=0.05$, $r=0.005$.

In order to investigate the disorder-induced localization we have considered
that the additive noise is present, without loss of generality, only in the
central laser of the array. If all the lasers are coupled with identical
coupling strength ($K_{i,i+1}=K_{0}$ $\forall$ $i$), then an excitability wave
propagates from the centre of the array both towards left and right (see Fig.
4). We have considered a chain of 150 lasers with an average coupling
$K_{0}=0.05$, a noise on the central laser with strength $D=0.1$ and a random
coupling distribution varying in the interval $r\in[0,0.002]$. We have
verified that in this condition, with noise acting only on one laser, we have
regular emission of pulses that propagate to all the lasers in the chain for
$r=0$.

Figure 5: The temporal average of 100 traces characterized by different
realizations of disorder with the same value of $r$ is plotted here together
with the fit to Eq. 5 used to extract the value of the localization exponent,
$\lambda$. To extract the correct localization exponent we discard the central
laser, whose average intensity is systematically too high due to the presence
of additive noise, and its nearest neighbors. Parameters used are the same of
Fig. 6, with $\alpha=\beta=0$ and $r=0.002$.

For each value of $r$, 100 numerical experiments with a time duration of
$T=50000$ has been carried out. For each disordered coupling configuration we
have computed the average intensity trace across the array. We have then
fitted the averaged intensity distribution after the 100 realizations of the
disorder with the same strength $r$, with the following exponential function:

$\displaystyle f=b+\exp\left(-\lambda|i-i_{0}|\right)$ (5)

where $i_{0}$ denotes the position of the central laser (see Fig. 5 for an
example). This procedure has been repeated 8 times, so that the average value
and standard deviation of the fit parameters can be extracted; the average
localization exponent $\langle\lambda\rangle$ versus the disorder strength $r$
is plotted in Fig. 6, together with error bars indicating the standard
deviation.

Figure 6: The average localization exponent $\langle\lambda\rangle$ has been
plotted versus randomness strength $r$ for 3 different values of the linewidth
enhancement factors $\alpha$ and $\beta$. For vanishing $\alpha$ and $\beta$ a
transition from ballistic to localized regime occurs, low values of $\alpha$
and $\beta$ require an higher threshold for localization. At large values of
the linewidth enhancement factor localization is lost. Parameters used are
$D=0.1$ and $K=0.05$.

Such disorder-induced localization occurs for coupling values such that if we
couple all the lasers with identical value of $K_{0}-r$ equal to the minimum
possible value generated by the random process, we would still have an
excitability wave spreading from the array center both towards left and right;
this ensures the fact that localization is not caused by a trivial breaking of
the laser chain but is indeed provoked by a dynamical effect. For vanishing
$\alpha$ and $\beta$ we can notice that localization occurs for coupling
strength much smaller than the ones needed in the real coupling case (about
0.3, see [19]).
We have furthermore investigated the impact of linewidth enhancement factor
both for the gain medium ($\alpha$) and for the saturable absorber ($\beta$)
on the disorder-induced localization of excitability. Physically $\alpha$ and
$\beta$ are responsible for coupling of amplitude to phase fluctuations. It is
interesting to appreciate that linewidth enhancement factor
($\alpha$,$\beta\neq 0$) causes a decrease (for $\alpha=\beta=1$) and even a
vanishing (for $\alpha=\beta=3$ or $5$) of the localization.

##  5 Conclusions

We have demonstrated with the help of numerical simulations that evanescently
coupled excitable lasers with saturable absorber synchronize in a more
efficient way compared to lossy coupled ones, but that they do not exhibit
array enhanced coherence resonance. We have furthermore demonstrated that
disorder-induced localization of excitability caused by randomness in the
coupling strength exists for imaginary coupling too. Taking into account the
linewidth enhancement factor both in the gain medium and in the saturable
absorber we have shown that localization is reduced and eventually quenched by
an increase of the linewidth enhancement factor. This fact may constitute a
serious obstacle in the experimental observation of the disorder-induced
localization of excitability. Our results shed further light on the collective
dynamics of coupled excitable lasers with saturable absorber and suggest
interesting directions for experimental studies.

##  6 Acknowledgements

The authors gratefully acknowledge useful discussion with Dr. Sylvain Barabay.

## References

  * [1] S. Furber, Journal of neural engineering 13, 051001 (2016)
  * [2] B.J. Shastri, A.N. Tait, T.F. de Lima, M.A. Nahmias, H.T. Peng, P.R. Prucnal, _Principles of neuromorphic photonics_ (2017), 1801.00016
  * [3] J.D. Murray, _Mathematical Biology: I. An Introduction_ , Vol. 1 (Springer-Verlag, 3rd Edition Berlin Heidelberg, 2002)
  * [4] E.M. Izhikevich, International Journal of Bifurcation and Chaos 10, 1171 (2000)
  * [5] I.K. Schuller, R. Stevens, Tech. rep. (2015), https://science.energy.gov/$\sim$/media/bes/ pdf/reports/2016/NCFMtSA_rpt.pdf
  * [6] P.R. Prucnal, B.V. Shastri, _Neuromorphic Photonics_ (Taylor & Francis, 1st Edition Boca Raton, 2017)
  * [7] J.L. Dubbeldam, B. Krauskopf, D. Lenstra, Physical Review E 60, 6580 (1999)
  * [8] S. Barbay, R. Kuszelewicz, A.M. Yacomotti, Optics letters 36, 4476 (2011)
  * [9] S. Barland, O. Piro, M. Giudici, J.R. Tredicce, S. Balle, Physical Review E 68, 036209 (2003)
  * [10] M. Giudici, C. Green, G. Giacomelli, U. Nespolo, J. Tredicce, Physical Review E 55, 6414 (1997)
  * [11] F. Plaza, M. Velarde, F. Arecchi, S. Boccaletti, M. Ciofini, R. Meucci, EPL (Europhysics Letters) 38, 85 (1997)
  * [12] M. Brunstein, A.M. Yacomotti, I. Sagnes, F. Raineri, L. Bigot, A. Levenson, Physical Review A 85, 031803 (2012)
  * [13] M. Turconi, M. Giudici, S. Barland, Physical Review Letters 111, 233901 (2013)
  * [14] B.J. Shastri, M.A. Nahmias, A.N. Tait, A.W. Rodriguez, B. Wu, P.R. Prucnal, Scientific Reports 6, 19126 (2016)
  * [15] B. Romeira, R. Avó, J.M. Figueiredo, S. Barland, J. Javaloyes, Scientific reports 6, 19510 (2016)
  * [16] B.J. Shastri, M.A. Nahmias, A.N. Tait, B. Wu, P.R. Prucnal, Optics express 23, 8029 (2015)
  * [17] B.J. Shastri, A.N. Tait, M. Nahmias, B. Wu, P. Prucnal, _Coincidence detection with graphene excitable laser_ , in _CLEO: Science and Innovations_ (Optical Society of America, 2014), pp. STu3I–5
  * [18] A.M. Perego, M. Lamperti, Physical Review A 94, 033839 (2016)
  * [19] M. Lamperti, A.M. Perego, Physical Review A 96, 041803(R) (2017)
  * [20] F. Selmi, R. Braive, G. Beaudoin, I. Sagnes, R. Kuszelewicz, S. Barbay, Opt. Lett. 40, 5690 (2015)
  * [21] M.G. Rosenblum, A.S. Pikovsky, J. Kurths, Circuits and Systems I: Fundamental Theory and Applications, IEEE Transactions on 44, 874 (1997)
  * [22] P.W. Anderson, Physical Review 109, 1492 (1958)
  * [23] J.B. Pendry, J. Phys. C: Solid State Phys. 15, 5773 (1982)